In McTavishLab/physcraperex: Physcraper Workflow Examples

source("../R/plots.R")
updated_sumtree <- ape::read.tree(file = "../data/pg_2827_tree6577/outputs_pg_2827tree6577/physcraper_pg_2827tree6577.tre")
updated_sumtree_drop <- ape::drop.tip(updated_sumtree, c("otu420083", "otu420099"))
updated_sumtree_drop <- ape::ladderize(updated_sumtree_drop)
internal_edges <- updated_sumtree_drop$edge[,2] > ape::Ntip(updated_sumtree_drop)
tip_edges <- updated_sumtree_drop$edge[,2] <= ape::Ntip(updated_sumtree_drop)
boots <- as.numeric(updated_sumtree_drop$node.label)
sum(boots>0.8)
# internal edges that are very small:
sum(updated_sumtree_drop$edge.length[internal_edges]<0.00001) # 30
# tip edges that are very small:
sum(updated_sumtree_drop$edge.length[tip_edges]<0.00001) # 77
updated_sumtree_drop$Nnode
ie <- updated_sumtree_drop$edge.length[internal_edges]
ie <- c(1, ie) # adding a root edge length
bb <- which(ie<0.00001 & boots>0.75)
ie[bb]
boots[bb]
updated_sumtree_drop$edge.length[bb+ape::Ntip(updated_sumtree_drop)]
ss <- sapply(bb+ape::Ntip(updated_sumtree_drop)-1, function(x)
              phytools::getDescendants(tree=updated_sumtree_drop, node=x))

Visualizing new tips and taxa

An easy way to visualize the new tips added into the updated tree compared to the original tree is coloring the tips and branches from those new tips with a highlighting color. We will use the function plot_branches defined in this package to do so.

The function takes as argument the paths to both the tree and the OTU info files:

# mytreefile = '../data/pg_2827_tree6577/run_pg_2827tree6577_run4/RAxML_bestTree.2020-07-31'
mytreefile = '../data/pg_2827_tree6577/outputs_pg_2827tree6577/physcraper_pg_2827tree6577.tre'
myotuinfofile = '../data/pg_2827_tree6577/outputs_pg_2827tree6577/otu_info_pg_2827tree6577.csv'

And uses phytools::plotBranchbyTrait, so it takes all arguments from ape::plot.phylo, too:

plot_branches(treefile = mytreefile, otufile = myotuinfofile, type = "phylogram", cex = 0.25, edge.width = 0.5, tip_label = "taxon", color = "goldenrod", edge_length = FALSE)
mtext(expression(Updated~italic(Ilex) ~ plain("gene tree - 4 BLAST cycles")), side= 3, line=-1, cex = 0.5)
# ape::nodelabels(node=bb+ape::Ntip(updated_sumtree_drop)-1, pch =8, cex = 0.1, col = "green")  # uncomment this if you want to show tiny-tiny branches that have a large bootstrap support!

The green nodes highlight branches with a length < 0.00001 and a bootstrap support > 0.75

Let's get the summary of new tips:

summ <- summarize(treefile = mytreefile, otufile = myotuinfofile)
updated_sumtree_drop <- ape::ladderize(updated_sumtree_drop)
internal_edges <- updated_sumtree_drop$edge[,2] > ape::Ntip(updated_sumtree_drop)
tip_edges <- updated_sumtree_drop$edge[,2] <= ape::Ntip(updated_sumtree_drop)
boots <- as.numeric(updated_sumtree_drop$node.label)
sum(boots>0.8)
# internal edges that are very small:
sum(updated_sumtree_drop$edge.length[internal_edges]<0.00001) # 30
# tip edges that are very small:
sum(updated_sumtree_drop$edge.length[tip_edges]<0.00001) # 77
updated_sumtree_drop$Nnode
ie <- updated_sumtree_drop$edge.length[internal_edges]
ie <- c(1, ie) # adding a root edge length
sum(boots<0.75)/ape::Ntip(updated_sumtree_drop)
bb <- which(ie<0.00001 & boots>0.75)
ie[bb]
boots[bb]
updated_sumtree_drop$edge.length[bb+ape::Ntip(updated_sumtree_drop)]
ss <- sapply(bb+ape::Ntip(updated_sumtree_drop)-1, function(x)
              phytools::getDescendants(tree=updated_sumtree_drop, node=x))

hist(updated_sumtree_drop$edge.length, breaks=1400, col = "skyblue3",
     border = "skyblue3", xlab= "Branch Length (substitution rate)",
     main = "Updated Gottlieb2005 tree \n branch length frequency distribution")

This branch length distribution excludes the outgroup branches, it contains branches from the ingroup only.

The updated tree is r round(updated_sumtree_drop$Nnode/ape::Ntip(updated_sumtree_drop), digits=2)*100% resolved, with r sum(updated_sumtree_drop$edge.length<0.01) branches that are <0.01 and r sum(updated_sumtree_drop$edge.length<0.00001) are < 0.00001, effectively negligible. The longest branch is r max(updated_sumtree_drop$edge.length) and the smallest branch is r min(updated_sumtree_drop$edge.length)

McTavishLab/physcraperex documentation built on April 10, 2021, 12:02 a.m.