Search code examples
rphylogenyggtree

Unable to save/write tree data that plotBS() results


I have aligned sequences in class of phyDat. Here is a minimal example:

library(seqinr)
seqs <- DNAMultipleAlignment( c(a1 = "------ATGTTCATTAACCGCTGACTATTCTCAACCA",
                                a2 = "------ATGTTCACCGACCGCTGACTATTCTCTACAA",
                                a3 = "---GTGACCTTCATCAACCGATGATTATTCTCAA---",
                                a4 = "---TCAGTCGTCACCAGGCGTTG-CAGGACCCGAC--",
                                a5 = "ATGGGGGTCTTCCTCA-TCGCCGTCGCCGCGT-----"))

library(phangorn)
phyDat_seqs <- as.phyDat(seqs)

I followed the manual of phangorn package to construct a phylogenic tree. (Below, I wrote a simplified code. You don't need to check manual, you can run the code directly.)

dm <- dist.ml(phyDat_seqs)
treeNJ <- NJ(dm)

fit <- pml(treeNJ, data=phyDat_seqs)
fitGTR <- optim.pml(fit, model="GTR")

bs <- bootstrap.pml(fitGTR, bs=5, optNni=TRUE) 

And I can simply plot this tree with plotBS :

plotBS(fitGTR$tree, bs, type="p")

There is no problem so far. My question starts from this point.

I want to customize this tree by using ggtree package but I couldn't figure out how to do that. When I use ggtree(fitGTR$tree), it gives me a different tree.

plotBS somehow uses bs values to change tree branching. If you compare plot(fitGTR$tree, type="p") and plotBS(fitGTR$tree, bs, type="p") difference can easily be seen. ggtree() and plot() construct the same trees, but not the tree that I want to customize.

I've tried different approaches such as saving tree with write.tree or write.nexus functions and calling the file again with ggtree package functions but I can't save the plotBS version of the tree. How can I get an output of plotBS tree?

Thank you in advance.


Solution

  • Problem was occuring because I misread the documentation. I did a mistake in this line:

    plotBS(fitGTR$tree, bs, type="p")
    

    Tree must be constructed through midpoint() function. Therefore, when I changed the line as this:

    plotBS(midpoint(fitGTR$tree), bs, type="p")
    

    no problem occurs. I can use ggtree and I can save/write the tree data without a problem.

    Even if I can't figure out what midpoint() function does to tree actually (you can comment about this), problem is solved.