Re: [R-sig-phylo] sort trees by constraint?
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 -- Klaus Schliep Université Paris 6 (Pierre et Marie Curie) 9, Quai Saint-Bernard, 75005 Paris ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] sort trees by constraint?
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
Re: [R-sig-phylo] reading nexus file from treebase?
Thanks for tracking this down, Emmanuel. I'm forwarding this to the TreeBASE developers list for consideration. BTW the phylobase implementation uses NCL (NEXUS Class Library, in C+ +) for parsing NEXUS. Given the difference, probably ape does not? -hilmar On Mar 25, 2011, at 12:00 AM, Emmanuel Paradis wrote: Scott, readNexus (phylobase) can read your tree but not read.nexus (ape). The problem is the two lines inserted within the TREES block: BEGIN TREES; TITLE Tb10793; LINK TAXA = M4787; TREE Fig._3c = [R] If you delete them, it's OK. FigTree also cannot read this file. Cheers, Emmanuel François Michonneau wrote on 25/03/2011 00:35: Hi Scott, Which version of phylobase are you using and which architecture? I can read the file on my machine. Cheers, -- François On Thu, Mar 24, 2011 at 12:39, Scott Chamberlain myrmecocys...@gmail.com wrote: Hello, I can't get read.nexus (ape) or readNexus (phylobase) to read nexus files downloaded from treebase with URLs parsed from xml files. I can't manually edit each file as I want to read a lot of these files. Is there an easy fix? One of the files is copied below. Thanks! Scott Chamberlain Rice University, EEB Dept. #NEXUS [!This data set was downloaded from TreeBASE, a relational database of phylogenetic knowledge. TreeBASE has been supported by the NSF, Harvard University, Yale University, SDSC and UC Davis. Please do not remove this acknowledgment from the Nexus file. Downloaded on March 24, 2011; 16:32 GMT TreeBASE (cc) 1994-2008 Study reference: Brown R., Yang Z. 2010. Bayesian Dating of Shallow Phylogenies with a Relaxed Clock. Systematic Biology, 59(2): 119-131. TreeBASE Study URI: http://purl.org/phylo/treebase/phylows/study/TB2:S10165] BEGIN TAXA; TITLE M4787; DIMENSIONS NTAX=16; TAXLABELS Chalcides_coeruleopunctatus_E2806.20 Chalcides_coeruleopunctatus_E2806.22 Chalcides_manueli_E2506.1 Chalcides_mionecton_mionecton_E2506.10 Chalcides_mionecton_mionecton_E2506.12 Chalcides_mionecton_trifasciatus_E2506.18 Chalcides_polylepis_E14124.1 Chalcides_polylepis_E14124.2 Chalcides_polylepis_E2506.21 Chalcides_sexlineatus_bistriatus_E2806.6 Chalcides_sexlineatus_sexlineatus_E2806.8 Chalcides_simonyi_E3007.2 Chalcides_sphenopsiformis_E8121.26 Chalcides_sphenopsiformis_E8121.27 Chalcides_viridanus_E2806.10 Chalcides_viridanus_E2806.14 ; END; BEGIN TREES; TITLE Tb10793; LINK TAXA = M4787; TREE Fig._3c = [R] ((Chalcides_sphenopsiformis_E8121.26,Chalcides_sphenopsiformis_E8121.27),(((Chalcides_viridanus_E2806.10,Chalcides_viridanus_E2806.14),((Chalcides_sexlineatus_bistriatus_E2806.6,Chalcides_sexlineatus_sexlineatus_E2806.8),(Chalcides_coeruleopunctatus_E2806.22,Chalcides_coeruleopunctatus_E2806.20))),(Chalcides_simonyi_E3007.2,((Chalcides_mionecton_trifasciatus_E2506.18,(Chalcides_mionecton_mionecton_E2506.12,Chalcides_mionecton_mionecton_E2506.10)),(Chalcides_manueli_E2506.1,(Chalcides_polylepis_E14124.1,(Chalcides_polylepis_E14124.2,Chalcides_polylepis_E2506.21))); [! TreeBASE tree URI: http://purl.org/phylo/treebase/phylows/tree/TB2:Tr6136] END; [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- Emmanuel Paradis IRD, Jakarta, Indonesia http://ape.mpl.ird.fr/ ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- === : Hilmar Lapp -:- Durham, NC -:- informatics.nescent.org : === ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] summarize lagrange results on a target tree?
Dear all, I am trying to figure out how to automate something in R. Any suggestions much appreciated. In general I am trying to summarize the results of ancestral area reconstructions done using Lagrange on 1000 trees from BEAST. These results would then be summarized and reported based on the nodes present in an optimal target tree. The trouble is, the nodes present in the target tree are not necessarily present in all trees in the 1000 BEAST trees due to phylogenetic uncertainly at some nodes. This is a similar approach to that used in Nylander et al (2008) but instead of using DIVA to do the ancestral area reconstructions I am using Lagrange. Here is what I am trying to figure out how to do in R. 1) Read in a target tree with nodes that are numbered (Results will be visualized on this tree) 2) List the MCRAs of each node. 3) Read in the 1000 Beast trees that vary in topology and branch lengths. 3) Starting with node 1 (or 0) in the target tree Search the distribution of 1000 BEAST trees for trees that contain the node of interest based on the MRCA list in the target tree. For nodes that are not maximally supported this will be a fraction of the 1000 trees. 4) Summarize the frequency of each reconstruction at that node in that subset of trees and save to a table ( for example: Node 1 30% area A, 30% Area B, 30% area AB) 5) repeat for each node in the target tree. 6) draw the target tree with pie charts at each node with the frequency of each reconstruction Looking forward to any thought on how to do any part of this! Best, -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
[R-sig-phylo] Possible Error in prop.part() ?
Hello all, I was using prop.part() to identify nodes shared between several topologies, when I discovered a strange error, where prop.part() was underestimating the number of partitions found in common. I was able to fix the error for some comparisons by outputting the trees as Newick strings and reading them back in with read.tree(text=write.tree()). However, this didn't work for all topologies, which actually allowed me to reproduce the error for you all, using the code below. tree1 and tree2 have five shared nodes, but prop.part() only reports two partitions as being found in common. prop.clades() seems to report the correct number. require(ape) tree1-read.tree(text=(((TAXt9,TAXt4),(((TAXt6,TAXt3),TAXt13),((TAXt1,(TAXt12,TAXt14)),TAXt11))),TAXt7);) tree2-read.tree(text=(((TAXt6,((TAXt13,TAXt3),(TAXt1,(TAXt12,TAXt14,(TAXt11,(TAXt9,TAXt4))),TAXt7);) prop.part(tree1,tree2) prop.clades(tree1,tree2) For this particular example, it turns out that prop.part() correctly estimated the number of partitions prior to me applying the read.tree() fix. But I'm doing this over a large number of pairs, so I need a solution that works for all. I'm using R v2.12.2 and ape v2.6-3. I haven't updated to 2.7 since I use read.nexus() quite a bit; let me know if it is reproducible in that version. -Dave -- David Bapst Dept of Geophysical Sciences University of Chicago 5734 S. Ellis Chicago, IL 60637 http://home.uchicago.edu/~dwbapst/ [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] Pagel's lambda greater than one?
Hmmm. I am not familiar with what this program is doing, but I thought that Pagel's lambda is supposed to be constrained in the range of zero to one. Note that Grafen's (1989) rho and the OU transform (d) or ACDC transform (g) of Blomberg et al. (2003) and Lavin et al. (2008) can exceed one. When those transform parameters exceed one, it puches the interior nodes up closer to the terminal taxa (tips) and hence makes the tree more hierarchical than the starter tree. (But things can get a little strange if your starter tree does not have contemporaneous tips [not ultrametric].) Cheers, Ted Original message Date: Fri, 25 Mar 2011 12:11:16 +0100 From: Christopher Turbill christopher.turb...@fiwi1.fiwi.at Subject: [R-sig-phylo] Pagel's lambda greater than one? To: r-sig-phylo@r-project.org I have a question regarding using the corPagel correlation structure in a gls model with package 'ape'. The model is simple, with body masss as a covariate and a single trait the predictor. Pagel's lambda is estimated to be 1.18422, which appears to turn a 'gunshot' relationship into a highly significant one. Forcing lamda to be = 1, makes the predictor very non-significant (as it is in the non-phylogenetically informed analysis). My question is: what does a lambda 1 really mean? Should I restrict the analysis to lambda = 1 when the 'best' lambda is estimated to be 1? Thanks for you help, Chris Turbill ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo