phylogeny.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Phylogenetic reconstruction using R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this course we will take a look at the phylogenetics of the bacterial representative OTU sequences obtained in the data processing part.\n",
    "\n",
    "* First we will import and filter the sequences.\n",
    "\n",
    "* After making the corresponding alignments we will perform different approaches to phylogenetic inference.\n",
    "\n",
    "* Finally, we will visualize and interpret the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the R packages we will need. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(msa)\n",
    "library(ggtree)\n",
    "library(phangorn)\n",
    "library(tidyverse)\n",
    "set.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequence import  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the taxonomy file corresponding to the representative sequences using the `read_delim` function with the parameter `delim = \"\\t\"` in order to fit the tab delimited format of the file. Assign the output of this function to a variable called `taxonomy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taxonomy <- read_delim(\n",
    "    \"otu_info/BV5.trim.contigs.good.renamed.unique.pick.dgc.unique_list.0.03.abund.0.03.cons.taxonomy\", \n",
    "    delim = \"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: type taxonomy in the code cell below and run it (by selecting the cell and pressing shift+enter) to inspect the object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: how many columns are there in the object? How many rows? What do they represent?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*type your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the fasta file from the representative sequences using the the function `readDNAStringSet`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bv5t <- readDNAStringSet(\n",
    "    \"otu_info/BV5.trim.contigs.good.renamed.unique.pick.dgc.unique_list.0.03.abund.0.03.rep.fasta\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: inspect what the function `readDNAStringSet` does by adding `?` before the name of the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## type here your answer "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: inspect the created object with the DNA sequences, called `bv5t`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## type your answer here "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Change fasta names to taxonomy\n",
    "Replace the names of the imported fasta object with the imported taxonomy. Since they are in the same order this can be simply done by replacing the `names` attribute from the `bv5t` fasta object with the `Taxonomy` column from the `taxonomy` dataframe. To prevent problems further on we also rename the duplicated taxonomy descriptions using the `make.unique` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names(bv5t) <- make.unique(taxonomy$Taxonomy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take only the sequences annotated as bacterial, that is, those that start with \"k__Bacteria(100)\" from the `bv5t` fasta object and name it as `bv5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bv5 <- bv5t[startsWith(names(bv5t), \"k__Bacteria(100)\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove bacteria from the name to be able to identify the phyla more easily. We do this by substracting the 21 first characters using the `substring` function with the parameter `first = 21`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names(bv5) <- substring(names(bv5), first = 21)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: how many OTU's are deleted?\n",
    "\n",
    "*hint: you may need to inspect both `bv5` and `bv5t` objects*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](figs/align.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Evolution and the true multiple sequence alignment.** The top sequence evolves into the bottom sequence via the deletion of the substring GGTG, the substitution of a T for a C, and the insertion of a T. This corresponds to the pairwise alignment on the right. Note that two letters are placed in the same column only when they have a common history. Thus, the substring GGTG in the top string is above dashes in the bottom string, and indicates that deletion event. Similarly, the red T is above the blue C, to indicate that they have a common history. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Align sequences using `msa` function (this may take a bit of time)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bv5_aligned = msa(bv5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: inspect the function `msa`. Which alignment algorithm was used? How would you use a different algorithm?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: rerun msa using a different algorithm, and name it as a different variable. How does it compare to the previous method?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: inspect the object containing the alignment. What information can you get from it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Phylogenetic tree reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse to phangorn format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alignment <- as.phyDat(bv5_aligned)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: inspect the new alignment object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distance based methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Distance methods try to fit a tree to a matrix of pairwise genetic distances.\n",
    "\n",
    "Advantage:\n",
    "- Fast speed of calculation\n",
    "\n",
    "Disadvantages:\n",
    "- Often not accurate (particularly in highly divergent sequences)\n",
    "- Loss of evolutionary information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Cluster analysis**\n",
    "**UPGMA** (unweighted pair group method with arithmetic means)\n",
    "![](figs/upgma1.png)\n",
    "![](figs/upgma2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Minimum evolution**\n",
    "**Neighbour Joining**\n",
    "![](figs/nj.svg.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct a rooted tree using UPGMA and an unrooted one using NJ (Neighbour-joining)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm <- dist.ml(alignment)\n",
    "treeUPGMA  <- upgma(dm)\n",
    "treeNJ  <- NJ(dm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the tree object `treeNJ` and `treeUPGM`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "treeUPGMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "treeNJ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: what information can you gather from the objects? what data is contained in it?\n",
    "\n",
    "*hint: you may need to use the str() function on the objects to see its data.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the tree with the tip names in small font using the `ggtree` package. The `geom_tip` adds the layer of leaf tips."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=10,repr.plot.res=500)\n",
    "ggtree(treeUPGMA) + geom_tiplab(size=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `layout` parameter in ggtree defines the shape of the tree. Other options are: \"slanted\", \"fan\", \"circular\", \"equal_angle\", \"daylight\". \n",
    "\n",
    "**Exercise**: plot the same rooted tree with the parameter \"fan\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=4,repr.plot.res=500)\n",
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: plot the NJ tree in the \"equal_angle\" layout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Character based methods\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of a distance matrix, the difference among the nucleotides is computed and fitted to a substitution model. A tree that explains this model is fitted to the data. \n",
    "\n",
    "Advantages: \n",
    "- Uses all the information available in sequences at each homologous site\n",
    "- Output measure of how good is the tree\n",
    "    - Possible to compare with other tree inferences for the same data\n",
    "\n",
    "Disadvantage:\n",
    "- Computationally demanding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum parsimony"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">The most believable or parsimonious phylogenetic tree will be the tree that invokes the smallest number of evolutionary changes during the divergence of the sequences it represents.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the parsimony score for the two trees computed so far. The output integer for the `parsimony` function, given a tree and its corresponding alignment, is the parsimony score. The parsimony score is defined as the number of changes which are at least necessary to describe the data for a given tree. To get more information you can inspect the help page for the `parsimony` function by running `?parsimony` in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parsimony(treeUPGMA, alignment)\n",
    "parsimony(treeNJ, alignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: which tree is more believable?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*type here your answer*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the function `optim.parsimony` which performs tree rearrangements to find trees with a lower parsimony score. The tree rearrangement implemented are\n",
    "nearest-neighbor interchanges (NNI) and subtree pruning and regrafting (SPR). See the image below for more information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](figs/nni.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "treePars  <- optim.parsimony(treeUPGMA, alignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the optimization but using the parsimony ratchet algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "treeRatchet  <- pratchet(alignment, trace = 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the parsimony score of both methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parsimony(c(treePars, treeRatchet), alignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: which optimization performs better? and which takes longer?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*type here your answer*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">The basic concept of likelihood is relatively simple to comprehend: given some data D (in our case, nucleotide or amino acid sequences), under a model of evolution, M (which is explicitly defined and describes the mutation process from one base to another), the likelihood of a set of parameters, θ (tree topology, tree branch lengths, substitution model parameters), corresponds to the probability of obtaining D under the model M with parameters θ. The maximum likelihood estimates of the parameter values included in θ correspond to the set of values that maximize this probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the maximum likelihood tree. The `pml` model contains the tree and different parameters of the model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit = pml(treeNJ, data = alignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: explore the newly created fit object. What information does it tell us?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform branch length optimization using Jukes-Cantor model (all base changes equally likely) with the function `optim.pml` and default values. First look at the help page of `optim.pml`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitJC  <- optim.pml(fit, TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: obtain the log likelihood of the optimizied tree using the function `logLik` in the object `fitJC`. What does \"df\" mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change fitting to GTR model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitGTR <- update(fit, k = 4, inv = 0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: explore the `fitGTR`object just created. How does the log likelihood compared to the previous models used?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimize all parameters to this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitGTR <- optim.pml(fitGTR, model=\"GTR\", optInv=TRUE, optGamma=TRUE,\n",
    "                    rearrangement = \"NNI\", control = pml.control(trace = 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: explore the object after optimization, how did it change?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: inspect the object fitGTR with the `attributes` and `str` function, where's the tree?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=9,repr.plot.res=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggtree(fitGTR$tree) + geom_tiplab(size=0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trim the name of the leave to the first word (the phylum). This is done using a regex expression that removes everything after the first `(` (substitutes it by nothing)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitGTRshortnames <- fitGTR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitGTRshortnames$tree$tip.label <- gsub('\\\\(.*$',\"\",fitGTRshortnames$tree$tip.label)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: explore the new names printing the object we just modified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type here your answer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the grouped OTUs, using the `groupOTU` function on some of the most abundant phyla before plotting the tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grouped_tree <- groupOTU(fitGTRshortnames$tree, c(\"Proteobacteria\", \"Actinobacteria\", \"Chloroflexi\", \"Acidobacteria\"))\n",
    "\n",
    "ggtree(grouped_tree, aes(color=group)) + geom_tiplab(size =1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: what would you expect the unclassified OTU's to be?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*type here your answer*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter clades by phylogeny."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "betaproteo <- fitGTR$tree$tip.label[startsWith(fitGTR$tree$tip.label,\"Proteobacteria(100);c__Bet\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the MRCA for nodes that start with a particular string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mrca_beta <- getMRCA(fitGTR$tree,betaproteo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise** what does MRCA mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type here your answer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the clade from the tree to make a new tree that only contains that clade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_beta <- extract.clade(fitGTR$tree,mrca_beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Highlight *betaproteobacteria* clade in green."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggtree(fitGTR$tree) + geom_hilight(node=mrca_beta, fill=\"darkgreen\", alpha=.6) ## change the node number to fit the actual clade"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: highlight also the *alphaproteobacteria* clade in blue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## type here your answer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](figs/models.png) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare nested models for JC and GTR using likelihood ratio statistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anova(fitJC, fitGTR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which nucleotide substitution model is more likely?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: read the help page for the `modelTest`function. What do the parameters multicore and m.core do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform a model test (this may take a bit of time)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = modelTest(alignment,fitGTR$tree,mc.cores=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: inspecfitGTR$treeults for the model test, which nucleotide substitution model is more likely?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OTU table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the OTU table file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "otu_table <- read_delim(\n",
    "    \"otu_info/BV5.trim.contigs.good.renamed.unique.pick.dgc.unique_list.0.03.abund.0.03.rep.count_table\",\n",
    "    \"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add first column of the table to rownames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "otu_table <- column_to_rownames(otu_table, var = \"Representative_Sequence\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: inspect the imported object, `otu_table`. What do the columns and rows mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: using the function `column_to_rownames` applied to the object `otu_table`, move the column first column called \"Representative_Sequence\" to the rownames. Name the result `otu_table` again. You may need to use the help page of `column_to_rownames` to find which parameter to add. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change to original sequence names back from the taxonomy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bv5nt <- readDNAStringSet(\"otu_info/BV5.trim.contigs.good.renamed.unique.pick.dgc.unique_list.0.03.abund.0.03.rep.fasta\")\n",
    "bv5nt <- bv5nt[startsWith(names(bv5t), \"k__Bacteria(100)\")]\n",
    "aln <- alignment\n",
    "names(aln) <- names(bv5nt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: how's the new object containing the original OTU names called?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a tree again using ML and the GTR substitution model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "treeNJ1  <- NJ(dist.ml(aln)) # constructs nj tree\n",
    "\n",
    "fit1 = pml(treeNJ1, data=aln) # fits tree to alignment\n",
    "\n",
    "fitGTR1 <- update(fit1, k=4, inv=0.2) \n",
    "\n",
    "fitGTR1 <- optim.pml(fitGTR1, model=\"GTR\", optInv=TRUE, optGamma=TRUE,\n",
    "                    rearrangement = \"NNI\", control = pml.control(trace = 0)) # fits tree to GTR substitution model\n",
    "\n",
    "fitGTR1$tree$tip.label <- gsub('\\t.*$',\"\",fitGTR1$tree$tip.label) # trims the names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: store tree figure with highlighted *alpha* and *proteobacteria clades*, this time using the newly created tree with annotated clades as a variable simply named `p`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heatmap abundances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the OTU table for the 12 samples as a heatmap of their abundance. Notice the log transformation on the OTU abundance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gheatmap(p,log(otu_table), \n",
    "         width =1,colnames_angle = 90, colnames_offset_y = -12) + \n",
    "scale_fill_viridis_c(option=\"A\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Add tip labels to the tree.\n",
    "\n",
    "*tip: you can adjust font with the `size` parameter*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: are the taxa abundances consistent along different types of sample? what are the more striking differences? \n",
    "\n",
    "How are *Pseudomonas*, *Corynebacterium* and *Acetobacter* distributed across the different samples?\n",
    "\n",
    "*tip: to make it easier to find the clades in the tree you can subset as previously in the `fitGTR` objects. Then subset the tree in fitGTR1 for the heatmap, as shown below for the betaproteobacteria.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=5,repr.plot.res=500)\n",
    "betaproteo <- fitGTR$tree$tip.label[startsWith(fitGTR$tree$tip.label,\"Proteobacteria(100);c__Bet\")]\n",
    "mrca_beta <- getMRCA(fitGTR$tree,betaproteo)\n",
    "tree_beta <- extract.clade(fitGTR1$tree,mrca_beta)\n",
    "p <- ggtree(tree_beta,branch.length = \"none\")\n",
    "gheatmap(p,log(otu_table), \n",
    "         width =1,colnames_angle = 90, colnames_offset_y = -2) + \n",
    "scale_fill_viridis_c(option=\"A\")+ ggtree::vexpand(.1, -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*type here your answer*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}