how to find the HPD 95% interval

[Huda Ghori]

Hi, 

I was interpreting my results and I dont know how to find the 95% HPD of a tMRCA. Specifically, i want to know how to find this:

" the majority of the Turkish sequences within clade A (pp = 0.51). The root of the tree dated back to the year 1910 (95% HPD 1875–1924). The two main clades A and B were respectively dated to 1924 (95% HPD 1895–1931) and 1917 (95% HPD 1893–1928). The origins of the Montenegrin groups included in subclade A and B were estimated to be 1932 (1906– 1942) and 1936 (1906–1938)."

Please help me

[Tim Vaughan]

Hi Huda,

If you're using BEAST 2, the Divergence Dating tutorial that can be found on the tutorials page of the website http://beast2.org/tutorials/ explains how to identify tMRCAs estimates and their 95% HPD intervals.

Cheers,

Tim

--
You received this message because you are subscribed to the Google Groups "beast-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to beast-users...@googlegroups.com.
To post to this group, send email to beast...@googlegroups.com.
Visit this group at http://groups.google.com/group/beast-users.
For more options, visit https://groups.google.com/d/optout.

[Huda Ghori]

Dear Tim ,

It is still not clear that how in many papers authors have mentioned the year range. From where will i get the year range, The tutorial asks the question regarding HPD and there is no answer. PLease help

--
You received this message because you are subscribed to a topic in the Google Groups "beast-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/beast-users/-wvgQdtek7E/unsubscribe.
To unsubscribe from this group and all its topics, send an email to beast-users...@googlegroups.com.

To post to this group, send email to beast...@googlegroups.com.
Visit this group at http://groups.google.com/group/beast-users.
For more options, visit https://groups.google.com/d/optout.

--

Noor ul Huda Ghori

M. Phil Scholar (Plant Biotechnology)

Atta-ur-Rahman School of Applied Biosciences

National University of Sciences & Technology (NUST)

Islamabad, Pakistan,

 

[Nick Matzke]

Hi Huda, this list is more for help with technical problems than for basic introductory training.  (you should really take a workshop or a course). But, since I have jet-lag, here's an basic outline of what you physically do:

(this is skipping all the Bayesian theory, which is very very important, and probably you shouldn't even run Beast until you've had enough training to understand what you are doing with Bayesian phylogenetics)

- Use BEAUTi to load your sequence data and set up your model. This produces an XML file.

- Use Beast2 to run the XML file. This produces a sampled trees file and a parameter log file.

- Use Tracer to look at the parameter logs, determine burnin and convergence. Tracer will also show you the 95% HPDs of the logged parameters.

- Use TreeAnnotator to take the post-burnin sampled trees and produce a Maximum Clade Credibility (MCC) tree.  This tree will have dates and also the 95% HPDs of the node dates

- Use FigTree to view the MCC tree, use the node bars option to view the 95% HPD of the node ages

Cheers,

Nick

[Huda Ghori]

Dear Nick,

I have followed exactly the same outline to generate trees. However, I am unable to find the HPD range in my tree. I have tried the node bar option but nothing displays and there is no change in the tree. I am sending you my MCC tree please see.

Thankyou

[Nick Matzke]

Ah - it wasn't clear in the initial email that the FigTree stage was where your problem was.

Anyhow -- FigTree1.4.0 gave me some problem loading your tree, but FigTree 1.4.2 seems to work.  Some of your taxon names start with "_", maybe that could cause an issue.

Files showing HPDs attached... Cheers, Nick

[Huda Ghori]

Dear Nick,

What i want is how to find this :

"(95% HPD 1895–1931)"

I cannot find this range in years. I do find the 95% HPD, values when i click on tree height but that is not in range of years. 

I hope you get what i am asking. Because i cant seem to find the answer anywhere

Regards,

[Nick Matzke]

Your questions are quite unclear!  I guess you are trying to get the dates *in calendar years*?  Are you trying to get them in a table, or on the plot?

