In McTavishLab/physcraperex: Physcraper Workflow Examples

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "../data-raw/malvaceae-images/"
)
knitr::opts_knit$set(
  global.par = TRUE
  )
options(rmarkdown.html_vignette.check_title = FALSE)

library(webshot)
getwd()
source("../R/plots.R")
library(ape)

knitr::include_graphics("malvaceae-plant.jpg")

I. Finding a tree to update

With the Open Tree of Life website

Go to the Open Tree of Life website and use the "search for taxon" menu to look up the taxon Malvaceae.

This is how the genus Malvaceae is represented on the Open Tree of Life synthetic tree at the middle of year 2020:

# webshot::webshot("https://tree.opentreeoflife.org/opentree/argus/ottol@279960/Malvaceae", file = "malvaceae-tree-opentree.png", delay = 1.5, zoom = 1, vwidth = 1400, vheight=1200)

knitr::include_graphics("malvaceae-tree-opentree.png")

---

Let's go to the study curator of OToL and look up for uploaded studies for the Malvaceae

Studies matching the word 'malvaceae' on the curator database, at the middle of year 2020. Some of these studies are not actually about the Malvaceae, but other taxa that have the keyword malvaceae.

# webshot::webshot("https://tree.opentreeoflife.org/curator/?match=malvaceae", file = "malvaceae-curation-opentree.png", delay = 1, zoom = 1, vwidth = 1400)
knitr::include_graphics("malvaceae-curation-opentree.png")

---

Finding a tree to update using the R package rotl

Explain what a focal clade is.

There is a handy function that will search a taxon among the focal clades reported across trees.

studies <- rotl::studies_find_studies(property="ot:focalCladeOTTTaxonName", value="Malvaceae")

knitr::kable(studies, caption = "Studies with the genus *Malvaceae* as focal clade.")

None of these trees are in the synthetic tree. Still, the family does not appear as a politomy and has phylogenetic information in the Open Tree of Life synthetic tree. Where is this information coming from?

It is probably coming from studies with focal clades that might be either above or below the particular clade we are interested in!

You can get supporting trees with the Open Tree of Life API, and you can currently do it now using a hidden rotl function:

res <- rotl:::.tol_subtree(ott_id=279960)
# studies <- rotl:::studies_from_otl(res) # not supported anymore as of April 2021
# studies

Updating study pg_55

# webshot::webshot("https://tree.opentreeoflife.org/curator/study/view/pg_55", file = "pg_55-curation-opentree.png", delay = 1, zoom = 1, vwidth = 1400)
knitr::include_graphics("pg_55-curation-opentree.png")

pg_55_meta = rotl::get_study_meta(study_id = "pg_55")
ls(pg_55_meta$nexml)
pg_55_meta$nexml$treesById
unlist(pg_55_meta$nexml$treesById)[1]
# rotl::taxonomy_taxon_info(pg_55_meta$nexml$"^ot:focalClade")

---

The focal clade of study is r unlist(unname(pg_55_meta$nexml$"^ot:focalCladeOTTTaxonName")). Is this taxon in the synthetic tree? r rotl::is_in_tree(pg_55_meta$nexml$"^ot:focalClade")

Let's get the original tree and plot it here:

par(mar=c(0,1,1,0), xpd=NA)  # adjust margins for all plots after this

original_tree <- rotl::get_study_tree(study_id = "pg_55", tree_id = "tree5864")
ape::plot.phylo(ape::ladderize(original_tree), type = "phylogram", cex = 0.3, label.offset = 1, edge.width = 0.5)
mtext("Original tree - original labels", side = 3, cex = 0.6)

---

Now, let's look at some properties of the tree:

ape::Ntip(original_tree)  # gets the number of tips
ape::is.rooted(original_tree)  # check that it is rooted
ape::is.binary(original_tree)  # check that it is fully resolved
datelife::phylo_has_brlen(original_tree)  # checks that it has branch lengths

The tree has r ape::Ntip(original_tree) tips, is rooted, has no branch lengths and is not fully resolved, as you probably already noticed. Also, labels correspond to the labels reported on the original study here. Other labels are available to use as tip labels. For example, you can plot the tree using other labels, like the unified taxonomic names, or the taxonomic ids.

II. Getting the underlying alignment

Using physcraper and the arguments -tb and -no_est

This will download the alignment directly from TreeBASE:

physcraper_run.py -s pg_55 -t tree5864 -tb -no_est -o data/pg_55

The alignment will be saved as pg_55tree5864.aln and physcraper_pg_55tree5864.fas in the inputs folder.

Downloading the alignment directly from a repository

We know that the alignments are in the TreeBASE repo at https://treebase.org/treebase-web/search/study/matrices.html?id=1361

Let's save them to a file named TB2:Tr240.nex in the data-raw/alignments folder:

