/examples/example_V2.5.R

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