Hi,
Klaus beat me to it, but since there are more ways to do it ;-)
Two functions in the APE package are handy: is.monophyletic(), and
all.equal.phylo(). See the help for options.
They could be used in wrapper functions, like this
filter.mono - function (x, grp) {
fun - function (x) all(is.monophyletic(x, grp))
mono - x[sapply(x, fun)]
invisible(mono)
}
filter.backbone - function (x, y) {
fun - function (x) all(all.equal.phylo(x, y, use.edge.length = FALSE))
equal - x[sapply(x, fun)]
invisible(equal)
}
# And to test
# Filter trees based on presence of monophyletic group
trees - rmtree(200, 5)
group - c(t1, t2)
monotrees - filter.mono(trees, group)
# Filter trees based on compatible with backbone tree
backbone - trees[[1]]
btrees - filter.backbone(trees, backbone)
Cheers
Johan
P.S. Don't know how efficient these approaches are on large data. Could be
worth comparing. D.S.
On Fri, 2011-03-25 at 09:47 +0100, Klaus Schliep wrote:
Dear Brad,
here are some solutions for your problems. I assume unrooted trees
here and ignore the position of the root, so you may additionaly check
for the position of the root.
On 3/24/11, Ruhfel Brad ruh...@gmail.com wrote:
Dear all,
I am wondering if there is some way using R to sort a file of multiple trees
(i.e., BEAST output) to either
1) return all trees that match a constraint topology (backbone or full
topology) or
Assume you have all the trees in given in a list of trees (phylo
object) or as multiPhylo object (ape package) and tree0 is the tree
you want them to be compared with:
library(ape)
library(phangorn)
# generate some toy trees
trees = rmtree(10, 5, rooted = FALSE)
tree0 = trees[[5]]
#
tmp = sapply(trees, RF.dist, tree0) # compute Robinson Foulds distance
result1 = trees[tmp==0] # the trees you want
2) return all trees that have a monophyletic clade of a list of taxa.
# Define your partition, here (t1,t2) go together
dat = matrix(c(1,1,0,0,0),dimnames=list(c(t1, t2, t3, t4,
t5),NULL) , ncol=1)
dat = phyDat(dat, type=USER, levels = c(0,1), ambiguity=NULL)
tmp2 = parsimony(trees, dat)
result2 = trees[tmp2==1] #
Any help much appreciated!
Just joined the list... looking forward to it!
-Brad
PhD Student
Harvard Herbaria
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Regards,
Klaus
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