The easiest method to just get the dates in a table is to click Scale Axis in FigTree to turn on the x-axis, click Reverse Axis, and then print it out, get a ruler, and write down the 95% HPDs.  Then transform years before present into calendar years using Excel.

It would be nice if the FigTree "Time Scale" tab would let you put in the appropriate "Offset by", you would need to put in -4000 or something to get the x-axis to be about right in years, but I can't get it to function on your tree beyond about -2000.

The automated method would be to use R to extract and transform all of the HPDs, and put them in a table. I could do it with BEASTmasteR but I am traveling.

You can also plot the phylogeny and HPDs in R with whatever timescale you want but this is also a custom R thing.  Email me off-list if you want to explore the R stuff.

Cheers, Nick

[Huda Ghori]

Dear Nick,

you mentioned this in your previous email:

"The easiest method to just get the dates in a table is to click Scale Axis in FigTree to turn on the x-axis, click Reverse Axis, and then print it out, get a ruler, and write down the 95% HPDs.  Then transform years before present into calendar years using Excel."

I have set my scale in calendar years. What i dont know is how to convert my 95%HPD in to calendar years..e.g. 1987 (95%HPD 1986-1990)

How do people find this value in calendar years. I dont seem to find this interval value in calendar years. 

Hope its clear now

[Nick Matzke]

Ages on the tree are in "age units before the present". The highest tip in the tree is considered to have age=0 units before the present.

So, ask yourself, what is the age of the highest tip in the tree?  If that sample was collected in year 2003, then year 2003 = 0 years before present.  If a node has an age of 10 units before the present, then that node has a calendar date of 2003-10 = 1993.

