[R-sig-phylo] Plotting ancestral reconstructions on subtrees
Dear all, I am attempting to plot ancestral reconstructions (using asr.marginal in Diversitree) of a discrete trait on a number of subtrees (subtree=T in ape), but I'm having trouble assigning node probabilites and tip states to the correct nodes and tips. I have been using the following code: #Extract state probabilitiesst - asr.marginal(lik.l, coef(fit.l))pie.probs-t(st) pr-apply(t(st), 1, which.max)tree$node.label-pr #Assign colours to extant statesstatesLabel - statesnames(statesLabel) - names(states)statesLabel[states==1] - redstatesLabel[states==2] - bluestatesLabel[states==3] - green #Extract subtreetrd - drop.tip(tree, 65:761, subtree=T) #Plot subtreeplot(trd, show.tip.label=T, edge.width =1)nodelabels(pie=pie.probs,piecol=c(red, blue, green))tiplabels(pch=22, bg=statesLabel) This works fine for my first subtree, i.e. tips 1-64, but when I try running the same code on subsequent parts of the tree, visual inspection confirms that tip colours and node probabilities are no longer correctly assigned to the correct tips and nodes. #I.e. running the following code plots the correct subtree, but states are incorrectly assigned to nodes and tips. trd - drop.tip(tree, c(1:64, 167:761), subtree=T)plot(trd, show.tip.label=T, edge.width =1)nodelabels(pie=pie.probs,piecol=c(red, blue, green))tiplabels(pch=22, bg=statesLabel) I'm guessing the solution is quite simple, but I haven't been able to work it out. Any suggestions would be greatly appreciated. Kind regards,Petter Zahl Marki [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] atomic vector for write.nexus
Hello, This is Elaine. I am using the command write.nexusin package (ape) to generate a nexus file after pruning a phylogeny tree. However, an error appeared, saying Error in obj[[i]]$tip.label : $ operator is invalid for atomic vectors I searched the archive but found no similar experience in this command. Please kindly help explain what it means and how to solve it. code library(ape) # read data and tree birddata - read.csv(H:/birddata_family_20130405.csv, header = TRUE) birdtree - read.nexus(H:/birddata_family_20130405.nexus) setwd(H:) species.to.keep -read.csv(H:/birddata_family_20130405.csv, header = TRUE) pruned.tree-drop.tip(birdtree,birdtree$tip.label[- na.omit(match(species.to.keep[,1],birdtree$tip.label))]) ## prune the data to include only those taxa in the pruned tree data.pruned - birddata[tree.pruned$tip.label, ] (Thanks to Klaus and others) write.nexus(tree.pruned, H:/pruneddata_20130405.nexus, translate = FALSE) Elaine [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] Fwd: atomic vector for write.nexus
I changed pruned.tree to tree.pruned. The error still existed. Elaine -- Forwarded message -- From: Elaine Kuo elaine.kuo...@gmail.com Date: Fri, Apr 5, 2013 at 6:12 PM Subject: atomic vector for write.nexus To: r-sig-phylo@r-project.org Hello, This is Elaine. I am using the command write.nexusin package (ape) to generate a nexus file after pruning a phylogeny tree. However, an error appeared, saying Error in obj[[i]]$tip.label : $ operator is invalid for atomic vectors I searched the archive but found no similar experience in this command. Please kindly help explain what it means and how to solve it. code library(ape) # read data and tree birddata - read.csv(H:/birddata_family_20130405.csv, header = TRUE) birdtree - read.nexus(H:/birddata_family_20130405.nexus) setwd(H:) species.to.keep -read.csv(H:/birddata_family_20130405.csv, header = TRUE) tree.pruned-drop.tip(birdtree,birdtree$tip.label[- na.omit(match(species.to.keep[,1],birdtree$tip.label))]) ## prune the data to include only those taxa in the pruned tree data.pruned - birddata[tree.pruned$tip.label, ] (Thanks to Klaus and others) write.nexus(tree.pruned, H:/pruneddata_20130405.nexus, translate = FALSE) Elaine [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] pruning a tree
Hello This Elaine. I tried to prune a phylogeny tree using two methods based on the code attached. Method 1 returned the tip.label sharing between birddata and birdtree. Method 2 returned nothing. Please kindly indicate why Method 2 failed to prune the tree. Also, please kindly explain the command birdtree$tip.label[-na.omit Thank you again. Code Method 1 Library(ape) birddata - read.csv(H:/ birddata_family_20130405.csv, header = TRUE) birdtree - read.nexus(H:/ birddata_family_20130405.nexus) setwd(H:) ## prune the phylogeny to include only species which are in the data set species.to.keep-read.csv(H:/ birddata_family_20130405.csv, header = TRUE) birdtree.pruned-drop.tip(birdtree,birdtree$tip.label[-na.omit(match(species.to.keep[,1],birdtree$tip.label))]) Method 2 Library(ape) birddata - read.csv(H:/ birddata_family_20130405.csv, header = TRUE) birdtree - read.nexus(H:/ birddata_family_20130405.nexus) setwd(H:) ## prune the phylogeny to include only species which are in the data set temp - name.check(birdtree, birddata) birdtree.pruned - drop.tip(birdtree, tip=temp$Tree.not.data) [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] ML ancestral state reconstruction using different softwares
Jingchun, Thanks to Daniel for giving a detailed answer to your first question. You can use ace() to fit a model with irreversible change. You need to build a custom model like this: Q - matrix(0, 2, 2) Q[2, 1] - 1 Q [,1] [,2] [1,]00 [2,]10 In this case there's one rate to be estimated and the other is fixed to zero (the 0's on the diagonal are ignored). This can be generalized to more complicated models (there's an example in my book). Then call ace with: ace(traits_data, tree, type = d, model = Q) This time you cannot compare this model with the ER one with an LRT: compare directly the log-liks since they have the same number of parameters or use AIC(). BTW, there is an anova methods to compare nested models fitted with ace: this avoids the pain to compute correctly df's. Cheers, Emmanuel -Original Message- From: Jingchun Li jingc...@umich.edu Date: Thu, 4 Apr 2013 22:25:09 To: emmanuel.para...@ird.fr; r-sig-phylo@r-project.org Subject: Re: [R-sig-phylo] ML ancestral state reconstruction using different softwares Sorry, I forgot to include the list. there we go. -- *From: * Jingchun Li jingc...@umich.edu *Date: *Thu, 4 Apr 2013 21:58:56 -0400 *To: *emmanuel.para...@ird.fr *Subject: *Re: [R-sig-phylo] ML ancestral state reconstruction using different softwares Thank you very much, Dr.Paradis. You are right that the two parameter model does not fit the data significantly better. I guess I'm more puzzled about why the three softwares gave me different likelihoods. I did fit the symmetrical model using BayesTraits later and it gave me about the same transition rates as the asymmetrical model and a similar likelihood value (-14.08 for symmetrical and -14.00 for asymmetrical). However I do have biological evidence that transition from 0 to 1 might be extremely hard. So maybe I can also try set q01 to 0 (which is actually what ace gave me under the asymmetrical model)? Cheers, Jingchun On Thu, Apr 4, 2013 at 9:21 PM, Emmanuel Paradis emmanuel.para...@ird.frwrote: Hi, You cannot really say that in both ace and Mesquite, it slightly favors the asymmetrical model: the increase in log-likelihood is less than 1 which is far from being significant with one additional parameter. So it seems that all three pieces of software agree well on the estimate of rate for the symmetrical model. You should fit also this model with BayeTraits to do the full comparison. Cheers, Emmanuel -Original Message- From: Jingchun Li jingc...@umich.edu Sender: r-sig-phylo-boun...@r-project.org Date: Thu, 4 Apr 2013 10:15:05 To: r-sig-phylo@r-project.org Subject: [R-sig-phylo] ML ancestral state reconstruction using different softwares Hi all, I am exploring different methods for ancestral state reconstruction for a small phylogeny I have (26 OTUs). I have one binary trait coded as 0s and 1s for all the taxa. And I am more interested in the transition rates between the two states. If I understand correctly, a maximum likelihood ancestral state reconstruction using a MK1 model or an asymmetrical rates model should give me the optimized estimated transition rates between the two states. The formal will assume the two rates are the same, the later will assume independent rates. So I firstly did this using the ape function in ace. ERreconstruction - ace(traits_data, tree, type=discrete, model=ER) ERreconstruction$rates ERreconstruction$loglik ARDreconstruction - ace(traits_data, tree, type=discrete, model=ARD) ARDreconstruction$rates ARDreconstruction$loglik Then I tried using Mesquite, using MK1 model and the Asymmetrical 2-param. Markov-k Model Then I tried using BayesTraits, with the MultiState - Maximum Likelihood option, with no restrictions on q10 and q01. What puzzled me is that I'm getting different results from the three softwares. -- aceER: transition rate: 0.0059, likelihood: -7.00 aceARD: transition rate: 0.019 and 0.00, likelihood: -6.11 Mesquite MK1: transition rate: 0.0059, likelihood: -7.70 Mesquite Asymmetrical: transition rate: 0.015 and 0.005, likelihood: -7.30 BayesTraits no restriction: transition rate: 0.0057 and 0.0057, likelihood: -14.09 I can see that in both ace and Mesquite, it slightly favors the asymmetrical model with two different rates, and the rates and likelihoods are more of less comparable. But BayesTraits seems to think the two rates should be equal, and it has a different likelihood. Does this has something to do with different searching algorithms? Or am I missing something here? I am aware that the scaled likelihoods in ace are the scaled conditional likelihoods, not the joint or marginal reconstructions. But this should not affect the estimated transition rates? Thank you very much! [[alternative HTML version deleted]]
Re: [R-sig-phylo] ML ancestral state reconstruction using different softwares
Thank you very much Daniel and Emmanuel! Now things start to make sense. I guess my last question is, if there is biological evidence that the transition rates between the two traits are different, is it even meaningful to compare the two models? Should I simply avoid using the equal rates model? Cheers, Jingchun On Fri, Apr 5, 2013 at 9:12 AM, Emmanuel Paradis emmanuel.para...@ird.frwrote: ** Jingchun, Thanks to Daniel for giving a detailed answer to your first question. You can use ace() to fit a model with irreversible change. You need to build a custom model like this: Q - matrix(0, 2, 2) Q[2, 1] - 1 Q [,1] [,2] [1,] 0 0 [2,] 1 0 In this case there's one rate to be estimated and the other is fixed to zero (the 0's on the diagonal are ignored). This can be generalized to more complicated models (there's an example in my book). Then call ace with: ace(traits_data, tree, type = d, model = Q) This time you cannot compare this model with the ER one with an LRT: compare directly the log-liks since they have the same number of parameters or use AIC(). BTW, there is an anova methods to compare nested models fitted with ace: this avoids the pain to compute correctly df's. Cheers, Emmanuel -- *From: * Jingchun Li jingc...@umich.edu *Date: *Thu, 4 Apr 2013 22:25:09 -0400 *To: *emmanuel.para...@ird.fr; r-sig-phylo@r-project.org *Subject: *Re: [R-sig-phylo] ML ancestral state reconstruction using different softwares Sorry, I forgot to include the list. there we go. -- *From: * Jingchun Li jingc...@umich.edu *Date: *Thu, 4 Apr 2013 21:58:56 -0400 *To: *emmanuel.para...@ird.fr *Subject: *Re: [R-sig-phylo] ML ancestral state reconstruction using different softwares Thank you very much, Dr.Paradis. You are right that the two parameter model does not fit the data significantly better. I guess I'm more puzzled about why the three softwares gave me different likelihoods. I did fit the symmetrical model using BayesTraits later and it gave me about the same transition rates as the asymmetrical model and a similar likelihood value (-14.08 for symmetrical and -14.00 for asymmetrical). However I do have biological evidence that transition from 0 to 1 might be extremely hard. So maybe I can also try set q01 to 0 (which is actually what ace gave me under the asymmetrical model)? Cheers, Jingchun On Thu, Apr 4, 2013 at 9:21 PM, Emmanuel Paradis emmanuel.para...@ird.fr wrote: Hi, You cannot really say that in both ace and Mesquite, it slightly favors the asymmetrical model: the increase in log-likelihood is less than 1 which is far from being significant with one additional parameter. So it seems that all three pieces of software agree well on the estimate of rate for the symmetrical model. You should fit also this model with BayeTraits to do the full comparison. Cheers, Emmanuel -Original Message- From: Jingchun Li jingc...@umich.edu Sender: r-sig-phylo-boun...@r-project.org Date: Thu, 4 Apr 2013 10:15:05 To: r-sig-phylo@r-project.org Subject: [R-sig-phylo] ML ancestral state reconstruction using different softwares Hi all, I am exploring different methods for ancestral state reconstruction for a small phylogeny I have (26 OTUs). I have one binary trait coded as 0s and 1s for all the taxa. And I am more interested in the transition rates between the two states. If I understand correctly, a maximum likelihood ancestral state reconstruction using a MK1 model or an asymmetrical rates model should give me the optimized estimated transition rates between the two states. The formal will assume the two rates are the same, the later will assume independent rates. So I firstly did this using the ape function in ace. ERreconstruction - ace(traits_data, tree, type=discrete, model=ER) ERreconstruction$rates ERreconstruction$loglik ARDreconstruction - ace(traits_data, tree, type=discrete, model=ARD) ARDreconstruction$rates ARDreconstruction$loglik Then I tried using Mesquite, using MK1 model and the Asymmetrical 2-param. Markov-k Model Then I tried using BayesTraits, with the MultiState - Maximum Likelihood option, with no restrictions on q10 and q01. What puzzled me is that I'm getting different results from the three softwares. -- aceER: transition rate: 0.0059, likelihood: -7.00 aceARD: transition rate: 0.019 and 0.00, likelihood: -6.11 Mesquite MK1: transition rate: 0.0059, likelihood: -7.70 Mesquite Asymmetrical: transition rate: 0.015 and 0.005, likelihood: -7.30 BayesTraits no restriction: transition rate: 0.0057 and 0.0057, likelihood: -14.09 I can see that in both ace and Mesquite, it slightly favors the asymmetrical model with two different rates, and the rates and likelihoods are more of less comparable. But BayesTraits seems to think the two rates should be
Re: [R-sig-phylo] pruning a tree
Hi Elaine. A simpler way to do method 1 that doesn't require na.omit or match is: # if there are species in the tree that are missing from the data: ss-species.to.keep[,1] # species in the data tree-birdtree tree.pruned-drop.tip(tree,setdiff(tree$tip.label,ss)) I would guess that method 2 isn't working because your species names are not row names in your data matrix. To have that, you could do: birddata-read.csv(H:/birddata_family_20130405.csv,header=T,row.names=1) All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 4/5/2013 8:54 AM, Elaine Kuo wrote: Hello This Elaine. I tried to prune a phylogeny tree using two methods based on the code attached. Method 1 returned the tip.label sharing between birddata and birdtree. Method 2 returned nothing. Please kindly indicate why Method 2 failed to prune the tree. Also, please kindly explain the command “birdtree$tip.label[-na.omit” Thank you again. Code Method 1 Library(ape) birddata - read.csv(H:/ birddata_family_20130405.csv, header = TRUE) birdtree - read.nexus(H:/ birddata_family_20130405.nexus) setwd(H:) ## prune the phylogeny to include only species which are in the data set species.to.keep-read.csv(H:/ birddata_family_20130405.csv, header = TRUE) birdtree.pruned-drop.tip(birdtree,birdtree$tip.label[-na.omit(match(species.to.keep[,1],birdtree$tip.label))]) Method 2 Library(ape) birddata - read.csv(H:/ birddata_family_20130405.csv, header = TRUE) birdtree - read.nexus(H:/ birddata_family_20130405.nexus) setwd(H:) ## prune the phylogeny to include only species which are in the data set temp - name.check(birdtree, birddata) birdtree.pruned - drop.tip(birdtree, tip=temp$Tree.not.data) [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] Seeking list of nucleotide substitution models
Hello, Is there a review or list of ~every specific nucleotide substitution model that has been proposed or used in the literature (with references)? I'm interested in reversible models. The most exhaustive list I have is from the jModeltest documentation, https://code.google.com/p/jmodeltest2/wiki/TheoreticalBackground#Models_of_ nucleotide_substitution I realise an enormous number of models is possible (and may have been used in model averaging). I'm keen to know which have some precedent in the literature, along the lines of those listed above. I'm finding it difficult to get beyond the standard sources. Thank you in advance, Daniel -- Daniel Barker http://bio.st-andrews.ac.uk/staff/db60.htm The University of St Andrews is a charity registered in Scotland : No SC013532 ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] read.nexus.data parser
Hi All, I am wondering if there is anyway to increase the speed of the read.nexus.data parser. Or if there is an alternative that is a faster nexus file data parser. THanks, Jess -- Jessica L. Sabo js...@ufl.edu sabo.j...@gmail.com PhD Student, The Edward Braun Lab Department of Biology, Zoology Program University of Florida P.O. Box 118525 Gainesville, FL 32611 [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] Exporting Ancestor Reconstruction from Mesquite to R for ouch
Dear forum, I am trying to conduct an ouch analysis and for it I must determine the ancestral selective regime at each node in my tree. To do so I need to conduct ancestral reconstruction on discrete unordered data. To my knowledge no package can do this in R (I hope I am wrong) and thus I used Mesquite's function to 'Trace Character History'. The problem is exporting this information to R. In Mesquite the only export option I see is 'Export Ancestral State Trace' but the only option is SIMMAP 1.5 format. I tried opening this file using phytools'read.simmap but it just hangs. I can open the file with FigTree and I can see the Nodel labels but when I export the tree as a nexus file it does not save the model labels. Any advice would be appreciated. thanks, Mart Martin Turcotte, Ph. D. Dept. of Biology University of Toronto at Mississauga 3359 Mississauga Road North, Mississauga, Ontario, Canada, L5L 1C6 http://individual.utoronto.ca/martinturcotte [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] Exporting Ancestor Reconstruction from Mesquite to R for ouch
Hi Martin. Without more information or your input file, it's hard for me to evaluate why read.simmap is hanging. Did you run: tree-read.simmap(filename,format=nexus,version=1.5) If you have your discrete character data in R, you can also use the phytools function make.simmap to generate stochastic maps. Make sure that you are using a recent version of phytools (=0.2-26; =0.2-33 for uncertain/unknown tip states; =0.2-37 for non-symmetric transition matrices) as earlier versions are buggy. My advice would be to use make.simmap to generate stochastic maps (or read them from file), and then use OUwie (in the package OUwie) to fit a multiple optima OU model to each mapped tree. For discrete character x and continuous trait y (each vectors in the same order with names(x)=names(y)=species names); and a phylogeny, tree, do: library(phytools) library(OUwie) mtrees-make.simmap(tree,x,nsim=100) XX-data.frame(names(x),x,y) result-lapply(mtrees,OUwie,data=XX,model=OUM,simmap.tree=TRUE) This just gives you a giant list of the results returned by OUwie for each tree. It would probably be a good idea to organize your results of interest more sensibly than this so that they can be easily aggregated across mappings. Good luck. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 4/5/2013 6:13 PM, Martin Turcotte wrote: Dear forum, I am trying to conduct an ouch analysis and for it I must determine the ancestral selective regime at each node in my tree. To do so I need to conduct ancestral reconstruction on discrete unordered data. To my knowledge no package can do this in R (I hope I am wrong) and thus I used Mesquite's function to 'Trace Character History'. The problem is exporting this information to R. In Mesquite the only export option I see is 'Export Ancestral State Trace' but the only option is SIMMAP 1.5 format. I tried opening this file using phytools'read.simmap but it just hangs. I can open the file with FigTree and I can see the Nodel labels but when I export the tree as a nexus file it does not save the model labels. Any advice would be appreciated. thanks, Mart Martin Turcotte, Ph. D. Dept. of Biology University of Toronto at Mississauga 3359 Mississauga Road North, Mississauga, Ontario, Canada, L5L 1C6 http://individual.utoronto.ca/martinturcotte [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] R-Heatmap--- node stack overflow
Dear All R-sig-phylo users, I was using R ploting a heatmap recently. The datatable consisting 2400 rows and 6 columns, which I do not think is too big for R to handle, however I keep receiving the error message node stack overflow. Does anyone have similar problem? Any help will be appreciated. Thanks. With regards, Yan Wang, Ph.D student Dr. Moncalvo's Group Dept. Ecology Evolutionary Biology, University of Toronto Dept. Natural History, Royal Ontario Museum, 100 Queens Park Toronto, ON M5S 2C6 Office: (416)586-8025 Skype: yan.wang.fungi Nothing in biology makes sense except in the light of evolution Theodosius Dobzhansky, 1977 [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