wget -O data-raw/alignments/TB2:Tr240.nex "http://purl.org/phylo/treebase/phylows/tree/TB2:Tr240?format=nexus"

III. Starting a physcraper run

We ran this example using the command:

physcraper_run.py -s pg_55 -t tree5864 -a data-raw/alignments/TB2:Tr240.nex -as nexus -o data/pg_55

IV. Reading the physcraper results

The physcraper tag

pg_55tree5864_ndhf

Input files

Physcraper rewrites input files for a couple reasons: reproducibility, taxon name matching, and taxon reconciliation. It writes the config file down if none was provided. Everything is written into the inputs_tag folder.

Run files

Files in here are also automatically written down by physcraper.

blast runs, alignments, raxml trees, bootstrap

The trees are reconstructed using RAxML, with tip labels corresponding to local ids (e.g., otu42009, otuPS1) and not taxon names (e.g., Ceiba), nor taxonomic ids (e.g., ott or ncbi). Branch lengths are proportional to relative substitution rates. The RAxML tree with taxon names as tip labels is saved on the outputs_tag folder.

updated_tree_otus <- ape::read.tree(file = "../data/pg_55/run_pg_55tree5864_ndhf/RAxML_bestTree.2020-06-09")
ape::plot.phylo(ape::ladderize(updated_tree_otus), type = "phylogram", cex = 0.25, label.offset = 0.001, edge.width = 0.5)
ape::add.scale.bar(cex = 0.3, font = 1, col = "black")
mtext("Updated tree - otu tags as labels", side = 3)

Output files

A number of files is automatically written down by physcraper.

A nexson tree with all types of tip labels is saved in here. From this tree, a tree with any kind of label can be produced. By default, the updtaed tree with taxon names as tip labels is saved in the output_tag folder as updated_taxonname.tre.

updated_tree_taxonname <- ape::read.tree(
  file = "../data/pg_55/outputs_pg_55tree5864_ndhf/updated_taxonname.tre")
updated_tree_taxonname$tip.label.original <- updated_tree_taxonname$tip.label

par(mar=c(1.8,1.7,1.5,1.7), xpd=NA)  # adjust margins for all plots after this

ape::plot.phylo(ape::ladderize(updated_tree_taxonname), type = "phylogram", cex = 0.2, label.offset = 0.001, edge.width = 0.5)
ape::add.scale.bar(cex = 0.3, font = 1, col = "black")
mtext("Updated tree - Taxon names as labels", side = 3)

IV. Analyzing the physcraper results

Compare the original tree with the pruned updated tree

original_tree_otus <- ape::read.tree(file = "../data/pg_55/inputs_pg_55tree5864_ndhf/physcraper_pg_55tree5864_ndhf.tre")
updated_tree_otus_pruned <- ape::read.tree(
  file = "../data/pg_55/pruned_updated.tre"
)
# This tree only has 31 tips because it is the pruned updated tree.

Lets get indices from the new tips added to the updated tree and see which taxa they correspond to:

new_tips <- grep("_otuPS[0-9]*", updated_tree_taxonname$tip.label.original)
updated_tree_taxonname$tip.label.original[new_tips]

There are r length(new_tips) new tips in the updated tree.

Now plot the trees face to face.

We can first prune the updated tree of the new tips, so it is a straight forward comparison:

cotree <-  phytools::cophylo(original_tree_otus, updated_tree_otus_pruned, rotate.multi =TRUE)
phytools::plot.cophylo(x = cotree, fsize = 0.5, lwd = 0.5, mar=c(.1,.2,2,.3), ylim=c(-.1,1), scale.bar=c(0, .05))
title("Original tr.   Pruned updated tr.", cex = 0.1)

But it is more interesting to plot it with all the new tips, so we see exactly where the new things are:

original_tree_taxonname <- ape::read.tree(file = "../data/pg_55/inputs_pg_55tree5864_ndhf/taxonname.tre")
cotree2 <-  phytools::cophylo(
  original_tree_taxonname,
  updated_tree_taxonname,
  rotate.multi =TRUE)
phytools::plot.cophylo(
  x = cotree2,
  fsize = 0.3,
  lwd = 0.4,
  mar=c(.1,.1,2,.5),
  ylim=c(-.1,1),
  scale.bar=c(0, .05),
  # link.type="curved",
  link.lwd=3,
  link.lty="solid",
  link.col=phytools::make.transparent("#8B008B",0.5))
title("Original tree      Updated tree", cex = 0.1)

We can also plot the updated tree against the synthetic subtree of Malvaceae, to visualize how it updates our current knowledge of the phylogenetic relationships within the family:

tolsubtree <- rotl::tol_subtree(ott_id = 279960)

ape::Ntip(tolsubtree)

head(tolsubtree$tip.label)
head(updated_tree_taxonname$tip.label)

