Assignment1_2023.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SCIE3100/BINF7000 Assignment 1\n",
    "\n",
    "## Probability, motif discovery, and ancestral sequence reconstruction\n",
    "\n",
    "* **Due:** 2PM Friday 18/8/2023 (Discussion board contributions), 2PM Friday 1/9/2023 (Part A and B solutions)\n",
    "* **Revision:** 2023 v1\n",
    "* **Marks:** 20% of course\n",
    "\n",
    "### Objectives\n",
    "\n",
    "Below are a number of exercises that aim to guide you through issues and help you understand concepts related primarily to:\n",
    "* Probability (introduced in week 1)\n",
    "* Motif discovery (introduced in week 2)\n",
    "* Graphical models and phylogenetics (introduced in week 4)\n",
    "\n",
    "\n",
    "### Format \n",
    "\n",
    "This assignment consists of two parts, each containing a number of problems of the two types discussed below. You are expected to work you way through this over a number of weeks in practicals and in your own time, supported by topics covered in course materials and tutorials. You are also able to interact with tutors in practicals and your fellow class mates.\n",
    "\n",
    "The assignment has a series of “A” problems to which short responses are assessed. Some responses are fixed-format and automatically marked as either correct (pass) or not (fail); these can be attempted multiple times before a final submission. Other responses are evaluated by a marker, after the submission. Solutions to “A” problems are based on individual work (such as research and experimentation, involving programming, data processing/analysis and interpretation). **Note:** Automatic marking is used in some instances to ensure an appropriate answer is reached which will allow completion of subsequent questions.\n",
    "\n",
    "The assignment also has more open-ended “B” problems for which short text responses are submitted and assessed. Solutions to “B” problems involve mostly *individual*, but also *collaborative* work based on exchange/discussion amongst members of the class. That said, your submission *must* report on your *own* work. So, you may *not* post your answers to the discussion board. You may *not* use somebody's text in your submission(see more on this below). The exchange and discussion are recorded via the online discussion forum. To acknowledge the exchange, posts that you \n",
    "have made and used information from *must* be listed in the submission. (A field will be available.)\n",
    "\n",
    "Note: the content to complete Part B has not been covered at the time of this assignment's release. Don't worry - we will be discussing the necessary concepts at length in upcoming lecture videos and tutorials!\n",
    "\n",
    "### Marking\n",
    "\n",
    "The assignment is worth 20% of the course, marked out of 20 marks. \n",
    "\n",
    "Marks are awarded as per the schedule below.\n",
    "\n",
    "#### “A” problems (10 marks total)\n",
    "\n",
    "| Marks | Criteria |\n",
    "| ----- | -------- |\n",
    "| 0 – 10| In proportion to number of *correct* responses|\n",
    "\n",
    "Marks are given for <span style=\"color:red\">*fixed-format questions*</span>, <span style=\"color:green\">*short text questions*</span>, and <span style=\"color:blue\">*submitted portions of your code*</span>. Note that it is *not* sufficient to attempt a question to get a mark.\n",
    "\n",
    "#### “B” problems (sum of two parts, capped at 10 marks in total)\n",
    "\n",
    "| Marks | Criteria | \n",
    "| ----- | -------- |\n",
    "| 0 – 3\t| Responses are inaccurate or absent |\n",
    "| 4 – 5\t| Responses are incomplete or unclear, contain some inaccuracies, lack evidence of research and experimentation |\n",
    "| 6 – 7\t| Responses are informative and reflect insights, contain limited inaccuracies, contain some evidence of research and experimentation |\n",
    "| 8 – 10| Responses are accurate, informative, insightful, and contain clear evidence of research and experimentation|\n",
    "|   +   | |\n",
    "| 0\t    | Listed no contributions to discussion forum |\n",
    "| 1\t    | Listed one or more well-informed questions posed in discussion forum |\n",
    "| 2\t    | Listed one or more contributions to discussion forum, in response to questions |\n",
    "\n",
    "As indicated, marks are given based on the quantity and quality of constructive interaction in class and in the online forum associated with the course: questions *and* answers. We learn from one another, and this should be acknowledged. To ensure the recording of forum activity, please refer to your posts in the submission, and posts that assisted you.\n",
    "\n",
    "Formative feedback on submissions should be actively sought in the timetabled practical and tutorial sessions from course staff. Awarded marks will be published on Blackboard Grade Centre.\n",
    "\n",
    "### Workflow and submission\n",
    "\n",
    "You will submit your responses to [Coder Quiz](https://coderquiz.scmb.uq.edu.au/). You may submit as many times as you would like, and your last response before the due date will be graded. Please ensure that everything has submitted correctly to Coder Quiz by clicking on the 'View Submissions' link and verifying that all of your answers and code display correctly.\n",
    "\n",
    "Coder Quiz has the ability to check correctness of provided answers to **some** <span style=\"color:red\">*fixed-format questions*</span> prior to submission. This means that as you progress through the practical, you can check whether you are on the right track or not by clicking *Check Answers* button. Auto-marked questions are found in Part A only, and are indicated on the submission form.\n",
    "\n",
    "*Some* questions may ask you to provide the code you used to reach your answer. In this case, you need to save the relevant section of code. Once you're ready to submit, condense all required code into a **single .py file**, and upload this file to Coder Quiz. **Please ensure that code is *only provided for questions where it is explicitly requested*, and that all questions are labelled appropriately. If the marker cannot determine which code corresponds to each question, marks will not be awarded for those problems.** Submitted code is for visual inspection of your attempt, so it does not need to run (i.e. you do not need to include import statements) but it should be appropriately commented and understandable by a tutor. Marking criteria will consider whether you demonstrate an understanding of the underlying concepts. A separate Coder Quiz form is provided for code submission.\n",
    "\n",
    "Coder Quiz **does not** save or retrieve partial attempts, so we recommend storing your work and answers in a separate file; we strongly recommend you use this Jupyter notebook with additional markdown cells to save your ongoing work, then use Coder Quiz to validate and submit once you are complete.\n",
    "\n",
    "Remember to use the discussion board if you are unsure about how to approach a question or you are not able to get the correct result. That said, your submission must be the result of *your* understanding; if your answer contains anything that you are unable to explain or reproduce without the help of somebody else, you *must* acknowledge this. There is a separate prompt at the end of your submission to list posts on the discussion board that you have made, and posts that you have benefitted from. (In Ed Discussion, \"... / Copy Link\" for each such post.)\n",
    "\n",
    "\n",
    "### Resources\n",
    "\n",
    "* Course materials are available via Blackboard; pay attention to the weekly Python notebooks, in particular\n",
    "* Quick link to [How-to install binfpy](#howto_install_binfpy)\n",
    "* The UQ Bioinformatics Python Guide (on Blackboard)\n",
    "* The [Python 3 documentation]. For those unfamiliar with Python the [official tutorial] is recommended\n",
    "* The Software Carpentry [novice Python lessons]\n",
    "* [IPython's own notebook tutorial](http://nbviewer.jupyter.org/github/ipython/ipython/blob/3.x/examples/Notebook/Index.ipynb)\n",
    "* [Markdown cheatsheet] (Markdown is the syntax you use to write formatted text into cells in a notebook.)\n",
    "\n",
    "[Practical 1 ECP]: https://course-profiles.uq.edu.au/student_section_loader/section_5/108015#407455 \n",
    "[Python 3 documentation]: https://docs.python.org/3/\n",
    "[official tutorial]: https://docs.python.org/3/tutorial/index.html\n",
    "[novice python lessons]: http://swcarpentry.github.io/python-novice-inflammation/\n",
    "[Markdown cheatsheet]: https://github.com/adam-p/markdown-here/wiki/Markdown-Here-Cheatsheet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part A: motif discovery in proteins\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "We will be using the custom `binfpy` Python modules for this assignment. Instructions for accessing these modules are found in [How-to install binfpy](#howto_install_binfpy). Data files required for this practical are found in the same folder as this notebook.\n",
    "\n",
    "Import the following modules into Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sequence\n",
    "import prob\n",
    "import sym"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The classes you will need to complete this assignment are introduced in the weekly Python notebooks. Some simple examples of their use are provied below, however you can refer back to these notebooks for more comprehensive demonstrations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequence data\n",
    "\n",
    "We will be working with biological sequences (`sequence.Sequence`), which are defined by a series of characters from an alphabet (which in turn specifies all valid characters; `sym.Alphabet`). To understand sequences, it is informative to analyse their composition, and do so *probabilistically*. The Python module `prob.py` in `binfpy` has a number of useful classes, e.g. `prob.Distrib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "< A=0.33 C=0.00 G=0.00 T=0.67  >\n"
     ]
    }
   ],
   "source": [
    "myDNA = sym.Alphabet('ACGT') # Define an alphabet, but many are pre-defined in sym.py, e.g. sym.Protein_Alphabet\n",
    "d = prob.Distrib(myDNA)      # Create a probability distribution for our alphabet\n",
    "d.observe('A')               # Count a single observation of a character\n",
    "d.observe('T', 2)            # Count an observation of a character twice\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background distributions\n",
    "\n",
    "Here, you will construct a “background distribution” of amino acids that is suited for scoring motifs in human protein sequences. First, consider the background used for constructing BLOSUM62 (a popular substitution matrix) by viewing `blosum62.distrib`. This file can be read using the method `readDistrib` in `prob.py`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.05688622754491018"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bg = prob.readDistrib('blosum62.distrib')\n",
    "bg['S']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file `up_bac.fa` contains a random sample of bacterial sequences from [Uniprot](https://www.uniprot.org/) for which protein-level evidence is available. Construct your own background distribution using the sequences in `up_bac.fa`.\n",
    "\n",
    "Read through the `Distrib` class (in `prob.py`) to see how to construct a `Distrib` object and how to use its methods. \n",
    "\n",
    "**<span style=\"color:red\"> Problem A1: Report the probability of Serine in both the BLOSUM62 and bacterial background distributions. Enter to three decimal places. </span>**\n",
    "\n",
    "Enter your code in the cell below. A few lines are provided to get you started.\n",
    "\n",
    "#### Tips:\n",
    "* Consider what a background distribution represents and how it could be generated\n",
    "\n",
    "* An instance of the class `Distrib` can be inspected like so:\n",
    "\n",
    "```python\n",
    "   print(bg) # refers to __str__ defined for Distrib in prob.py\n",
    "```\n",
    "\n",
    "* Once you've generated your bacterial background distribution, save it in a file named `bac_bg.distrib` using the `writeDistrib` method:\n",
    "\n",
    "```python\n",
    "   bac_distrib.writeDistrib('bac_bg.distrib')\n",
    "```\n",
    "\n",
    "* A function `prob.writeDistribs` is also available, and will be useful later.\n",
    "\n",
    "* Note: the BLOSUM62 background is provided to you in the file `blosum62.distrib`. \n",
    "\n",
    "* If you are not already familiar with the `binfpy` library of Python code, please browse `sequence.py`, in particular the code to read FASTA files, and then access sequences' contents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "406 sequences loaded\n",
      "< A=0.05 C=0.05 D=0.05 E=0.05 F=0.05 G=0.05 H=0.05 I=0.05 K=0.05 L=0.05 M=0.05 N=0.05 P=0.05 Q=0.05 R=0.05 S=0.05 T=0.05 V=0.05 W=0.05 Y=0.05  >\n"
     ]
    }
   ],
   "source": [
    "# Write code to solve Problem A1 here\n",
    "bac_seqs = sequence.readFastaFile('up_bac.fa', sequence.Protein_Alphabet)\n",
    "print(len(bac_seqs), 'sequences loaded')\n",
    "# The following lines will create and print a clean-slate distribution\n",
    "bac_distrib = prob.Distrib(sequence.Protein_Alphabet)\n",
    "print(bac_distrib) # print out the distribution BEFORE any data have been looked-at\n",
    "# Your code to construct bacterial background below\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "< A=0.09 C=0.01 D=0.06 E=0.06 F=0.04 G=0.08 H=0.02 I=0.06 K=0.06 L=0.10 M=0.03 N=0.04 P=0.04 Q=0.04 R=0.05 S=0.06 T=0.06 V=0.07 W=0.01 Y=0.03  >\n"
     ]
    }
   ],
   "source": [
    "for seq in bac_seqs:\n",
    "    for pos in seq:\n",
    "        bac_distrib.observe(pos)\n",
    "bac_distrib['S']\n",
    "print(bac_distrib)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gibbs sampling\n",
    "\n",
    "Two methods for motif discovery are discussed in week 2 content: Gibbs sampling and expectation maximisation. This assignment makes use of the former. A key distinction of Gibbs sampling is its *stochastic* nature, as opposed to expectation maximisation which is *deterministic*.\n",
    "\n",
    "An overly-simplistic example of Gibbs sampling is provided below. This is to ensure that you undertand the concepts before we introduce some real data. For further clarification, refer to the lecture videos and relevant Python notebook.\n",
    "\n",
    "Here, we imagine a sequence alphabet which corresponds to the standard English alphabet.\n",
    "\n",
    "```python\n",
    "my_alpha = sym.Alphabet('ABCDEFGHIJKLMNOPQRSTUVWXYZ')\n",
    "```\n",
    "\n",
    "The file `example_seqs.fa` contains 150 sequences of length 80. We are told that a *motif* of length 7 is expected to appear in each of these sequences. In this case, the motif is a word relevant to this assignment's context. Given the relatively small length and number of sequences, you could probably identify this word by eye with relative ease. However, lets see if we can instead recover it using Gibbs sampling.\n",
    "\n",
    "For a more in-depth discussion of this approach in the context of biological sequence analysis, take a look at the publicaiton which first described it: [Lawrence et al. 1993](http://dx.doi.org/10.1126/science.8211139)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "import gibbs\n",
    "\n",
    "seqs = sequence.readFastaFile('example_seqsa.fa', my_alpha)\n",
    "W = ? # the width of the motif sought\n",
    "g = gibbs.GibbsMotif(seqs, W)\n",
    "q = g.discover(niter = 1000)\n",
    "```\n",
    "\n",
    "The above code can be found in a cell below. Extra print statements have been included to display the consensus sequence and the background distribution.\n",
    "\n",
    "**<span style=\"color:blue\">Problem A2: In the cell below, add comments *in your own words* explaining what the variables `seqs`, `g`, `q`, `p` and `a` represent (as an example, a comment has already been added to the variable `W`). Export this cell or copy and paste the commented code into a separate Python file and upload to Coder Quiz.</span>**\n",
    "\n",
    "**Note: For an example on how to easily save a Jupyter notebook cell to a Python file, [click here](#export).**\n",
    "\n",
    "Run Gibbs sampling several times using different window sizes (`W`), bearing in mind that we know the motif is 7 characters long. The size and number of iterations will influence the ability of the algorithm to discover anything significant. Compare the outcomes of different runs. You should do this by visually comparing motifs using [WebLogo], which takes an alignment like that printed by the example code. Determine approximately, on average, how many iterations are required for convergence. (Note: 'convergence' in the case of a stochastic algorithm does not mean that the log-likelihood score will not change between iterations. Rather, the score will gradually rise and then appear to 'level off' and fluctuate around a (possibly local) optima). Additionally, examine the log-likelihood of the final model and alignment. Can you distinguish between runs where Gibbs sampling has found the global optima, and those where it gets stuck in a local optima?\n",
    "\n",
    "The main motif discovery method is in `gibbs.discover`. By adjusting positions at which sequences are aligned, it essentially tries to maximise the (log) ratio of the foreground over background probabilities for observed sequences. It prints a sum of these log ratios. (The sum of log $x_1, x_2, ..., x_n$, is equal to the product of $x_1, x_2, ..., x_n$, where $x_1, x_2, ..., x_n$ are the ratios of foreground to background probabilities.)\n",
    "\n",
    "Naively your comparison could be based on the sums above but if you think about it, and as discussed in the paper by [Lawrence et al.](http://dx.doi.org/10.1126/science.8211139), they will inevitably be greater for longer motifs (at least as long as there is a modicum of conservation in included positions). So this score should not be relied on in isolation.\n",
    "\n",
    "[WebLogo]: http://weblogo.threeplusone.com\n",
    "\n",
    "**<span style=\"color:red\">Problem A3: Identify the 7-residue mystery sub-motif in these example sequences. </span>**\n",
    "\n",
    "**<span style=\"color:green\">Problem A4: Examine the background disribution obtained from Gibbs sampling, and suggest how these example sequences were likely generated. Is it likely that a similar phenomenon would be observed in real biological sequences?</span>**\n",
    "\n",
    "Due to the stochastic nature of this approach, the path the algorithm takes through the sequence-probability space varies between runs. In this case, a perfect copy of the 'motif' is present in each of the input sequences, and hence Gibbs sampling converges on this intended sub-sequence in a decent proportion of runs. In subsequent sections of this assignment, consider whether such strict conservation is biologically realistic, and how this will impact the algorithm's behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LL @     0=\t87.34\n",
      "LL @   100=\t759.74\n",
      "LL @   200=\t1453.37\n",
      "LL @   300=\t1785.69\n",
      "LL @   400=\t2134.16\n",
      "LL @   500=\t2318.80\n",
      "LL @   600=\t2318.96\n",
      "LL @   700=\t2314.48\n",
      "LL @   800=\t2321.83\n",
      "LL @   900=\t2320.87\n",
      "LL @  1000=\t2345.58\n",
      "LL @  1100=\t2343.25\n",
      "LL @  1200=\t2344.55\n",
      "LL @  1300=\t2379.35\n",
      "LL @  1400=\t2370.42\n",
      "LL @  1500=\t2379.85\n",
      "LL @  1600=\t2380.27\n",
      "LL @  1700=\t2380.62\n",
      "LL @  1800=\t2380.25\n",
      "LL @  1900=\t2378.25\n",
      "Consensus: CHPROTE\n",
      "Background distribution: < A=0.03 B=0.04 C=0.04 D=0.04 E=0.04 F=0.04 G=0.04 H=0.04 I=0.05 J=0.04 K=0.03 L=0.03 M=0.04 N=0.05 O=0.04 P=0.04 Q=0.04 R=0.04 S=0.04 T=0.04 U=0.04 V=0.04 W=0.04 X=0.04 Y=0.04 Z=0.04  >\n",
      "\n",
      "Identified sub-sequences of length:  7 in input sequences.\n",
      "Seq_1\tDFPROTE\n",
      "Seq_2\tAGPROTE\n",
      "Seq_3\tLIPROTE\n",
      "Seq_4\tAUPROTE\n",
      "Seq_5\tRWPROTE\n",
      "Seq_6\tYKPROTE\n",
      "Seq_7\tCGPROTE\n",
      "Seq_8\tTJPROTE\n",
      "Seq_9\tPNPROTE\n",
      "Seq_10\tKVPROTE\n",
      "Seq_11\tNQPROTE\n",
      "Seq_12\tEHPROTE\n",
      "Seq_13\tDJPROTE\n",
      "Seq_14\tHXPROTE\n",
      "Seq_15\tAAPROTE\n",
      "Seq_16\tOHPROTE\n",
      "Seq_17\tZEPROTE\n",
      "Seq_18\tMMPROTE\n",
      "Seq_19\tJXPROTE\n",
      "Seq_20\tYDPROTE\n",
      "Seq_21\tCKPROTE\n",
      "Seq_22\tTYPROTE\n",
      "Seq_23\tZKPROTE\n",
      "Seq_24\tUZPROTE\n",
      "Seq_25\tCAPROTE\n",
      "Seq_26\tHZPROTE\n",
      "Seq_27\tIZPROTE\n",
      "Seq_28\tNJPROTE\n",
      "Seq_29\tNXPROTE\n",
      "Seq_30\tSGPROTE\n",
      "Seq_31\tCMPROTE\n",
      "Seq_32\tCPPROTE\n",
      "Seq_33\tCHPROTE\n",
      "Seq_34\tDWPROTE\n",
      "Seq_35\tLTPROTE\n",
      "Seq_36\tNWPROTE\n",
      "Seq_37\tFVPROTE\n",
      "Seq_38\tKQPROTE\n",
      "Seq_39\tTHPROTE\n",
      "Seq_40\tGFPROTE\n",
      "Seq_41\tSQPROTE\n",
      "Seq_42\tQNPROTE\n",
      "Seq_43\tOKPROTE\n",
      "Seq_44\tCZPROTE\n",
      "Seq_45\tDZPROTE\n",
      "Seq_46\tMHPROTE\n",
      "Seq_47\tEWPROTE\n",
      "Seq_48\tUGPROTE\n",
      "Seq_49\tJJPROTE\n",
      "Seq_50\tYUPROTE\n",
      "Seq_51\tRNPROTE\n",
      "Seq_52\tMWPROJE\n",
      "Seq_53\tOVPROTE\n",
      "Seq_54\tMUPROTE\n",
      "Seq_55\tKHPROTE\n",
      "Seq_56\tRVPROTE\n",
      "Seq_57\tKKPROTE\n",
      "Seq_58\tQUPROTE\n",
      "Seq_59\tWEPROTE\n",
      "Seq_60\tNQPROTE\n",
      "Seq_61\tONPROTE\n",
      "Seq_62\tPGPROTE\n",
      "Seq_63\tUNPROTE\n",
      "Seq_64\tJUPROTE\n",
      "Seq_65\tXDPROTE\n",
      "Seq_66\tRLPROTE\n",
      "Seq_67\tDAPROTE\n",
      "Seq_68\tVQPROTE\n",
      "Seq_69\tLVCRESQ\n",
      "Seq_70\tAJPROTE\n",
      "Seq_71\tNRPROTE\n",
      "Seq_72\tAZPROTE\n",
      "Seq_73\tKQPROTE\n",
      "Seq_74\tLVPROTE\n",
      "Seq_75\tOFPROTE\n",
      "Seq_76\tAZPROTE\n",
      "Seq_77\tBXPROTE\n",
      "Seq_78\tVXPROTE\n",
      "Seq_79\tUDPROTE\n",
      "Seq_80\tRNPROTE\n",
      "Seq_81\tJSPROTE\n",
      "Seq_82\tNQPROTE\n",
      "Seq_83\tHXPROTE\n",
      "Seq_84\tJBPROTE\n",
      "Seq_85\tABPROTE\n",
      "Seq_86\tUBPROTE\n",
      "Seq_87\tIAPROTE\n",
      "Seq_88\tDUPROTE\n",
      "Seq_89\tUGPROTE\n",
      "Seq_90\tUTPROTE\n",
      "Seq_91\tOLPROTE\n",
      "Seq_92\tCOPROTE\n",
      "Seq_93\tWAPROTE\n",
      "Seq_94\tBVPROTE\n",
      "Seq_95\tNZPROTE\n",
      "Seq_96\tOFPROTE\n",
      "Seq_97\tHLPROTE\n",
      "Seq_98\tYOPROTE\n",
      "Seq_99\tZSAROTS\n",
      "Seq_100\tECPROTE\n",
      "Seq_101\tTKPROTE\n",
      "Seq_102\tMZPROTE\n",
      "Seq_103\tIMPROTE\n",
      "Seq_104\tLWPROTE\n",
      "Seq_105\tJWPROTE\n",
      "Seq_106\tXNBRUNT\n",
      "Seq_107\tVUPROTE\n",
      "Seq_108\tRPPROTE\n",
      "Seq_109\tROPROTE\n",
      "Seq_110\tJIPROTE\n",
      "Seq_111\tLWPROTE\n",
      "Seq_112\tUEPROTE\n",
      "Seq_113\tFJPROTE\n",
      "Seq_114\tZFPROTE\n",
      "Seq_115\tWPPROTE\n",
      "Seq_116\tFEPROTE\n",
      "Seq_117\tYXLRYBE\n",
      "Seq_118\tILPROTE\n",
      "Seq_119\tWBPROTE\n",
      "Seq_120\tPDPROTE\n",
      "Seq_121\tOQPROTE\n",
      "Seq_122\tPUPROTE\n",
      "Seq_123\tVRPROTE\n",
      "Seq_124\tIKPROTE\n",
      "Seq_125\tJHPROTE\n",
      "Seq_126\tNYPROTE\n",
      "Seq_127\tLIPROTE\n",
      "Seq_128\tVQPROTE\n",
      "Seq_129\tXKPROTE\n",
      "Seq_130\tIGPROTE\n",
      "Seq_131\tMAPROTE\n",
      "Seq_132\tFLPROTE\n",
      "Seq_133\tOIPROTE\n",
      "Seq_134\tSXPROTE\n",
      "Seq_135\tBEPROTE\n",
      "Seq_136\tCUPROTE\n",
      "Seq_137\tISPROTE\n",
      "Seq_138\tWDPROTE\n",
      "Seq_139\tPTPROTE\n",
      "Seq_140\tFTPROTE\n",
      "Seq_141\tCVPROTE\n",
      "Seq_142\tWMPROTE\n",
      "Seq_143\tYMPROTE\n",
      "Seq_144\tXPPROTE\n",
      "Seq_145\tCBPROTE\n",
      "Seq_146\tYOPROTE\n",
      "Seq_147\tTHPROTE\n",
      "Seq_148\tVKPROTE\n",
      "Seq_149\tRHPROTE\n",
      "Seq_150\tDPPROTE\n"
     ]
    }
   ],
   "source": [
    "import gibbs\n",
    "import sym\n",
    "import sequence\n",
    "\n",
    "my_alpha = sym.Alphabet('ABCDEFGHIJKLMNOPQRSTUVWXYZ')\n",
    "\n",
    "seqs = sequence.readFastaFile('example_seqs.fa', my_alpha)\n",
    "W = 7 # the width of the motif sought\n",
    "g = gibbs.GibbsMotif(seqs, W)\n",
    "q = g.discover(niter = 2000)\n",
    "print('Consensus: ', end='')\n",
    "for pos in q:\n",
    "    print(pos.getmax(), end='')\n",
    "print()\n",
    "p = g.getBackground()\n",
    "print('Background distribution:', p)\n",
    "print()\n",
    "a = gibbs.getAlignment(seqs, q, p)\n",
    "k = 0\n",
    "print('Identified sub-sequences of length: ', W, 'in input sequences.')\n",
    "for seq in seqs:\n",
    "    #print(\"%s \\t%s\" % (seq.name, seq[a[k]:a[k]+W]))\n",
    "    print(f\"{seq.name}\\t{seq[a[k]:a[k]+W]}\")\n",
    "    k += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before moving on, run Gibbs sampling with W=7 until you find the desired motif. Save the foreground and background probability distributions in case you need to re-load them later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save your foreground (q) and background (p) distributions\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Motif searching\n",
    "\n",
    "Above, we identified our mystery 'motif'. Now, you are going to construct position specific scoring matrices (log-odds matrices), otherwise called “position weight matrices” (PWMs) which represent this motif and can be used to search for occurrences of the motif. We need to be able to search for a motif because it may not always be in the same position of a sequence.\n",
    "\n",
    "One of the variables that you commented on above contains the probabilities defining your motif that you need to construct PWMs. We also previously defined two backgrounds, as well as the background generated in Gibbs sampling (based on the\" training sequences excluding the motifs). Both the probability defining the motif and a background are required to generate a PWM. The class `PWM` is defined in the `sequence.py` module. An example of constructing a PWM is included below, but their use is also demonstrated in weekly notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A\t +0.31  +0.16  -1.60  -4.90  -4.90  -4.90  -4.90\n",
      "B\t -0.60  -0.10  -1.67  -4.97  -4.97  -1.67  -4.97\n",
      "C\t +0.63  -1.73  -1.73  -5.03  -5.03  -5.03  -5.03\n",
      "D\t +0.25  -0.09  -4.96  -4.96  -4.96  -4.96  -4.96\n",
      "E\t -0.59  -0.08  -4.96  -4.96  -1.66  -4.96  +3.28\n",
      "F\t -0.19  -0.19  -5.06  -5.06  -5.06  -5.06  -5.06\n",
      "G\t -1.71  +0.21  -5.00  -5.00  -5.00  -5.00  -5.00\n",
      "H\t -0.27  +0.53  -4.92  -4.92  -4.92  -4.92  -4.92\n",
      "I\t -0.29  -0.69  -5.35  -5.35  -5.35  -5.35  -5.35\n",
      "J\t +0.40  +0.12  -4.94  -4.94  -4.94  -1.64  -4.94\n",
      "K\t -0.02  +0.44  -4.90  -4.90  -4.90  -4.90  -4.90\n",
      "L\t +0.34  +0.00  -1.58  -4.87  -4.87  -4.87  -4.87\n",
      "M\t +0.07  -0.11  -4.99  -4.99  -4.99  -4.99  -4.99\n",
      "N\t +0.16  -0.09  -5.30  -5.30  -5.30  -2.00  -5.30\n",
      "O\t +0.51  -0.29  -4.95  -4.95  +3.30  -4.95  -4.95\n",
      "P\t -0.18  -0.18  +3.18  -5.06  -5.06  -5.06  -5.06\n",
      "Q\t -1.00  +0.49  -4.97  -4.97  -4.97  -4.97  -1.67\n",
      "R\t +0.22  -1.15  -5.12  +3.14  -5.12  -5.12  -5.12\n",
      "S\t -0.67  -0.67  -5.04  -5.04  -5.04  -1.74  -1.74\n",
      "T\t -0.13  -0.35  -5.01  -5.01  -5.01  +3.23  -1.71\n",
      "U\t +0.39  +0.51  -4.95  -4.95  -1.66  -4.95  -4.95\n",
      "V\t +0.08  +0.37  -4.97  -4.97  -4.97  -4.97  -4.97\n",
      "W\t +0.10  +0.38  -4.96  -4.96  -4.96  -4.96  -4.96\n",
      "X\t -0.29  +0.40  -4.95  -4.95  -4.95  -4.95  -4.95\n",
      "Y\t +0.26  -0.98  -4.95  -4.95  -1.65  -4.95  -4.95\n",
      "Z\t -0.39  +0.42  -5.04  -5.04  -5.04  -5.04  -5.04\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pwm = sequence.PWM(q, p)\n",
    "print(pwm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PWM we've just created uses the background from the proteins in `example_seqs.fa`. In the following exercise, assume that the new sequences' backgrounds are similar to that of `example_seqs.fa`. Why this is this an important consideration?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a couple of ways of using a PWM to score sequences as exemplified below. Here, we define four sequences in which we'll search for our motif."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_alpha = sym.Alphabet('ABCDEFGHIJKLMNOPQRSTUVWXYZ')\n",
    "seq1 = sequence.Sequence('XZGUIZMVHVEVZRHJJUCXKBOCBXEUMIXQZITGTAMLPTJSNWRFDDOBANZLGQQZHONFVHKIKJYXCWCUIYOV', my_alpha)\n",
    "seq2 = sequence.Sequence('SMATNQFJYWFHPRUTYINWLNAAOZXWUQKWCVOZSIPPLCONJRYVKBOBTIAIUVOBUNUANFVSEQGMHMIORBRC', my_alpha)\n",
    "seq3 = sequence.Sequence('OQUPKMLEZXELCMIVMVUDBPEAAPFJPROTEINFNBYNLDXTNYLCLUCHROCOOHFWQYCIHVVKLAEASAZXUINI', my_alpha)\n",
    "seq4 = sequence.Sequence('JUYROTEMNHKREGNRPVSHLKCUUSWIMWRLRATIPROTEINFMEHMFMXIGFUVJEEWQLGHUOLHLVALELRGDYMS', my_alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`seq1` does not contain a copy of the motif. Using the PWM class method `PWM.search` won't return any hits, as the default threshold score is 0. Given the highly specific nature of our motif, we would expect random sub-sequences of length W to score poorly (<0). This can be understood by viewing the PWM (printed above). Characters which are not the consensus character in a given position have very negative log-odds, providing little flexibility to variation in the motif. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 'JUYROTE', 8.909350821056016), (34, 'TIPROTE', 15.304611078951343)]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Search for the motif in the above sequences\n",
    "pwm.search(seq4,0)  # A search in seq1 returns no hits for our motif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Try searching the other three sequences using your PWM, and understand what `PWM.search` outputs. You should also see examples of hits where a slightly 'mutated' version of the motif is present.\n",
    "\n",
    "**<span style=\"color:red\">Problem A5: Which sequence defined above appears to contain two copies of the mystery motif? </span>**\n",
    "\n",
    "A second method, `PWM.maxscore`, will return the position and score of the highest scoring window of length W, regardless of the score. See below for an example. Observe the output for each of the sequences above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum score -13.170229717869052 occurs at position 46\n",
      "All matches above 0 are:\n"
     ]
    }
   ],
   "source": [
    "myseq = seq1\n",
    "result = pwm.maxscore(myseq)\n",
    "print('The maximum score', result[0], 'occurs at position', result[1])\n",
    "threshold = 0\n",
    "result = pwm.search(myseq, threshold)\n",
    "print('All matches above', threshold, 'are:')\n",
    "for r in result:\n",
    "    print('\\t', r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The metallo-beta-lactamase superfamily\n",
    "\n",
    "The following information is relevant to the sequence data you will use throughout this assignment.\n",
    "\n",
    "The term 'homologous protein' is often associated with different versions of the same gene across species - distinct in sequence, yet ultimately undertaking a highly similar function. This over-simplification conflates 'homologs' and 'orthologs', and fails to consider modes of molecular evolution beyond just speciation, such as gene duplication and subsequent functional specialistaion. \n",
    "\n",
    "Studying the composition of large and diverse protein families - often termed 'superfamilies' - is further complicated by the phenomenon of molecular 'promiscuity', whereby some proteins are seemingly capable of catalysing multiple distinct reaction types, often on different substrates. Tracing the evolution of function within such families is a complex undertaking, particulary where the level of sequence divergence is high and the last common ancestor is billions of years old.\n",
    "\n",
    "The metallo-beta-lactamase (MBL) superfamily is an ancient group of related enzymes with an extraordinarily broad catalytic repertoire. Named for the first such activity to be discovered (degradation of beta-lactam antibiotics), these enzymes possess a highly conserved $\\alpha\\beta\\beta\\alpha$ fold.\n",
    "\n",
    "**Figure 1: Conserved fold of the MBL superfamily**\n",
    "\n",
    "<img src=\"mbl_fold.png\" width=\"400\" height=\"400\" />\n",
    "\n",
    "Despite this structural conservation, pairwise sequence similarities can be as low as 5%, meaning that alignment, let alone phylogenetic inference, is anything but straightforward. How then, can we identify MBL superfamily members without undertaking the laborious process of structural determination?\n",
    "\n",
    "Incredibly, the sequence and structural positions implicated in catalysis are the same throughout the entire superfamily. These positions localise to two distinct sites in 3D space upon folding, each of which can coordinate metal ions (hence *metallo*-beta-lactamase) that participate in catalysis.\n",
    "\n",
    "**Figure 2: Metal binding residues of the MBL superfamily**\n",
    "\n",
    "<img src=\"mbl_metals.png\" width=\"250\" height=\"250\" />\n",
    "\n",
    "Residues present at these sites are well conserved, and invariably include Histidine. Subtle differences occur at some positions which has been linked to specific metal preference and/or catalytic function.\n",
    "\n",
    "All protein entries in Uniprot are scanned for the presence of MBL domains, among many others, using a complementary database, [Interpro](https://www.ebi.ac.uk/interpro/). Take a look at entry [P52700](uniprot.org/uniprotkb/P52700/entry), and visualise the annotation of various protein domains by Interpro.\n",
    "\n",
    "**<span style=\"color:red\">Problem A6: Provide the start and end positions of the Interpro-annotated MBL domain for P52700.</span>**\n",
    "\n",
    "Interpro also attempts to classify proteins into homologous superfamilies.\n",
    "\n",
    "**<span style=\"color:red\">Problem A7: What alternative name does Interpro use for the MBL superfamily?</span>** \n",
    "\n",
    "Several different functions are attributed to proteins in the MBL superfamily. Those for which beta-lactamase is the primary activity are often referred to as *True MBLs* or *Class B MBLs*. These can be further divided into types B1, B2 and B3. The marked ability of *promiscuous* MBL superfamily members to catalyse multiple reaction types is of growing interest to researchers.\n",
    "\n",
    "Evolutionarily distinct from these enzymes are Classes A, C and D, collectively termed *Serine beta-lactamases* (SBLs), which form part of another protein superfamily. While both capable of degrading beta-lactam antibiotics, MBLs and SBLs gained this ability via *convergent evolution.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Datasets\n",
    "\n",
    "You are provided with a number of datasets for use in the following exercises. These are briefly described here:\n",
    "\n",
    "* `mbl_seqs.fa` - All bacterial sequences with an annotated MBL domain were extracted from Uniprot and clustered at 40% sequence identity. A representative from each cluster was placed in this file.\n",
    "* `positives.fa` - A set of 20 proteins verified as Class B beta-lactamases.\n",
    "* `negatives.fa` - A set of 20 random proteins from various model species known **not** to contain an MBL domain.\n",
    "* `active.fa` - A set of proteins found to degrade beta-lactam antibiotics in a high-throughput screening experiment (further context provided below)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Motif discovery in MBLs\n",
    "\n",
    "The above examples provide the essential tools for linear motif discovery with Gibbs sampling, and can be generalised to any biological sequence type. You will now perform motif discovery and searching on MBL sequences. Note that there are no single correct solutions to the following problems, however your processes must be sufficiently justified.\n",
    "\n",
    "#### Perform Gibbs sampling on the full set of bacterial MBL sequences\n",
    "\n",
    "`mbl_seqs.fa` contains a representative sample of all bacterial MBL superfamily members in Uniprot. Using your knowledge of Gibbs sampling and MBL proteins, you will now carry out motif discovery on this set of sequences. While the seven highly conserved metal-binding residues do not all occur together in the *linear* sequence, you will notice in Figure 2 that four of these residues are in close proximity. Assume that the location of these positions represents the most highly conserved sequence region.\n",
    "\n",
    "You should consider several factors while implementing your motif discovery method.\n",
    "* What window size will you use?\n",
    "* How many iterations are required?\n",
    "* How will you determine the 'best' motif if different results are obtained between runs?\n",
    "* Given your knowledge of conserved residues in MBLs, what may you expect to see in your consensus motif?\n",
    "\n",
    "You will want to run Gibbs multiple times when considering the above parameters. Once happy with your motif, ensure that you save any relevant distributions for subsequent searching."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code for motif discovery in MBLs here...\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:green\">Problem A8: State and justify the window size used in your motif discovery process.</span>**\n",
    "\n",
    "**<span style=\"color:green\">Problem A9: State and justify the number of iterations used.</span>**\n",
    "\n",
    "**<span style=\"color:red\">Problem A10: Provide the consensus sequence of your final motif.</span>**\n",
    "\n",
    "**Please ensure Coder Quiz accepts your motif before moving on to subsequent sections.**\n",
    "\n",
    "**<span style=\"color:blue\">Problem A11: Provide the code used for motif discovery from the above cell.</span>** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Construct a PWM for motif searching "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code for constructing your PWM here...\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:green\">Problem A12: Describe and justify the background distribution you used.</span>**\n",
    "\n",
    "`print` the PWM object to view the log-odds matrix.\n",
    "\n",
    "**<span style=\"color:green\">Problem A13: Which positions appear to have the highest levels of conservation? What is the significance of these residues?</span>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Determine a score threshold for MBL domain-containing proteins\n",
    "\n",
    "In order to search for MBL domains among unknown sequences, an appropriate score threshold is required. Use the two sequence sets, `positives.fa` and `negatives.fa` (described above) to select this threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code for treshold determination here...\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:green\">Problem A14: State and justify the score threshold you have chosen, making reference to relevant observations from the positive and negative sequences.</span>**\n",
    "\n",
    "**<span style=\"color:green\">Problem A15: Given your observations and chosen threshold, comment generally on the likely sensitivity and specificity of your PWM for identifying MBL superfamily proteins. Calculations of these values are *not* required. </span>**\n",
    "\n",
    "**<span style=\"color:blue\">Problem A16: Provide the code used for threshold determination from the above cells.</span>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Search a set of unknown proteins for your motif\n",
    "\n",
    "A biochemist purified several proteins for which beta-lactamse activity was suspected, and high-throughput screening was performed to identify those capable of degrading at least one beta-lactam substrate. The genes corresponding to these active proteins were sequenced, and subsequently identified by BLAST search of the Uniprot database. These records were downloaded in FASTA format and are available in the file `active.fa`.\n",
    "\n",
    "Check for the existence of your motif in each of these beta-lactamase sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code for searching the beta-lactam degrading enzymes here...\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:red\">Problem A17: How many of these proteins score above your chosen threshold? </span>**\n",
    "\n",
    "**<span style=\"color:red\">Problem A18: Provide the Uniprot accession for one of the highest scoring proteins. </span>**\n",
    "\n",
    "**<span style=\"color:blue\">Problem A19: Provide the code used for motif searching from the above cell(s).</span>**\n",
    "\n",
    "The Sequence class has an attribute `.info` which may be helpful in answering the following questions.\n",
    "\n",
    "**<span style=\"color:green\">Problem A20: Explain why several of the proteins capable of beta-lactamase activity do not appear to contain your motif. </span>**\n",
    "\n",
    "Identify at least one sequence among your hits with an unusual description, then take a look at this paper by [Fröhlick *et al.*](https://academic.oup.com/peds/article/doi/10.1093/protein/gzab013/6294778) \n",
    "\n",
    "**<span style=\"color:green\">Problem A21: Provide a likely explanation for the presence of your motif in this sequence, as well as the protein's apparent beta-lactamase activity.</span>**\n",
    "\n",
    "**<span style=\"color:green\">Problem A22: Suggest how the motif discovery processes undertaken in this section could be altered to more specifically identify *true* beta-lactam degrading enzymes from the MBL superfamily.</span>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part B: inferring ancestral protein sequences using graphical models\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "\n",
    "Part B problems are less guided and require some choices to be made by you, informed by your understanding of the relevant probabilistic and biological concepts. Activity on the discussion board around these problems is encouraged, however ensure that no code or solutions are shared publicly. Marks are available for discussion board participation (details in *Marking* section above).\n",
    "\n",
    "Import the following modules into Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: sympy in /opt/conda/lib/python3.10/site-packages (1.11.1)\n",
      "Requirement already satisfied: mpmath>=0.19 in /opt/conda/lib/python3.10/site-packages (from sympy) (1.2.1)\n"
     ]
    }
   ],
   "source": [
    "!pip install sympy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import asr\n",
    "import sequence\n",
    "import prob\n",
    "import sym"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`asr.py` makes use of the package, `sympy` for matrix operations. If you find that it is not available to import, it is easily installable via standard methods (see https://pypi.org/project/sympy/)\n",
    "\n",
    "Given a set of evolutionarily related present day ('extant') biological sequences, we are often interested in the paths by which said sequences diverged from a common ancestor. In the absence of ancient DNA from fossils (or a time machine), such queries can only be answered with predictions. Techniques for doing so are commonly termed ancestral sequence reconstruction (ASR). Various paradigms exist for performing ASR. One such paradigm, Maximum Likelihood (ML), will be discussed and used in this assignment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assumptions in ancestral sequence reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several key assumptions underly ML-ASR, and are discussed below.\n",
    "\n",
    "The first of these is the phylogenetic topology describing the relationships between extant sequences. Inferring these relationships is a separate (but related) matter which we won't concern ourselves with in this assignment. Some are discussed in SCIE2100/BINF6000, so refer back to those notes if you'd like a refresher; however, all you need to know here is that ML-ASR takes a fixed phylogenetic tree as input. Extant sequences are represented as terminal, or leaf, nodes, while ancestral sequences are represented as internal nodes. Given a bifurcating tree with *m* leaf nodes, the total number of nodes, *n*, is *2m-1*. The below example shows a tree with 4 leaf nodes (A, B, C, and D), and 3 internal nodes (X, Y, and Z).\n",
    "\n",
    "<img src=\"example_tree.png\" width=\"600\" height=\"600\" />\n",
    "\n",
    "In general, we have access to the sequences at extant nodes. A common motivation of ASR is to determine the most likely *joint* assignment of sequences to all internal nodes. The fixed nature of our input phylogenetic tree allows us to make some helpful assumptions of *conditional independence* between nodes. In fact, the structure of phylogenetic trees can be intuitively imagined as Bayesian networks in which we assume that sequence characters at nodes are dependent only on the characters at their direct parents, and the corresponding branch length. If the concept of Bayesian network representations for phylogenetic trees is unfamiliar, refer to week 4 lecture materials. \n",
    "\n",
    "**<span style=\"color:green\">Problem B1: Assuming a Bayesian network representation of the phylogenetic tree in the figure above, provide an expression for the joint probability of some arbitrary assignment of characters, $\\{\\alpha_A,\\alpha_B,\\alpha_C,\\alpha_D,\\alpha_X,\\alpha_Y,\\alpha_Z\\}$, in terms of all relevant conditional and/or prior probabilities. </span>**\n",
    "\n",
    "**<span style=\"color:green\">Problem B2: Briefly explain why the assumptions of conditional independence underlying Bayesian networks are sensbile for inferring ancestral sequence characters. </span>**\n",
    "\n",
    "Another major assumption often (although not always) made in ASR is that of *column independence*. It specifies that all aligned columns are considered separately when inferring ancestral states. Effectively, this means that for a given phylogenetic tree built from an alignment of length *L*, determining the most maximally likely set of ancestral sequences involves inference on *L* independent Bayesian networks.\n",
    "\n",
    "**<span style=\"color:green\">Problem B3: Briefly explain why the assumption of column independence for ASR may be considered a statistical simplification, but biologically unrealistic. </span>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time-dependent amino acid substitution models\n",
    "\n",
    "In week 4 lectures, we discussed the concept of sequence evolution along branches as continuous-time Markov chains. Specifically, we considered the Jukes-Cantor model for DNA, which makes an assumption of equal mutation rates to any character other than the current one. While this model requires only one rate parameter, the key assumption is not well supported in the reality of evolution.\n",
    "\n",
    "Time-dependent substitution models are generally represented by an *instantaneous rate matrix*, or IRM. This matrix contains parameters relating to the likelihood of transitioning between two characters over time. The DNA sequence alphabet contains only four characters, and hence only 12 character 'transitions' (**excluding** cases where the character does not change) are possible. Many biological substitution models are *time-reversible*, which means that the direction of evolution on a tree is unimportant to the model. This is useful when the *root* of a phylogenetic tree is not confidently known. For time-reversible models, only one rate parameter is required to determine both transition probabilities for a pair of characters. For a DNA alphabet, this reduces the number of required rate parameters to 6.\n",
    "\n",
    "**<span style=\"color:red\">Problem B4: Assuming a time-reversible substitution model, provide a mathematical expression for the number of rate parameters required given a sequence alphabet of size $N$. How many rate parameters are therefore required for a standard amino acid alphabet (excluding stop and gap characters)?</span>**\n",
    "\n",
    "Substitution models also require an equilibrium probability distribution over the alphabet of characters. This can be thought of as the *prior* probability distribution of observing a character at any node.\n",
    "\n",
    "For a given substitution model, probabilities of transitioning from one state to another over a specific time $t$ are determined through some basic linear algebra on the equilibrium and instantaneous rate matrices. You aren't required to understand the mathematics underlying these, however if you're interested, it is demonstrated in the `asr.py` code. This results in a matrix, $P(t)$, which directly states probabilities of a particular sequence character at a child node conditional on its immediate ancestor and preceding branch length.\n",
    "\n",
    "Luckily, several pre-defined substitution models for amino acid sequences are available. Currently, only the Jones-Taylor-Thornton (JTT) model - which we will use in this section - is implemented in `asr.py`. Let's take a look at this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "JTT = asr.MODELS['JTT']   # The pre-defined JTT model in asr.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will investigate transition probabilities from a parent node to a single child node for a given sequence position. Hypothetically, we know that the ancestral character at this position was alanine (A). Let's take a look at the *equilibrium* frequency of A in the JTT model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0768620000000000\n"
     ]
    }
   ],
   "source": [
    "print(JTT.priorProb('A'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's consider how likely it is (based on our model) that this A will transition to various amino acids over a branch length of $t=10$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A to V:\n",
      "0.0474905331230159\n",
      "\n",
      "A to L:\n",
      "0.0104986827604496\n",
      "\n",
      "A to K:\n",
      "0.00944150861945459\n"
     ]
    }
   ],
   "source": [
    "print('A to V:')\n",
    "print(JTT.condProb('A', 'V', 10))\n",
    "print()\n",
    "print('A to L:')\n",
    "print(JTT.condProb('A', 'L', 10))\n",
    "print()\n",
    "print('A to Y:')\n",
    "print(JTT.condProb('A', 'Y', 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:green\">Problem B5: Suggest a biological explanation for the differences in transition probabilities observed.</span>**\n",
    "\n",
    "In addition to the probability of transitioning between two characters, we're also often interested in the probability of a character being maintained. Investigate the probability of an alanine residue being maintained over branch lengths of $t=1$, $t=5$, $t=10$, $t=50$, $t=100$, $t=500$, and $t=1000$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t=1: 0.950237226392922\n",
      "t=5: 0.778061310669179\n",
      "t=10.5: 0.598272508900055\n",
      "t=100: 0.0888551064076305\n",
      "t=500: 0.0768625237487430\n",
      "t=1000: 0.0768619232195663\n"
     ]
    }
   ],
   "source": [
    "# Write your code here\n",
    "branch_lengths = [1,5,10.50,100,500,1000]\n",
    "\n",
    "for t in branch_lengths:\n",
    "    print(f\"t={t}: {JTT.condProb('A','A',t)}\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:green\">Problem B6: Describe what you observe with increasing time. What is happening to the transition probabilities in terms of the model?</span>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ancestral inference in phylogenetic trees\n",
    "\n",
    "We're generally faced with far more complex situations than a single parent-child relationship. When performing evolutionary inference on actual phylogenetic trees, we must consider many dependencies - and hence many conditional probabilities - simultaneously. \n",
    "\n",
    "To demonstrate some further functionality of `asr.py`, we'll first consider a simple phylogenetic tree. It is represented below as a Bayesian network.\n",
    "\n",
    "<img src=\"simple_bn.png\" width=\"400\" height=\"400\" />\n",
    "\n",
    "Network nodes can be represented as `PhyloBNode` objects which hold information such as the node's parent (if one exists), the preceding branch length, and the alphabet of allowable sequence characters. If the character state at the node is known, such is the case for the leaf nodes in the figure above, then these can also be annotated. From the code below, identify which variables represent each node in the example Bayesian network. Ensure that you understand each input to the PhyloBNode constructors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "JTT = asr.MODELS['JTT']\n",
    "\n",
    "ancestor_1 = asr.PhyloBNode(JTT, label='ancestor_1')\n",
    "ancestor_2 = asr.PhyloBNode(JTT, parent=ancestor_1, distance=1, label='ancestor_2')\n",
    "\n",
    "child_1 = asr.PhyloBNode(JTT, parent=ancestor_1, distance=2, label='child_1', annot='K')\n",
    "child_2 = asr.PhyloBNode(JTT, parent=ancestor_2, distance=1, label='child_2', annot='A')\n",
    "child_3 = asr.PhyloBNode(JTT, parent=ancestor_2, distance=1, label='child_3', annot='V')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our nodes constructed, we now need to connect them in a `PhyoBNet` object which is representative of the network shown above. We can do this by first initialising our network with the root node, and subsequently adding a list of all other nodes. It's important that the parent fields for all nodes are correct, and that no duplicate labels exist! Note that the instantiation of the network adds the root node to the network. You don't need to (and should not) add it again - this will cause unusual behaviour of the object's methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "bn = asr.PhyloBNet(root=ancestor_1)\n",
    "bn.addNodes([ancestor_2, child_1, child_2, child_3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Provided that node instantiations are consistent with a valid phylogenetic tree, the network's structure is automatically resolved using the information provided in the nodes. To demonstrate this, we can extract the children of a given node, even though we haven't explicitly provided this information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['child_2', 'child_3']\n"
     ]
    }
   ],
   "source": [
    "print([child.label for child in bn.getChildrenOf('ancestor_2')])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Joint maximum likelihood ancestral inference\n",
    "\n",
    "Ancestral reconstruction can be used to ask different questions about the evolution of related sequences. In some cases, we may be interested in the most likely sequence character at a specific node, without consideration of the assignments to all others. Such inference is referred to as a *marginal* reconstruction at the node in question. Conversely, we may want to determine the most likely joint assignment to all unknown nodes; that is, perform a *joint* reconstruction of the phylogenetic tree. This task demonstrates the latter.\n",
    "\n",
    "The `PhyloBNet` method `getMLJoint` performs joint reconstruction on the tree. The implemented algorithm takes a 'brute force' approach. This means that it exhaustively calculates the likelihood of all joint assignments to unknown nodes, and finds that which gives rise to the most likely network. A 'shortcut' version of the algorithm can be used to dramatically reduce the number of joint assignments tested. By default the `reduced` parameter, which instructs the algorithm to take this shortcut, is set to `True`. Take a look at the code for `getMLJoint`.\n",
    "\n",
    "**<span style=\"color:green\">Problem B7: In your own words, describe the shortcut which is taken when `getMLJoint` is run in 'reduced' mode.</span>**\n",
    "\n",
    "**<span style=\"color:green\">Problem B8: For the example network above, calculate the number of joint assignments which `getMLJoint` would test using both the reduced and non-reduced modes. </span>**\n",
    "\n",
    "Finally, we will perform joint reconstruction for our example network. As you should appreciate from Problem B8, this will take significantly less time than it would without running in 'reduced' mode!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "({'ancestor_1': 'A', 'ancestor_2': 'A'}, -13.978427822548813)\n"
     ]
    }
   ],
   "source": [
    "print(bn.getMLJoint())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will notice that two values are returned.\n",
    "* A dictionary indicating the most likely sequence character at each unknown node\n",
    "* The log-likelihood of the network when annotated with these characters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ancestral inference in the MBL protein superfamily"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will now perform some basic inference on MBL proteins. Candidates from some functionally characterised sub-families of the MBL superfamily were aligned, and a phylogenetic tree was subsequently inferred. The tree, including branch lengths, is shown below.\n",
    "\n",
    "<img src=\"mbl_tree.png\" width=\"800\" height=\"800\" />\n",
    "\n",
    "A contiguous 9-residue segment of the alignment was extracted and is shown in the figure below. It is known that some conserved residues in this site are involved in the coordination of metal ion cofactors. The equivalent sequences segments are provided to you in the file, `mbl_subseqs.fa`.\n",
    "\n",
    "<img src=\"subseq_aln.png\" width=\"600\" height=\"600\" />\n",
    "\n",
    "**On the basis of the tree and sub-sequence alignment provided, perform joint maximum likelihood ancestral reconstruction. Use the result(s) to infer the most likely 9-residue subsequence of these sequences' last common ancestor.**\n",
    "\n",
    "Some general pointers:\n",
    "* All nodes (*including internal*) require a unique label\n",
    "* The total log-likelihood of a subsequence reconstruction is simply the sum of the log-likelihood values for the individual columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cells provided here for your code. Add more as required\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:green\">Problem B9: Briefly describe your approach to the problem in words.</span>**\n",
    "\n",
    "**<span style=\"color:blue\">Problem B10: Provide the code you used to determine the maximum likelihood ancestral sequence.</span>**\n",
    "\n",
    "**<span style=\"color:red\">Problem B11: Provide your inferred subsequence of the last common ancestor. What is the  log-likelihood of this joint assignment of characters?</span>**\n",
    "\n",
    "Much about the true evolutionary history of the MBL superfamily is unknown to date due to its vast sequence and functional diversity. Research of this nature requires comprehensive sequence curation to ensure sufficient coverage of the complete phylogenetic space. Such coverage is not achieved by five representative sequences - in fact, the number of identified sub-families is at least triple that which is represented here. It is always important to consider the limitations - both statistical and biological - of any bioinformatics workflow.\n",
    "\n",
    "**<span style=\"color:green\">Problem B12: Explain how insufficient sampling of sequences in the superfamily will bias ancestral inference.</span>**\n",
    "\n",
    "As discussed above, an alternative approach to joint reconstruction is marginal reconstruction, where a single node is targetted by *marginalising* all others. Marginalisation in Bayesian networds was discussed in the week 4 material. These two reconstruction types can sometimes yield different inference results at a given node in a Bayesian network.\n",
    "\n",
    "**<span style=\"color:green\">Problem B13: Explain how the process of node marginalisation occurs by describing operations on a joint probability table.</span>**\n",
    "\n",
    "**<span style=\"color:purple\">Discussion Board Contributions: Provide a list of links to contributions you've made on the Ed discussion board.</span>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix: How-to install binfpy\n",
    "<a id='howto_install_binfpy'></a>\n",
    "\n",
    "\n",
    "The `binfpy` library is available from our local git server [GitLab_binfpy]. You can use git to store the binfpy directory on your computer where you can easily update it if any changes are made (see instructions below). You can also use the link to download the files in binfpy and place them in a folder of your choice. If any changes are made you will have to repeat the download process. **You will have to update your path regardless of which method you choose.**\n",
    "\n",
    "[GitLab_binfpy]: http://bioinf1.biosci.uq.edu.au/opensource/binfpy.git\n",
    "\n",
    "To install binfpy on your own computer, you need to have a git client installed on your computer or retrieve files from the web site above. \n",
    "\n",
    "**Mac OS X or Linux**\n",
    "\n",
    "If you're on Linux or Mac OS X, you should be set already. Open the terminal and change to a directory of choice and type \n",
    "```\n",
    "git clone http://bioinf1.biosci.uq.edu.au/opensource/binfpy.git\n",
    "```\n",
    "That will create a new directory called `binfpy` with a bunch of Python files.\n",
    "Add this directory, e.g. `/Users/johndoe/binfpy` to your PYTHONPATH, by adding the following line to your start-up file, e.g. `.profile`\n",
    "```\n",
    "export PYTHONPATH=/Users/johndoe/binfpy\n",
    "```\n",
    "This will be read next time you start a new shell, or you can activate immediately by \n",
    "```\n",
    "source .profile\n",
    "```\n",
    "\n",
    "**Windows**\n",
    "\n",
    "Install [git_for_windows]. Open the windows command prompt, navigate to a directory of your choice and type\n",
    "```\n",
    "git clone http://bioinf1.biosci.uq.edu.au/opensource/binfpy.git\n",
    "```\n",
    "That will create a new directory called `binfpy` with a bunch of Python files.\n",
    "Add this directory, e.g. `/Users/johndoe/binfpy` to your PYTHONPATH using the following instructions:\n",
    "\n",
    "Hit start and search for 'Environment variables'\n",
    "\n",
    "Add or edit the variable PYTHONPATH to include the `binfpy` directory\n",
    "\n",
    "[git_for_windows]: https://git-for-windows.github.io/\n",
    "\n",
    "**Updating**\n",
    "\n",
    "With all the above installed, you should be able to fire up your Python environment of choice, and `import` statements will be able to find the `binfpy` files. If there is an update, you can return to the binfpy directory and type\n",
    "```\n",
    "git pull\n",
    "```\n",
    "This will keep your files up-to-date. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<a id='export'></a>\n",
    "## How to export a single cell as a Python file\n",
    "If you are asked to save a cell as a Python file and upload it to Coder Quiz, you can do this easily by adding the line \n",
    "\n",
    "`%%writefile test.py`\n",
    "\n",
    "to the start of the Python cell, where test can be any name.\n",
    "\n",
    "Try it below to see how you can save this simple Python cell to your current working directory\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing test.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile test.py\n",
    "\n",
    "animal = \"dog\"\n",
    "\n",
    "if animal == \"dog\":\n",
    "    print (\"Bark\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}