Hello again everyone,
forwarded is some code provided by Dave Bapst for simulating a certain
number of trait shifts in nodes AND tracking those changes as you go.
This code will simulate "Nchanges" shifts and give both tip (tipChar)
and node (intNodeChar) states at the end of the run.
Dave correctly suggested this might be useful for other people, so I am
forwarding his code with my apologies for discussing this with him offline.
Dave, thanks again for the help. I think this can be really helpful for
some simulation purposes.
Cheers,
Antigoni
On 10/13/2011 02:53 PM, David Bapst wrote:
Hello Antigoni,
I think most people might actually want to know that, but I didn't
code it to make it easy to figure it out. So, I can understand why you
had trouble getting my simple brutish kludge to do it. The following
code should do what you are looking for. If you find this helpful,
please forward this email to R-Sig-Phylo for posterity.
-Dave
#a different version which will track the character changes
library(phangorn)
tree<-rtree(30) #you can use your own tree here or rcoal() or whatever
Nchanges<-10 #number of shifts at nodes
nodeShifts<-sample((2:Nnode(tree))+Ntip(tree),Nchanges)
DescList<-Descendants(tree,nodeShifts,type="all")
nodeChar<-rep(T,Ntip(tree)+Nnode(tree))
for(i in 1:length(DescList)){
nodeChar[DescList[[i]]]<-!nodeChar[DescList[[i]]]
nodeChar[nodeShifts[i]]<-!nodeChar[nodeShifts[i]]
}
nodeChar<-as.numeric(nodeChar)
tipChar<-nodeChar[1:Ntip(tree)]
intNodeChar<-nodeChar[-(1:Ntip(tree))]
plot(tree,show.tip.label=F)
tiplabels(pch=21,bg=tipChar,cex=2)
nodelabels(pch=21,bg=intNodeChar,cex=2)
nodeShifts #which nodes the changes occurred at
tipChar #characters at tips
nodeChar #character at each node, including tips
intNodeChar #character at each internal node
On Thu, Oct 13, 2011 at 1:50 PM, Antigoni Kaliontzopoulou
<[email protected]> wrote:
Hi Dave,
another question (off-R-sig-phylo because I don't think it's particularly
interesting for people): do you by any chance also have any code that will
track the final node character states after simulating the changes?
I've been banging my head on that for some hours now, so I'd really
appreciate any help with that...
Thanks a lot!
Cheers,
Antigoni
On 10/13/2011 07:02 AM, David Bapst wrote:
Hello Rob and Antigoni,
Would the following meet your purposes?
Cheers,
-Dave Bapst, UChicago
library(ape)
tree<-rtree(30) #you can use your own tree here or rcoal() or whatever
Nchanges<-10 #number of shifts at nodes
nodeShifts<-sample((2:Nnode(tree))+Ntip(tree),Nchanges)
DescList<-prop.part(tree)[nodeShifts-Ntip(tree)]
tipChar<-rep(T,Ntip(tree))
for(i in 1:length(DescList)){
tipChar[DescList[[i]]]<-!tipChar[DescList[[i]]]
}
tipChar<-as.numeric(tipChar)
plot(tree,show.tip.label=F)
tiplabels(pch=21,bg=tipChar,cex=2)
nodeShifts #which nodes the changes occurred at
tipChar #characters at tips
On Wed, Oct 12, 2011 at 10:33 PM, Rob Lanfear<[email protected]>
wrote:
Hi Antigoni,
I had been thinking about some of these kinds of simulations recently
too.
One challenge I encountered (which is similar to your problem here, if
I've
understood it) was how to make my simulations actually relevant.
Specifically, given a particular tree I needed to simulate data with N
changes in a trait on the phylogeny, and I needed to make sure that what
I
simulated was something I would have recovered using the simulated data
at
the tips. This won't always be the case for Liam's method (if I've
understood it correctly), because even if I fix the tree (trivial) I
could
end up with simulated data (e.g. multiple changes on a single branch)
which
I would not then reconstruct using the data from the tips. Note that this
isn't a problem with Liam's method, but would be a problem if I wanted to
use it for my purposes.
I gave up in my attempt to solve this, but this thread, and Liam's
approach,
made me wonder if this problem could be solved using a related approach.
In
particular, one could use a fixed tree, and change the loop in Liam's
method
to randomise the tip labels, re-do the ancestral state reconstruction,
and
keep going until you get a reconstruction with the desired number of
changes. That would make sure that you only simulate datasets which you
would have recovered with extant data, and it also has the advantage (for
me
at least) of holding the topology and the number of extant species with
each
trait value constant across replicates.
I'd be interested to know if anyone thinks this is a reasonable approach?
Cheers,
Rob
On 13 October 2011 02:20, Antigoni
Kaliontzopoulou<[email protected]
wrote:
Hi Liam,
thank you very much for this. Yes, the function works beautifully for
what
we wanted to do and the code you sent greatly facilitates
implemantation.
I was, however, wondering whether you can somehow "limit" the shifts to
happen in nodes, to avoid character reversals across a single branch.
Again, thanks a lot for the help.
Cheers,
Antigoni
On 10/11/2011 08:34 PM, Liam J. Revell wrote:
Hi Antigoni.
I have a function in the new version of "phytools" called sim.history()
that might help with this. This version is not available yet on CRAN,
but
it can be downloaded from my webpage: http://anolis.oeb.harvard.edu/**
~liam/R-phylogenetics/<http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/>.
The function generates stochastic histories according to some rate
matrix,
and stores the histories along the states at all terminal nodes.
This code might be helpful in your task. It simulates random trees,
for a
given number of changes it computes the instantaneous rate matrix that
will
produce that number of changes, and then it simulates (repeatedly)
under
that rate until the desired number of changes is obtained. The code
fragment uses functions in "ape", "geiger", and "phytools":
nchanges<-5 # desired number of changes
ntrees<-100 # desired number of trees
ntaxa<-100 # desired number of taxa
require(phytools); require(geiger) # required packages
mtrees<-list()
for(i in 1:ntrees){
# simulate a stochastic B-D tree using birthdeath.tree
tree<-drop.tip(birthdeath.**tree(b=1,d=0,taxa.stop=(ntaxa+**1)),as.character(ntaxa+1))
# compute the rate that results in (on average) nchanges changes
q<-nchanges/sum(tree$edge.**length)
# simulate a stochastic history
mtree<-sim.history(tree,Q=**matrix(c(-q,q,q,-q),2,2))
# loop until desired number of changes
while((length(unlist(mtree$**maps))-nrow(mtree$edge))!=**nchanges)
mtree<-sim.history(tree,Q=**matrix(c(-q,q,q,-q),2,2))
mtrees[[i]]<-mtree
# for fun, plot the history
plotSimmap(mtrees[[i]])
}
class(mtrees)<-"multiPhylo"
# tip states for the kth tree are in mtrees[[k]]$states
Please let us know if this is useful or if you have any questions or
difficulties implementing.
All the best, Liam
--
Antigoni Kaliontzopoulou
CIBIO, Centro de Investigacao em Biodiversidade e Recursos Geneticos
Campus Agrario de Vairao, 4485-661 Vairao
PORTUGAL
Department of Ecology, Evolution, and Organismal Biology
Iowa State University, Ames,
Iowa 50011, USA
tel: +351 91 3086188
mail to: [email protected]
[email protected]
______________________________**_________________
R-sig-phylo mailing list
[email protected]
https://stat.ethz.ch/mailman/**listinfo/r-sig-phylo<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
--
Rob Lanfear
Postdoc,
Centre for Macroevolution and Macroecology,
Research School of Biology,
Australian National University
Tel: +61 2 6125 7270
www.robertlanfear.com
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
--
Antigoni Kaliontzopoulou
CIBIO, Centro de Investigacao em Biodiversidade e Recursos Geneticos
Campus Agrario de Vairao, 4485-661 Vairao
PORTUGAL
Department of Ecology, Evolution, and Organismal Biology
Iowa State University, Ames,
Iowa 50011, USA
tel: +351 91 3086188
mail to: [email protected]
[email protected]
--
Antigoni Kaliontzopoulou
CIBIO, Centro de Investigacao em Biodiversidade e Recursos Geneticos
Campus Agrario de Vairao, 4485-661 Vairao
PORTUGAL
Department of Ecology, Evolution, and Organismal Biology
Iowa State University, Ames,
Iowa 50011, USA
tel: +351 91 3086188
mail to: [email protected]
[email protected]
_______________________________________________
R-sig-phylo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