(I believe it's possible to set up the BEAST xml so that ages will be in years above the root, but I have not experimented with that, and it would involve re-doing the setup and then re-running everything else.)

Cheers,

Nick

[Huda Ghori]

so does this mean i have to calculate the HPD manually?

[chantal agbemabiese]

Hello Huda,

As Nick explained above, that's exactly how I go about it. Unfortunately I have had to calculate my HPDs manually all the time in a couple of publications, it's quite a task, I have been wondering if there's any easier way out.

Will try what Nick suggested and see how it goes!

All the best, cheers!

Chantal.

[Huda Ghori]

Dear Chantal,

Can you tell me exactly how you did that step wise. 

Thanks a ton

Regards,

[David Buckley]

If you define the nodes/clades in BEAUTi (in the “Taxa" option), the TMRCAs statistics for these nodes will be specified in the xml files and you will get the mean/median/stdev/HPDs etc for every node defined, which can be easily viewed and exported from TRACER.

dsvid

[Huda Ghori]

Dear David,

I am attaching my files can you please check and tell me what i need to do

[Santiago Sanchez]

What about the simple option in FigTree that shows 95% HPD on nodes? I'm not referring to the "display 95% HPD node height bars", but just 95% HPD node ages/heights on the node labels tab in FigTree. You would get a value like "(1920,1970)", which represents the node height range. You can get this range for any node your MCC tree without the need of any additional calculations/measurements.

Cheers,

Santiago

Santiago Sanchez-Ramirez

Ecology and Evolutionary Biology, University of Toronto

Natural History (Mycology), Royal Ontario Museum

100 Queen's Park

Toronto, ON

M5S 2C6

Canada

+1 416 586 8025

Website: www.sites.google.com/site/santiagosnchezrmirez/

<1b+NS5B.log.txt>

<1b+NS5B.tree>

[chantal agbemabiese]

Hello Santiago,

THanks for the suggestion but there is no "95% HPD node ages" in the lost of options in the dropdown menu of the Display of node labels tab. As you know, there is node age, height 95% HPD, etc. We can only display one at a time. 

Node age gives the age eg. 1999, the height 95% HPD gives the the interval Huda is looking for, however, these are not quoted as: (1980-2000) which Huda is interested in.

Huda, this is what I usually do, I just hope someone has a better way of doing this.

1. Display the height 95% HPD e.g. [53.45, 36.55], etc appears at each node.

2. Subtract these from the age of your youngest taxon e.g. if the latest taxon was detected in 2012, then you will have: (2012-53.45) to (2012-36.55). 

This is now your height 95% HPD in years as you wanted. THis was described by Nick.

3. I usually edit my tree which was saved in pdf format from figtree using adobe illustrator, so the tree with the node ages and node bars are edited by adding the intervals i calculated. 

You can convert your pdf to ppt using a free online help if you don't have illustrator, that way you can add the intervals to the ages in your figure.

I hope this helps. 

Cheers!

[Huda Ghori]

Dear Chantal,

Thankyou so much for the answer

Your answer is very very useful. You completely got my problem. Now just one thing i am confused about and want to double check. Do I have to subtract each HPD from the youngest node always?

Regards,

[chantal agbemabiese]

Dear Huda,

The 95% HPD interval is given in reference to the age of the youngest taxon, not the node. The nodes serve as the times of divergence from common ancestors. so the youngest node will definitely be older than the age of the youngest taxon. subtracting the HPD interval from the age of the youngest node will give you a wrong estimate of your HPD interval.

Answer: Deduct the HPD values from the age of the youngest taxon.

I hope this helps and some expert somewhere will comment on this if there is a better way of doing this. This will help all of us.

Cheers! 

[Huda Ghori]

Dear Chantal,

Thank you.  you solved my problem

Best regards

[Figueroa, Modesto]

Fffff
Sent from my Verizon Wireless 4G LTE Smartphone

Huda Ghori <huda...@gmail.com> wrote:

[Nick Matzke]

It looks like OutbreakTools has a solution to this also, but I have been getting questions about the below, so I have posted an example script for extracting Beast MCC tree statistics to a table, along with node ages, APE node numbers, posterior probabilities, clade memberships, etc.:

Example script: Extracting statistics from a Beast2 MCC tree

http://phylo.wikidot.com/beastmaster#extract_stats_script 

Cheers!

Nick

[Keith Barker]

FYI, Christoph Heibl's package "phyloch" also has much of this functionality (e.g., the function read.beast()).  It hasn't been updated since 2013, but it still seems to work, at least with BEAST 1.X trees (not sure if there are format details from some BEAST 2 output that might cause this to choke).

Keith

--
You received this message because you are subscribed to the Google Groups "beast-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to beast-users...@googlegroups.com.
To post to this group, send email to beast...@googlegroups.com.
Visit this group at http://groups.google.com/group/beast-users.
For more options, visit https://groups.google.com/d/optout.

-- 
F. Keith Barker, Ph.D.
Associate Professor, Department of Ecology, Evolution and Behavior
Curator of Genetic Resources, Bell Museum of Natural History
University of Minnesota
100 Ecology
1987 Upper Buford Circle
St. Paul, MN 55108
612.624.2737 (phone)
612.624.6777 (fax)
bark...@umn.edu
http://www.tc.umn.edu/~barke042

[Nick Matzke]

Yeah, part of the above script modifies stuff from phyloch, as stated in various places. I haven't checked the original phyloch in awhile, I don't think the MCC tree format has changed though. Part of the point in my script is to link up the MCC tree stats with the data.frame you can get from BioGeoBEARS::prt(tree), which puts everything in one nice big table corresponding to the R node numbering, default node order, ancestor/descendant nodes, etc. Mostly I did this because, most days, I am too dumb to parse the APE tree structure in my head.

Cheers,

Nick

[Chris Buddenhagen]

I was having a problem if I drop.tips and then ladderize, in terms of keeping support values and under BEAST parameters in the right order. Any ideas on how to implement that?

[Nick Matzke]

The function drop.tip() would lose all of the old node numbers etc. so it would take custom code to keep track of posterior probability values etc. for plotting. It would all depend on exactly what the goal is, so email me off list if you like. Cheers, Nick

Visit this group at https://groups.google.com/group/beast-users.