MSCquartets (version 3.1)

estimateEdgeLengths: Estimate edge lengths on a species tree from gene tree quartet counts

Description

Estimate edge lengths, in coalescent units, on an unrooted species tree from a table of resolved quartet counts from a collection of gene trees.

Usage

estimateEdgeLengths(tree, rqt, terminal = 1, method = "simpleML", lambda = 1/2)

Value

an unrooted metric tree with the same topology as tree, of type phylo

Arguments

tree

a phylo object, giving a resolved tree on which to estimate edge lengths

rqt

a resolved quartet table, as from quartetTableResolved, in which all taxa on tree appear

terminal

an edge length to assign to terminal edges, whose lengths cannot be estimated

method

"simpleML" or "simpleBayes"

lambda

a positive parameter for the "simpleBayes" method

Details

While the argument tree may be supplied as rooted or unrooted, metric or topological, only its unrooted topology will be used for determining new metric estimates.

Counts of quartets for all those quartets which define a single edge on the tree (i.e., whose internal edge is the single edge on the unrooted input tree) are summed, and from this an estimate of the branch length is computed. If method= "simpleML" this is the maximum likelihood estimate. If method="simpleBayes" this is the Bayesian estimate of Theorem 2 of SayMir16;textualMSCquartets, using parameter lambda. Using lambda=1/2 gives a flat prior on [1/3,1] for the probability of the quartet displayed on the species tree.

These methods are referred to as `simple' since they use only the quartets defining a single edge of the species tree. Quartets with central edges composed of several edges in the species tree are ignored.

Note that branch length estimates may be 0 (if the count for the quartet displayed on the input tree is not dominant), positive, or Inf (if the counts for quartet topologies not displayed on the input tree are all 0, and method="simpleML").

References

SayMir16MSCquartets

Examples

gtrees=read.tree(file=system.file("extdata","dataGeneTreeSample",package="MSCquartets"))
taxanames=taxonNames(gtrees)
QT=quartetTable(gtrees,taxanames[1:6])
RQT=quartetTableResolved(QT)
DQT=quartetTableDominant(RQT)
tree=QDS(DQT)
write.tree(tree)
plot(tree)
metricMTree=estimateEdgeLengths(tree,RQT,method="simpleML")
write.tree(metricMTree)
plot(metricMTree)
metricBTree=estimateEdgeLengths(tree,RQT,method="simpleBayes")
write.tree(metricBTree)
plot(metricBTree)