machu.tree.unc: Account for uncertainty in divergence time estimation

In wxguillo/machuruku: Ancestral niche reconstruction

Account for uncertainty in divergence time estimation

Description

Create a list of three trees that represent uncertainty in divergence times for downstream analysis

Usage

machu.tree.unc(tree, burnin = 0.1, conf = 0.95, inc = 10000, verbose = T)

Arguments

tree	Primary input. Can be a single treedata object (loaded via treeio::read.beast) or a string indicating a nexus-style treefile containing multiple trees (i.e. a Bayesian posterior distribution).
burnin	Percentage of trees to skip at the start of the multiple-trees file since MCMC algorithms take a while to converge. Must be between 0 and 1. Pertains only when the input is a Bayesian posterior of multiple trees. Default = 0.1.
conf	Confidence level from which to calculate HPD limits for the tree heights in a multiple-trees file. Must be between 0 and 1. Pertains only when the input is a Bayesian posterior of multiple trees. Default = 0.1.
inc	Increment at which to report progress, i.e. when every 'inc' trees are processed. Pertains only when the input is a Bayesian posterior of multiple trees, and verbose = TRUE. Default = 10000.
verbose	Report progress and checks to screen. Default = TRUE.

Details

This function generates a list of three trees that characterize the uncertainty of divergence time estimation. These trees can then be used in machu.2.ace() to explore scenarios in which varying numbers of taxa may have been present at a certain timeslice, given divergence time uncertainty.

The function has two modes, depending on what input is provided. When a single tree is provided, the function outputs the input tree, plus two additional trees, constructed from the upper and lower confidence limits provided to characterize divergence time uncertainty by time calibration software such as BEAST. The 'phylo' object (from ape::read.tree) usually used to represent phylogenetic trees in R does not store this information, so the input must be in the 'treedata' format (from treeio::read.beast). Currently the function is compatible with output from BEAST and BEAST-adjacent software such as SNAPP, as well as output from MEGA calibrated with the RelTime method.

The other mode characterizes uncertainty directly from a Bayesian posterior of trees. The function identifies the trees in the posterior with total tree heights closest to the median and the highest posterior density (HPD) limits of the distribution of tree heights. Because a tree posterior can consist of millions of trees, rather than overwhelm R by loading them all into an object (i.e. 'treedataList' from treeio), the function reads the trees directly from a file one at a time. Thus this mode is activated by specifying a filename; the file should be output from BEAST or BEAST-adjacent software, generally a .trees file in NEXUS format (open the file in a text editor to check, and examine the example file provided by the tutorial). The method consists of three "passes" through the posterior distribution of trees. In the first pass, the number of trees is simply counted, and the number of samples to skip as burn-in, given the burnin percentage specified by the user, is calculated. In the second pass, total tree heights are calculated for each tree, and the median and HPD limits at the confidence level specified are calculated for this distribution. In the third pass, the trees with total heights closest to the median and HPD limits are found, and sent to the output list.

Because different software or even versions of the same software may output subtly variable versions of the same formats, I can't guarantee that this function can read in your multiple-tree nexus file. In some cases the issue may be easily fixed with a patch - contact me (WXG) with bugs.

Value

A list of three 'phylo' objects characterizing uncertainty in divergence times.

Examples

# load a tree and store divergence time uncertainty info
tree <- treeio::read.beast("tree.treefile")
# characterize uncertainty
trees <- machu.tree.unc(tree)
# plot trees to visualize uncertainty
machu.treeplot(trees)

# characterize uncertainty from a posterior of trees and treat the first 20% as burn-in
trees <- machu.tree.unc("posterior.trees", burnin=0.2)

wxguillo/machuruku documentation built on Jan. 23, 2025, 3:25 p.m.