/examples/example_V2.5.R

http://github.com/sbotond/phylosim · R · 60 lines · 17 code · 15 blank · 28 comment · 0 complexity · 45e21a8dd0399b09922f26f467cb43d4 MD5 · raw file

  1. #!/usr/bin/env Rscript
  2. ##
  3. ## Example V2.5: Simulating partitions
  4. ## See also the package vignette (vignette("PhyloSim",package="phylosim")).
  5. ##
  6. # load PhyloSim
  7. library("phylosim")
  8. # The following example demonstrates how to use the processes and site- and process-specific
  9. # parameters to simulate ``partitions'' with different properties.
  10. # We will simulate four partitions:
  11. # * Partition 1: sites in range \code{1:25} evolving by JC+d$\Gamma$ with a shape parameter
  12. # alpha=1
  13. #
  14. # * Partition 2: sites in range \code{26:50} evolving by JC+d$\Gamma$ with a shape parameter
  15. # alpha=0.5
  16. #
  17. # * Partition 3: sites in range \code{51:75} evolving by HKY+d$\Gamma$ with a shape parameter
  18. # alpha=1
  19. #
  20. # * Partition 4: sites in range \code{76:100} evolving by HKY+d$\Gamma$ with a shape parameter
  21. # alpha=0.5
  22. # First construct two substitution process objects:
  23. jc69<-JC69()
  24. hky<-HKY(rate.params=list( "Alpha"=5,"Beta"=2),
  25. base.freqs=c(4,3,2,1)/10
  26. )
  27. # Construct a root sequence object of length 100:
  28. root.seq<-NucleotideSequence(length=100)
  29. # Attach process jc69 to range 1:50:
  30. attachProcess(root.seq,jc69,1:50)
  31. # Attach process hky to range 51:100:
  32. attachProcess(root.seq,hky,51:100)
  33. # Sample rate multipliers in the four partitions:
  34. plusGamma(root.seq,jc69,1,1:25)
  35. plusGamma(root.seq,jc69,0.5,26:50)
  36. plusGamma(root.seq,hky,1,51:75)
  37. plusGamma(root.seq,hky,0.5,76:100)
  38. # Construct the PhyloSim object, sample states, set root sequence, set the phylo object (random
  39. # coalescent tree for three taxa) and run the simulation:
  40. sim<-Simulate(PhyloSim(
  41. root.seq=sampleStates(root.seq),
  42. phylo=rcoal(3)
  43. ))
  44. # Plot the alignment alongside the tree, skip sequences at ancestral nodes:
  45. plot(sim,plot.ancestors=FALSE)