/examples/example_V2.7.R

 1#!/usr/bin/env Rscript
 2
 3##
 4## Example V2.7: Simulating many replicates in parallel
 5## See also the package vignette (vignette("PhyloSim",package="phylosim")).
 6##
 7
 8# load PhyloSim
 9library("phylosim")
10
11# The the speed of the above method for simulating replicates can be improved on a multicore machine by running
12# many replicates in parallel by using the mclapply method from the multicore package. 
13
14# Under default settings, the mclapply method launches one replication per core and this approach needs enough memory to run all
15# of them in parallel.
16
17# The following code illustrates how to simulate many replicates in parallel under the JC69+dG model.
18
19# Construct the root sequence object and attach the substitution process:
20p<-JC69();
21root.seq<-NucleotideSequence(length=50)
22attachProcess(root.seq,p)
23
24# Read the required phylogeny from file (this will remain fixed in the simulated replicates):
25tree<-read.tree("data/3taxa.nwk");
26
27# Function to simulate a single replication: 
28sim.replicate<-function(i){
29        name<-paste("replication_",i,sep="")
30        clearStates(root.seq)
31        plusGamma(root.seq,p,0.25)
32        sampleStates(root.seq)
33
34        sim<-Simulate(PhyloSim(
35                name=name,
36                root.seq=root.seq,
37                phylo=tree,
38        ),
39            quiet=TRUE
40        )
41
42        saveAlignment(sim,file=paste("aln_",i,".fas",sep=""))
43        return(sim)
44        # return(TRUE)
45}
46
47# Note that the states are cleared and resampled; the rate multipliers 
48# are resampled as well. The resulting alignments are stored in files aln_1.fas, aln_2.fas, aln_3.fas.
49
50# Memory can be saved by throwing away the objects generated by the replication by returning TRUE from the
51# sim.replicate function.
52
53# Load the multicore package:
54library(multicore)
55
56# Run replicates in parallel:
57nr.replicates <-3
58res.objects<-mclapply(1:nr.replicates, sim.replicate)
59
60# Print the resulting PhyloSim objects:
61print(res.objects)
62
63# Plot one of the resulting alignments:
64plot(res.objects[[1]])
65