ape::Ntip tells us that there are almost 6k taxa in the Malvaceae OpenTree synthetic tree.

Tip labels in both trees have taxon identifier tags that will make name matching impossible. Lets remove them:

updated_tree_taxonname$tip.label <- gsub("_otuPS[0-9]*", "", updated_tree_taxonname$tip.label.original)
updated_tree_taxonname$tip.label <- gsub("_otu[0-9]*", "", updated_tree_taxonname$tip.label)

tolsubtree$tip.label.original <- tolsubtree$tip.label
tolsubtree$tip.label <- gsub("_ott[0-9]*", "", tolsubtree$tip.label)

Eliminate ott and otu tags before comparison, otherwise nothing matches

# cotree3 <-  phytools::cophylo(tolsubtree, updated_tree_taxonname, rotate.multi =TRUE)
# phytools::plot.cophylo(x = cotree3, fsize = 0.2, lwd = 0.5, mar=c(.1,.1,2,.3), ylim=c(-.1,1), scale.bar=c(0, .05))
# title("Tol subtree      Updated tree", cex = 0.1)

The whole Malvaceae tree might be too large for a cotree plot. Lets try something else.

First, get indices of new tips from updated tree present in the Malvaceae synthetic tree:

all_updated_tree_in_tol <- match(updated_tree_taxonname$tip.label, tolsubtree$tip.label)
new_tips_in_tol <- match(updated_tree_taxonname$tip.label[new_tips], tolsubtree$tip.label)
# one tip seems to be absent frm the Malvaceae subtree
updated_tree_taxonname$tip.label[47]
# going to OpenTree it looks like it atually corresponds to Firmiana simplex, not sure why it was not standardized, probably bc it was in NCBI
# I'm changing it by hand:
updated_tree_taxonname$tip.label[47] <- "Firmiana_simplex"
# And running match again
new_in_tol <- match(updated_tree_taxonname$tip.label, tolsubtree$tip.label)

Then, use the `phytools package to create a data object that will contain info from our tolsubtree tips and edges:

tolsubtree <- ape::ladderize(tolsubtree)

tip_colors <- rep("gray", ape::Ntip(tolsubtree))
tip_colors[all_updated_tree_in_tol] <- "black"
tip_colors[new_tips_in_tol] <- "red"

tip_values <- rep(0, ape::Ntip(tolsubtree))
tip_values[tip_colors == "black"] <- 1
tip_values[tip_colors == "red"] <- 100
names(tip_values) <- tolsubtree$tip.label
head(tip_values)
all_values <- sum_tips(tolsubtree, tip_values)

all_cols <- rep("gray", length(all_values))
all_values[all_values >= 100]
all_values[all_values >= 1]

all_cols[all_values >= 1] <- "black"
all_cols[all_values >= 100] <- "red"
edge_colors <- all_cols[tolsubtree$edge[,2]]
length(edge_colors) == nrow(tolsubtree$edge)
# tip_cols <- all_cols[seq(ape::Ntip(tolsubtree))]

namess <- rep("Synth tips", length(edge_colors))
namess[which(edge_colors == "black")] <- "Original tips"
namess[which(edge_colors == "red")] <- "Updated tips"

tolsubtree$maps <- as.list(rep(1, length(edge_colors)))
for (i in seq(length(edge_colors))){
  names(tolsubtree$maps[[i]]) <- namess[i]
}
# fam_tree_brlen3$maps[1]
# str(fam_tree_brlen3$maps)
cols <- setNames(c("red","black", "gray"), c("Updated tips", "Original tips", "Synth tips"))
tolsubtree_wbl <- ape::compute.brlen(tolsubtree)
tolsubtree_wbl$maps <- tolsubtree$maps

par(fg="transparent")
phytools::plotSimmap(tolsubtree_wbl,type="fan",part=0.95, ftype="i", colors = cols, lwd = 0.1, fsize = 0.05) #lwd = 0.5, fsize=0.1
lastPP<-get("last_plot.phylo",env=.PlotPhyloEnv)
par(fg="black")
phytools::add.simmap.legend(colors=cols,vertical=FALSE,x=-0.9, y=-1.2,prompt=FALSE, fsize = 0.5)
tt<-gsub("_"," ",tolsubtree_wbl$tip.label)
#tt[tip_colors == "gray"] <- ""
text(lastPP$xx[1:length(tt)],lastPP$yy[1:length(tt)],
  tt,cex=0.2,col=tip_colors,pos=4,offset=0.1,font=3)

References

Trick for the cophylo titles and margins from https://cran.r-project.org/web/packages/phangorn/vignettes/IntertwiningTreesAndNetworks.html

Trick for colord tip labels: http://blog.phytools.org/2013/09/coloring-tip-labels-by-character-state.html

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