How do I estimate branch length of a maximum parsimony phylogeny in R?

Question

I am trying to estimate a phylogenetic tree in R with ape and phangorn. With pratchet() and optim.parsimony(), I can estimate a phylogenetic tree, but branch lengths are not estimated. Is there a way to estimate a phylogeny with branch lengths with the following made-up data:

x <- data.frame("L1" = c(0,0,1,0,0,0,1,1,1,1,0,0,0,0,1), "L2" = c(0,1,1,0,1,1,0,1,1,0,1,1,0,1,1),"L3" = c(1,1,1,0,0,0,1,1,1,0,0,0,1,1,1), "L4"=c(0,0,1,1,0,0,0,0,1,1,0,0,0,0,1)) 

Answer 1

@Nathan S. Watson-Haigh's answer does give branch lengths, but NOT parsimony branch lengths, which may be what the OP wanted. nnls.phylo is giving distances calculated from a pairwise distance matrix. To get the parsimony branch lengths (i.e. number of changes mapped on to each branch), score the tree against the dataset using acctran:

treeRatchet <- pratchet(primates)
treeRatchet.BL <- acctran(treeRatchet,primates)
plot(treeRatchet.BL,main='Phylogeny (Parsimony branch lengths)')

Note that the branch lengths in the parsimony phylogram are 'chunkier' than those in the nnls.phylo phylogram -- that's because

1. they are all integers, and
2. they aren't corrected with reference to a model of substitution that considers multiple hits or base composition.

Answer 2

Terminology

First, a couple of corrections to your use of terminology. A tree must have branch lengths associated with it to be called a phylogenetic tree as it represents both the relationships between OTU/taxa/leaves AND the distances between them (i.e. the branch lengths. A tree where branch lengths have no meaning is called a cladograms or a dendrogram. You can tell if a tree is a cladogram since all the OUTs/taxa/leaves end at the same position and internal branch lengths are usually (not always) of equal length. Whereas a phylogeny the OTUs/taxa/leaves will typically look "staggered" and internal branches will be of different lengths.

Therefore, If the tree you have estimated, truly doesn't have branch lengths estimated and associated with it, it is actually a cladogram.

Example Derived from phangorn Vignette

library(phangorn)

# Get example data
fdir <- system.file("extdata/trees", package = "phangorn")
primates <- read.phyDat(file.path(fdir, "primates.dna"), format = "phylip")

# Calculate pairwise distances between OTUs using Maximum Likelihood
dm <- dist.ml(primates)

# Use the distances to cluster the OTUs using UPGMA algorithm
# This creates a dendrogram
treeUPGMA <- upgma(dm)
plot(treeUPGMA, main="UPGMA")
dev.off()

# Search for different, similar trees to the initial UPGMA tree
# looking for a tree with the lowest parsimony score
# There are 2 search algorithms implemented:
#   optim.parsimony()
#   pratchet()
# These are still cladograms
treePars <- optim.parsimony(treeUPGMA, primates)
treeRatchet <- pratchet(primates, trace = 0)
plot(treePars, main="optim.parsimony()")
plot(treeRatchet, main="pratchet()")
dev.off()

# Calculate branch lengths using the most parsimonious tree identified by
# the pratchet() searchh algorithm
phylogeny <- nnls.phylo(treeRatchet, dm)
plot(phylogeny, main="Phylogeny")
dev.off()

Choices, Choices, Choices

As you can see from the above example, there are lots of steps to calculating a phylogenetic tree. At each of these steps there may be multiple algorithms to choose from. You choice will depend on the type of data you are analysing, the size of the data set, your available compute resources and the algorithms implemented and available through the tools you use.

Ultimately, you will have to know a fair bit about these algorithms to make an informed choice.