Skip to content
This repository
tree: 1c51cf90e4
Fetching contributors…

Octocat-spinner-32-eaf2f5

Cannot retrieve contributors at this time

executable file 104 lines (75 sloc) 2.283 kb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
#!/usr/bin/env Rscript

##
## Example V3.2: Evolving codon sequences
## See also the package vignette (vignette("PhyloSim",package="phylosim")).
##

library(phylosim)

# Enable "fast & careless" mode:
PSIM_FAST<-TRUE;

# Construct a GY94 codon substitution model:
p<-GY94();

# Set the transition/transverion rate ratio:
p$kappa=2

# Sample codon frequencies from a normal distribution:
p$equDist<-abs(rnorm(61,mean=10,sd=3))

# Get object summary for p:
summary(p)

# Get a bubble plot of p:
plot(p,scale=0.8)

# Construct a discrete deletor process:
d<-DiscreteDeletor(
rate=1,
sizes=1:4,
probs=c(4,3,2,1)/10
);

# Construct a discrete insertor process inserting neutrally evolving sites:
i<-DiscreteInsertor(
rate=1.5,
sizes=1:4,
probs=c(4,3,2,1)/10,
template.seq=CodonSequence(length=4,processes=list(list(p)))
);

# Construct root sequence and attach process p:

s<-CodonSequence(length=30,processes=list(list(p)))

# Sample omegas from a discrete model:
omegaVarM3(s,p,omegas=c(0,1,2),probs=c(2/4,1/4,1/4))

# Plot the omega values across sites:
plotParametersAtSites(s,p,"omega");

# Sample states:

sampleStates(s)

# Construct the simulation object:
sim<-PhyloSim(
root.seq=s,
phylo=read.tree("data/smalldemotree.nwk")
);

# Create a node hook function and attach to node 9:
node.hook<-function(seq){

# Set all omegas to 1 (neutral):
setOmegas(seq,p,1);
# attach the deletion process:
attachProcess(seq,d)
# attach the insertion process:
attachProcess(seq,i)

        return(seq);
}

attachHookToNode(
                sim, # PhyloSim object.
                node=9, # the node
                fun=node.hook # the node hook function
);

# Disable fast mode just before simulation in order to preserve branch statistics:
rm(PSIM_FAST)

# Run the simulation:
Simulate(sim)

# Plot the resulting alingment alongside the tree:
plot(sim);

# Export the nonsynonymous substitution counts as a phylo object:
nsyn.subst<-exportStatTree(sim,"nr.nsyn.subst")

# Plot the exported phylo object:
plot(nsyn.subst)
nodelabels()

# Save the resulting alignment:
saveAlignment(
                sim, # the phylo object
                file="example_V3.2_aln.fas", # filename for alignment
);

Something went wrong with that request. Please try again.