/examples/example_V2.2.R

 1#!/usr/bin/env Rscript
 2
 3##
 4## Example V2.1.1: Simulating substitutions under the HKY model
 5## See also the package vignette (vignette("PhyloSim",package="phylosim")).
 6##
 7
 8# load PhyloSim
 9library("phylosim")
10
11# Construct an HKY substitution process object:
12p<-HKY(rate.params=list( "Alpha"=10,"Beta"=2),
13                     base.freqs=c(4,3,2,1)/10
14   )
15
16# Get a plot of the instantaneous substitution matrix and equilibrium distribution of the p process (bubble plot):
17plot(p,scale=0.5)
18
19# Construct the root sequence, attach the substitution
20# process to root sequence via the constructor and sample states:
21root.len50.seq<-NucleotideSequence(length=50,processes=list(list(p)))
22sampleStates(root.len50.seq)
23
24# Print out root sequence:
25print(root.len50.seq)
26
27# Construct a PhyloSim object, set the phylo object and the root sequence:
28sim<-PhyloSim(
29        root.seq=root.len50.seq,
30        phylo=read.tree("data/3taxa.nwk")
31);
32
33# Run simulation:
34Simulate(sim)
35
36# Save the resulting alignment (fasta format):
37saveAlignment(sim,file="HKY_sim.fas")
38
39