Closed marwa38 closed 1 year ago
See the documentation on Estimating phylogenetic trees with phangorn.
There is an additional step you need to perform (bootstrap.pml
) to get bootstrap confidence estimates.
Thanks again @benjjneb
It takes about 5 days to run bootstrap (bold) in the code as follows
seqs <- getSequences(seqtab.nochim2) #extract sequences from dada2 output
names(seqs) <- seqs # This propagates to the tip labels of the tree
alignment <- AlignSeqs(DNAStringSet(seqs), anchor=NA)
#change sequence alignment output into a phyDat structure
phang.align <- phyDat(as(alignment, "matrix"), type="DNA")
dm <- dist.ml(phang.align) #create distance matrix
#reconstruct unrooted tree
treeNJ <- NJ(dm) # Note, tip order != sequence order #perform neighbor joining
fit <- pml(treeNJ, data=phang.align) #internal maxiumum likelihood
fit <- pml(treeNJ, data=phang.align)
fitJC <- optim.pml(fit, TRUE)
**bs = bootstrap.pml(fitJC, bs=100, optNni=TRUE,
control = pml.control(trace = 0))**
I need to run this on fitGTR, is there a difference between bootstrapping fitJC vs fitGTR?
fitGTR <- update(fit, k=4, inv=0.2)
fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
rearrangement = "stochastic", control = pml.control(trace = 0))
Thanks very much
is there a difference between bootstrapping fitJC vs fitGTR?
I would presume yes.
For larger trees, the phangorn implemantations get slow enough to be unwieldy to unusuable for microbiome datasets. You may want to consider using an external program such as RaxML to do the tree-ing part of the analysis.
Thanks very much for your paper on the workflow I am trying to find bootstrap in my tree so that I can have it in my tree graph but I haven't been successful yet, your advice is appreciated about which of the below 24 list is considered the bootstrap? I created the tree based on the paper as follows:
a quick example of the tree I viewed using plot_tree
here the structure of the final tree created before merging it into the phyloseq object dismiss the $data because it takes a space out of the 24 list shown below I don't know which one of these is considered the bootstrap, if this is to be figured out, I think bootstrap would be easily integrated into the tree plot as I need to do that.
> sessionInfo()