Matrix representation parsimony (MRP) supertree estimation

Are any of my readers aware of an existing function to do Matrix Representation Parsimony (MRP) supertree estimation (Baum 1992; Ragan 1992; also see Sanderson et al. 1998; Bininda-Emonds et al. 2007) in R? I'm not, but I was thinking it would be a fun - and easy - thing to program. Consequently I have just posted a function for this, called mrp.supertree(), on my R-phylogenetics page (direct link to code here). This function could also be used to get the MRP consensus tree (that is, if all the input phylogenies have the same set of taxa; e.g., see Felsenstein 2004 pp. 527-528).

I programmed this by using two handy functions: the function prop.part() in {ape} (which returns all the non-trivial bipartitions for a given tree) and pratchet() in {phangorn} (Schliep 2010), which estimates the MP tree using the parsimony ratchet method of Nixon (1999).

The way MRP supertree estimation works is by first constructing N (for N trees) n_i x m_i binary matrices (for n_i species and m_i non-trivial bipartitions in the ith tree). We score the tree matrices by putting a 1 in the j,kth position if the jth taxon is found in the kth bipartition, for each bipartition, and a 0 otherwise. We then accumulate these matrices into a supermatrix which contains union(species_i) rows and sum(m_i) columns. We put a "?" into the supermatrix if a species is missing from the tree. Finally, our MRP supertree is the MP tree estimated for this binary supermatrix.

To see how well this function does, let's first generate a random reference tree with 50 species:

> true.tree<-rtree(n=50,br=NULL)

Now let's randomly prune 40-60% of the species from this tree ten times:

> trees<-list()
> class(trees)<-"multiPhylo"
> for(i in 1:10)
trees[[i]]<- drop.tip(true.tree,sample(true.tree$tip.label,
   size=(20+round(runif(1)*10))))

Now let's load our MRP supertree code from source and run it:

> source("mrp.supertree.R")
> super.tree<-mrp.supertree(trees)
The MRP supertree, optimized via pratchet(),
has a parsimony score of 247 (minimum 247)

Now we can compare the estimated supertree to the true tree using dist.topo() in {ape}:

> dist.topo(true.tree,super.tree)
[1] 4

In this case, it is close. If we have high overlap this will generally be the case. For instance if we only drop 10-20 species from the each tree at random:

> for(i in 1:10)
trees[[i]]<- drop.tip(true.tree,sample(true.tree$tip.label,
   size=(10+round(runif(1)*10))))
> super.tree<-mrp.supertree(trees)
The MRP supertree, optimized via pratchet(),
has a parsimony score of 315 (minimum 315)
> dist.topo(true.tree,super.tree)
[1] 0

whereas if we drop 30-40 at random:

> for(i in 1:10)
trees[[i]]<-drop.tip(true.tree,
   sample(true.tree$tip.label,
size=(30+round(runif(1)*10))))
> super.tree<-mrp.supertree(trees)
The MRP supertree, optimized via pratchet(),
has a parsimony score of 208 (minimum 141)
> dist.topo(true.tree,super.tree)
[1] 82

A few notes on this implementation. Firstly, I set pratchet() to return all equally parsimonious trees (if several are found, mrp.supertree() will return a "multiPhylo" object) - but this is not the same as a guarantee to find all most parsimonious trees. This could only be guaranteed with an exhaustive search. Secondly, pratchet() does not always find the MP tree - so you might run mrp.supertree() multiple times if the parsimony score is greater than the theoretically minimal score.

Feedback welcome.