parsimony and maximum likelihood tree comparison in R

I don't use R to build trees, I use established subprocessing stripts in Python and ETE3. So I've not run this code, but I do understand the critical parts of the code.

The hard bit is bootstrapping a maximum likelihood tree, first is the ML model 1.

fitGTR <-optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE, + rearrangement = "NNI", control = pml.control(trace = 0))

The above model is perfect, general time reversible model with a gamma correction and an estimate of the invariant sites. This is the most important part of tree estimation. The NNI for tree rearrangements is a good one. So THIS is the model you want. You can emperical test this, but these days we don't bother.

The information given for bootstrapping,

2.

BStree <- bootstrap.pml(x, bs = 500, trees = TRUE, multicore = FALSE,
  mc.cores = NULL, ...)

The critical parameter is bs=500. By default it uses bs=100 which I don't think is enough. I would set multicore = TRUE to speed the calculation (it will speed it up alot) and I assume there must be a setting for the number of cores (?).

I assume you can then stack all the parameters from 1/ into code 2. So in other words the model you used in part 1. (GTR PINVAR and Gammma) is same model used in the bootstrapping.

This will produce a bootstrap file, probably of 500 bootstrapped ML trees, so you need to produce a consensus and represent the values on the optimal ML tree.

Thus you need to overlay the bootstrap onto the tree that is output in 1. This says fitGTR and I suspect that will be a tree file with a likelihood, you might need to set a tree output and that would be 'tree' in the code below. 3.

plotBS(tree, BStrees, type = "myroot", bs.col = "black",
  bs.adj = NULL, p = 50, frame = "none", ...)

The type = myroot is your defined outgroup for the tree. This should be your nicely bootstrapped tree.

There was also this code in the help manual and I'm not sure what this is doing. I suspect it is just producing bootstrapped datasets, so its not needed here.

bootstrap.phyDat(x, FUN, bs = 100, multicore = FALSE,
      mc.cores = NULL, jumble = TRUE, ...)

To do the same with parsimony then just replace pml for parsimony. Comparing 2 bootstrap trees is a node is robust when it scores >80%