Assignment1_answers.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sequence\n",
    "import prob\n",
    "import sym"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.05688622754491018\n",
      "0.057644892481508066\n"
     ]
    }
   ],
   "source": [
    "# Problem A1\n",
    "\n",
    "myDNA = sym.Alphabet('ACGT') # Define an alphabet, but many are pre-defined in sym.py, e.g. sym.Protein_Alphabet\n",
    "d = prob.Distrib(sequence.Protein_Alphabet) \n",
    "bacterial_sequences = sequence.readFastaFile('up_bac.fa')\n",
    "bg = prob.readDistrib('blosum62.distrib')\n",
    "\n",
    "# First part\n",
    "print(bg['S'])\n",
    "\n",
    "\n",
    "# Second part\n",
    "for bac_seq in bacterial_sequences:\n",
    "    for nucleotide in bac_seq:\n",
    "        d.observe(nucleotide)\n",
    "\n",
    "print(d['S'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'my_alpha' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[3], line 4\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[39m# Problem A2\u001b[39;00m\n\u001b[1;32m      2\u001b[0m \n\u001b[1;32m      3\u001b[0m \u001b[39m# Read sequences from the FASTA file 'example_seqsa.fa' using a specific alphabet 'my_alpha'\u001b[39;00m\n\u001b[0;32m----> 4\u001b[0m seqs \u001b[39m=\u001b[39m sequence\u001b[39m.\u001b[39mreadFastaFile(\u001b[39m'\u001b[39m\u001b[39mexample_seqsa.fa\u001b[39m\u001b[39m'\u001b[39m, my_alpha)\n\u001b[1;32m      6\u001b[0m \u001b[39m# Define the width of the motif we are looking for (this value needs to be assigned)\u001b[39;00m\n\u001b[1;32m      7\u001b[0m W \u001b[39m=\u001b[39m \u001b[39m7\u001b[39m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'my_alpha' is not defined"
     ]
    }
   ],
   "source": [
    "# Problem A2\n",
    "\n",
    "# Read sequences from the FASTA file 'example_seqsa.fa' using a specific alphabet 'my_alpha'\n",
    "seqs = sequence.readFastaFile('example_seqsa.fa', my_alpha)\n",
    "\n",
    "# Define the width of the motif we are looking for (this value needs to be assigned)\n",
    "W = 7\n",
    "\n",
    "# Initialize the Gibbs Motif discovery process for the sequences 'seqs' with the motif width 'W'\n",
    "g = gibbs.GibbsMotif(seqs, W)\n",
    "\n",
    "# Start the motif discovery process using the Gibbs algorithm for 1000 iterations\n",
    "q = g.discover(niter = 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LL @     0=\t97.29\n",
      "LL @   100=\t332.64\n",
      "LL @   200=\t1488.32\n",
      "LL @   300=\t2082.11\n",
      "LL @   400=\t2839.24\n",
      "LL @   500=\t3137.14\n",
      "LL @   600=\t3257.63\n",
      "LL @   700=\t3298.69\n",
      "LL @   800=\t3380.76\n",
      "LL @   900=\t3377.13\n",
      "LL @  1000=\t3376.79\n",
      "LL @  1100=\t3376.72\n",
      "LL @  1200=\t3378.45\n",
      "LL @  1300=\t3378.76\n",
      "LL @  1400=\t3379.86\n",
      "LL @  1500=\t3380.05\n",
      "LL @  1600=\t3375.67\n",
      "LL @  1700=\t3379.28\n",
      "LL @  1800=\t3376.35\n",
      "LL @  1900=\t3377.13\n",
      "Consensus: PROTEIN"
     ]
    }
   ],
   "source": [
    "# Problem A3\n",
    "import gibbs\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",
    "    \n",
    "    \n",
    "# Answer: PROTEIN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Background distribution: < A=0.04 B=0.04 C=0.04 D=0.04 E=0.04 F=0.04 G=0.04 H=0.04 I=0.04 J=0.04 K=0.04 L=0.03 M=0.04 N=0.04 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"
     ]
    }
   ],
   "source": [
    "# Problem A4\n",
    "p = g.getBackground()\n",
    "print('Background distribution:', p)\n",
    "# Answer: < A=0.04 B=0.04 C=0.04 D=0.04 E=0.04 F=0.04 G=0.04 H=0.04 I=0.04 J=0.04 K=0.04 L=0.03 M=0.04 N=0.04 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",
    "# The background was generated to be uniform across all amino acids\n",
    "# This is likely not ot happen in realy biological sequences as there are evolutionary pressures being applied to DNA, making the background\n",
    "# distribution to conform to certain, non-uniform probabilities\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A5\n",
    "# Answer: seq4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A6\n",
    "# Answer: 58-246"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A7\n",
    "# Answer: Ribonuclease Z/Hydroxyacylglutathione hydrolase-like \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'readFastaFile' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[14], line 7\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[39mfrom\u001b[39;00m \u001b[39mbinfpy\u001b[39;00m\u001b[39m.\u001b[39;00m\u001b[39mgibbs\u001b[39;00m \u001b[39mimport\u001b[39;00m GibbsMotif\n\u001b[1;32m      4\u001b[0m \u001b[39mfrom\u001b[39;00m \u001b[39mbinfpy\u001b[39;00m\u001b[39m.\u001b[39;00m\u001b[39msym\u001b[39;00m \u001b[39mimport\u001b[39;00m Protein_Alphabet\n\u001b[0;32m----> 7\u001b[0m seqs \u001b[39m=\u001b[39m readFastaFile(\u001b[39m'\u001b[39m\u001b[39mmbl_seqs.fa\u001b[39m\u001b[39m'\u001b[39m,Protein_Alphabet)\n\u001b[1;32m      8\u001b[0m W \u001b[39m=\u001b[39m \u001b[39m12\u001b[39m \u001b[39m# the width of the motif sought\u001b[39;00m\n\u001b[1;32m      9\u001b[0m \u001b[39m# Check for convergence\u001b[39;00m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'readFastaFile' is not defined"
     ]
    }
   ],
   "source": [
    "# Motif Discovery of MBLS\n",
    "from  pprint import pprint\n",
    "from binfpy.gibbs import GibbsMotif\n",
    "from binfpy.sym import Protein_Alphabet\n",
    "\n",
    "\n",
    "seqs = readFastaFile('mbl_seqs.fa',Protein_Alphabet)\n",
    "W = 12 # the width of the motif sought\n",
    "# Check for convergence\n",
    "\n",
    "g = GibbsMotif(seqs,W)\n",
    "q = g.discover(niter=4000)\n",
    "print('Consensus: ', end='')\n",
    "for pos in q:\n",
    "    print(pos.getmax(), end='')\n",
    "\n",
    "print()\n",
    "p = g.getBackground()\n",
    "print('Background distribution:', p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A8\n",
    "# Answer: The window sized used was 50 as this is the residues which are close to each other at the active site which are conserved.\n",
    "# By setting the window size to 50, it incorporates a balance between the foreground and background (signal to noise) and also allows for \n",
    "# variation around those conserved sites to account for the metal preference and/or catalytic function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A9\n",
    "# Answer: The number of iterations used was 4000 because upon each run of this niter paramater, \n",
    "# this contained the relevant histidine residues within the subsequence \"THLHGDH\" in the motif. \n",
    "# This pattern of histidiens being located close to each other was also showing up in niters lower  \n",
    "# than the selected one although less often, indicating that the gibbs algorithm was hovering around a likely global optima"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A\t -0.84  -1.61  -2.05  +0.32  -3.33  +0.15  -2.36  -2.36  +0.20  -0.09  +0.19  +0.59\n",
      "C\t +0.28  -5.49  -2.45  +0.22  -5.49  -0.70  -1.38  -1.38  +1.57  -0.09  -0.01  -0.41\n",
      "D\t -1.59  -2.88  -1.79  -0.52  +1.05  -0.25  +2.76  -2.19  -0.94  -0.16  +0.23  -0.66\n",
      "E\t -1.72  -2.92  -2.23  -0.29  -1.14  -0.25  -1.63  -2.23  -1.08  -0.81  -0.98  -1.95\n",
      "F\t -0.03  -1.46  -1.57  +0.45  -1.20  +0.43  -1.82  -1.68  +0.39  -0.33  -1.06  -0.08\n",
      "G\t -1.52  -1.52  -2.56  +0.56  -1.72  +0.05  -1.58  -1.41  -1.79  +1.70  +1.82  -0.72\n",
      "H\t -0.30  -0.90  +3.45  +0.55  +3.25  +0.03  -0.52  +3.44  -1.35  +0.17  -1.35  -0.41\n",
      "I\t +1.31  -1.76  -1.67  -0.98  -1.67  -0.16  -2.36  -1.45  +1.09  -0.98  -1.27  +0.33\n",
      "K\t -1.73  -2.20  -1.86  -1.73  -3.10  -0.48  -2.70  -1.17  -2.20  -1.41  -1.51  -1.61\n",
      "L\t +1.06  -1.39  -2.60  +0.02  -1.98  -0.25  -1.98  -2.12  -0.21  -0.59  -0.72  +1.06\n",
      "M\t -0.38  -2.97  -6.02  +0.11  -6.02  -0.53  -6.02  -6.02  +0.24  -0.25  -1.07  +0.35\n",
      "N\t +0.60  -1.90  -2.41  -0.09  -1.03  -0.03  -2.41  -2.12  -1.57  +0.39  +0.19  +0.73\n",
      "P\t -1.56  -1.65  -1.39  +0.59  -1.65  +0.79  -2.56  -1.65  -1.87  -1.18  -0.87  -1.01\n",
      "Q\t -1.81  -2.31  -1.22  -0.48  -1.22  -0.17  -3.38  -1.22  -2.71  -0.13  -0.78  -0.72\n",
      "R\t -2.47  -2.13  -1.78  -1.31  -2.97  -0.95  -1.68  -2.47  -1.37  -0.99  -1.52  -1.44\n",
      "S\t -0.66  +1.55  -2.17  +0.05  -1.41  +0.27  -1.56  -2.51  -0.15  +0.08  -0.39  -1.18\n",
      "T\t -1.20  +2.41  -1.43  +0.09  -1.83  -0.42  -1.61  -1.71  +0.65  -1.20  -1.08  -0.70\n",
      "V\t +0.52  -1.68  -1.15  -1.11  -1.94  -0.63  -2.04  -1.30  +0.76  -1.60  -1.68  +0.26\n",
      "W\t -2.66  -2.66  -1.99  +0.34  -5.71  +1.03  -1.60  -2.66  -0.51  -0.62  +0.07  -5.71\n",
      "Y\t -0.28  -2.70  -2.70  +0.47  -3.37  +0.68  -2.70  -2.30  +0.72  -0.42  -0.28  -0.93\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Construct PWM\n",
    "from binfpy.sequence import PWM\n",
    "\n",
    "pwm = PWM(q,p)\n",
    "print(pwm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A12\n",
    "# Answer: Background distribution: \n",
    "# < A=0.09 C=0.01 D=0.05 E=0.06 F=0.04 G=0.08 H=0.03 I=0.06 K=0.05 L=0.11 M=0.02 N=0.03 P=0.05 Q=0.03 R=0.06 S=0.06 T=0.05 V=0.07 W=0.02 Y=0.03  >\n",
    "# The background distribution is based off the gibbs sampling algorithm \n",
    "# The background distribution is the result of sampling a specific dataset that is \n",
    "# The reason why we did not use the background distribution based on a whole genome is because it would include too much noise\n",
    "# When its more general, its not representative of the domain we're trying to represent. These active sites exist within specific regions\n",
    "# of genomes and therefore the background distribution should be based on amino acids "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A13\n",
    "# Answer: The positions which stand out and have the highest level of conservation are the histidines that are within 2-3 AA's of each other\n",
    "# The reasons for their conservation are due to their involvement with the active site of the protein. \n",
    "# Histidine is the highest scoring residue indicating that it has particular importantce in the active site."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.099999999999987\n",
      "0.85\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "# Construct threshold discovery code\n",
    "positive_seqs = readFastaFile('positives.fa', Protein_Alphabet)\n",
    "negative_seqs = readFastaFile('negatives.fa', Protein_Alphabet)\n",
    "\n",
    "best_threshold = None\n",
    "best_combined_value = 0\n",
    "threshold = 0\n",
    "while threshold <= 10:\n",
    "    \n",
    "    # True Positives: Sequences in positive set that match\n",
    "    true_positives = sum([1 for sequence in positive_seqs if pwm.search(sequence, threshold)])\n",
    "    \n",
    "    # False Positives: Sequences in negative set that match\n",
    "    false_positives = sum([1 for sequence in negative_seqs if pwm.search(sequence, threshold)])\n",
    "    \n",
    "    false_negatives = 20 - true_positives # If there are less than 20 true positives, then false_negatives exist\n",
    "    \n",
    "    sensitivity = true_positives / 20 \n",
    "    \n",
    "    # Specificity calculation\n",
    "    specificity = 1 - (false_positives / 20)\n",
    "    \n",
    "    combined_value = sensitivity + specificity\n",
    "    if combined_value > best_combined_value:\n",
    "        best_combined_value = combined_value\n",
    "        best_threshold = threshold\n",
    "        \n",
    "    threshold += 0.1\n",
    "    \n",
    "print(best_threshold)\n",
    "print(sensitivity)\n",
    "print(specificity)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A14\n",
    "# Answer: The reason why this threshold was chosen is because it optimizes for having 20  matches in the positive.fa file\n",
    "# and 0 matches in the negatives.fa file.\n",
    "# it results in a sensitivity of 0.85 and a specifcity of 1\n",
    "# This threshold incorporates the scores needed from the conserved residues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A15\n",
    "# Answer: Given my threshold and observations, the likely sensitivity and specificity of the PWM are both high (close to 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of proteins 40\n",
      "Number of matches: 20\n",
      "Highest scoring protein: Q988B9\n",
      "Sequences identified by the motif\n",
      "sp|Q9KJB0|BLAB6_ELIME Metallo-beta-lactamase type 2 OS=Elizabethkingia meningoseptica OX=238 GN=blaB6 PE=3 SV=1\n",
      "tr|Q9K2N0|Q9K2N0_PSEAI Beta-lactamase class B VIM-2 OS=Pseudomonas aeruginosa OX=287 GN=blaVIM-2 PE=1 SV=1\n",
      "tr|Q5U807|Q5U807_PSEAI IMP-18 metallo-beta-lactamase OS=Pseudomonas aeruginosa OX=287 GN=IMP-18 PE=1 SV=2\n",
      "tr|Q7BJM5|Q7BJM5_PSEPU Class B beta-lactamase VIM-2 OS=Pseudomonas putida OX=303 GN=VIM-2 PE=1 SV=1\n",
      "sp|Q988B9|PDLA_RHILO 4-pyridoxolactonase OS=Mesorhizobium japonicum (strain LMG 29417 / CECT 9101 / MAFF 303099) OX=266835 GN=mlr6805 PE=1 SV=1\n",
      "sp|P26918|BLAB_AERHY Metallo-beta-lactamase type 2 OS=Aeromonas hydrophila OX=644 GN=cphA PE=1 SV=1\n",
      "tr|Q9K578|Q9K578_9GAMM FEZ-1 metallo-beta-lactamase OS=Fluoribacter gormanii OX=464 GN=blaFEZ-1 PE=1 SV=1\n",
      "tr|G9JVE6|G9JVE6_ECOLX NDM metallo-beta-lactamase OS=Escherichia coli OX=562 GN=blaNDM-6 PE=1 SV=1\n",
      "sp|P0A8V0|RBN_ECOLI Ribonuclease BN OS=Escherichia coli (strain K12) OX=83333 GN=rbn PE=1 SV=1\n",
      "tr|D1MEN9|D1MEN9_PSEAI Metallo-beta-lactamase VIM-2 OS=Pseudomonas aeruginosa OX=287 GN=blaVIM-2 PE=1 SV=1\n",
      "tr|Q8G9Q0|Q8G9Q0_PSEAI Metallo-beta-lactamase IMP-1 OS=Pseudomonas aeruginosa OX=287 GN=spm-1 PE=1 SV=1\n",
      "tr|X5CZV5|X5CZV5_KLEPN Metallo-beta-lactamase NDM-7 OS=Klebsiella pneumoniae OX=573 GN=blaNDM-7 PE=1 SV=1\n",
      "sp|A8ABX8|Y1254_IGNH4 UPF0173 metal-dependent hydrolase Igni_1254 OS=Ignicoccus hospitalis (strain KIN4/I / DSM 18386 / JCM 14125) OX=453591 GN=Igni_1254 PE=1 SV=1\n",
      "sp|O24495|GLO2M_ARATH Hydroxyacylglutathione hydrolase 1, mitochondrial OS=Arabidopsis thaliana OX=3702 GN=GLX2-1 PE=2 SV=2\n",
      "sp|P52700|BLA1_STEMA Metallo-beta-lactamase L1 type 3 OS=Stenotrophomonas maltophilia OX=40324 PE=1 SV=1\n",
      "tr|Q704V1|Q704V1_PSEAI GIM-1 metallo-beta-lactamase OS=Pseudomonas aeruginosa OX=287 GN=blaGIM-1 PE=1 SV=1\n",
      "sp|P14488|BLAB_BACCE Metallo-beta-lactamase type 2 OS=Bacillus cereus OX=1396 PE=1 SV=1\n",
      "tr|M1VE66|M1VE66_ECOLX Metallo-beta-lactamase NDM-8 OS=Escherichia coli OX=562 GN=blaNDM-8 PE=1 SV=1\n",
      "tr|A0A1L5BQA7|A0A1L5BQA7_SPHIB Subclass B3 metallo-beta-lactamase OS=Sphingobium indicum (strain DSM 16412 / CCM 7286 / MTCC 6364 / B90A) OX=861109 GN=SIDU_11385 PE=1 SV=1\n",
      "tr|B2FTM1|B2FTM1_STRMK Putative metallo-beta-lactamase l1 (Beta-lactamase type ii) (Ec 3.5.2.6) (Penicillinase) OS=Stenotrophomonas maltophilia (strain K279a) OX=522373 GN=Smlt2667 PE=1 SV=1\n",
      "\n",
      "Sequences not identified by the motif\n",
      "sp|A5U493|BLAC_MYCTA Beta-lactamase OS=Mycobacterium tuberculosis (strain ATCC 25177 / H37Ra) OX=419947 GN=blaC PE=1 SV=1\n",
      "tr|Q3XY33|Q3XY33_ENTFD Class C beta-lactamase OS=Enterococcus faecium (strain ATCC BAA-472 / TX0016 / DO) OX=333849 GN=HMPREF0351_10906 PE=4 SV=1\n",
      "tr|L7Y0F6|L7Y0F6_9BACT Class C beta-lactamase OS=uncultured bacterium OX=77133 GN=blaUPO01 PE=4 SV=1\n",
      "sp|P9WKD3|BLAC_MYCTU Beta-lactamase OS=Mycobacterium tuberculosis (strain ATCC 25618 / H37Rv) OX=83332 GN=blaC PE=1 SV=1\n",
      "tr|S3ZZU2|S3ZZU2_9ACTN Putative Beta-lactamase OS=Streptomyces aurantiacus JA 4570 OX=1286094 GN=STRAU_2968 PE=4 SV=1\n",
      "tr|F8V181|F8V181_9ZZZZ Putative beta-lactamase class C OS=uncultured microorganism OX=358574 PE=4 SV=1\n",
      "sp|P67918|BLA1_HAEIF Beta-lactamase ROB-1 OS=Haemophilus influenzae OX=727 GN=rob1 PE=3 SV=1\n",
      "tr|A0A1G9MHT8|A0A1G9MHT8_9FIRM Beta-lactamase class D OS=Romboutsia lituseburensis DSM 797 OX=1121325 GN=SAMN04515677_103211 PE=4 SV=1\n",
      "tr|A0A8F7JPT8|A0A8F7JPT8_9PSEU Class A beta-lactamase OS=Amycolatopsis sp. TNS106 OX=2861750 GN=CVV72_16010 PE=4 SV=1\n",
      "tr|A0A089NGP0|A0A089NGP0_9BACT Beta-lactamase class C OS=bacterium enrichment culture clone cep-03 OX=1548703 GN=blacep03 PE=4 SV=1\n",
      "sp|P54427|YBXI_BACSU Probable beta-lactamase YbxI OS=Bacillus subtilis (strain 168) OX=224308 GN=ybxI PE=1 SV=2\n",
      "tr|A0A0S2ZW34|A0A0S2ZW34_FUSNP Class D beta-lactamase OS=Fusobacterium nucleatum subsp. polymorphum OX=76857 GN=RN93_09935 PE=4 SV=1\n",
      "tr|A0A1C4A7I3|A0A1C4A7I3_9BACT Beta-lactamase class D OS=Chitinophaga costaii OX=1335309 GN=GA0116948_1022 PE=4 SV=1\n",
      "tr|Q66VX4|Q66VX4_9ZZZZ Class D beta-lactamase OS=Plasmid pB10-2 OX=291586 GN=oxa-2 PE=3 SV=1\n",
      "tr|A0A0Q8B7Q8|A0A0Q8B7Q8_9HYPH Class D beta-lactamase OS=Rhizobium sp. Root564 OX=1736563 GN=ASD74_08945 PE=4 SV=1\n",
      "tr|Q66VX4|Q66VX4_9ZZZZ Class D beta-lactamase OS=Plasmid pB10-2 OX=291586 GN=oxa-2 PE=3 SV=1\n",
      "tr|A0A8I1PET0|A0A8I1PET0_9PSEU Class A beta-lactamase OS=Pseudonocardia sp OX=60912 GN=bla PE=4 SV=1\n",
      "tr|A0A8F5NP02|A0A8F5NP02_NOCIO Class A beta-lactamase OS=Nocardia iowensis OX=204891 GN=bla PE=4 SV=1\n",
      "sp|Q5YXD6|BLAC_NOCFA Beta-lactamase FAR-1 OS=Nocardia farcinica (strain IFM 10152) OX=247156 GN=bla PE=1 SV=2\n",
      "sp|P0A9Z8|BLA2_KLEPN Beta-lactamase SHV-2 OS=Klebsiella pneumoniae OX=573 GN=bla PE=3 SV=1\n"
     ]
    }
   ],
   "source": [
    "# Write code for searching the beta-lactam degrading enzymes here...\n",
    "from sequence import readFastaFile\n",
    "active_seqs = readFastaFile('active.fa', Protein_Alphabet)\n",
    "# Search and store both the sequence and its best result\n",
    "results_with_seq = [(sequence, pwm.search(sequence, best_threshold)) for sequence in active_seqs]\n",
    "print(f\"Total number of proteins {len(active_seqs)}\")\n",
    "print(f\"Number of matches: {sum(len(matches[1]) for matches in results_with_seq)}\")\n",
    "\n",
    "# Extract the best score for each sequence\n",
    "best_score_with_seq = [(seq, max(res, key=lambda x: x[2])) for seq, res in results_with_seq if len(res) > 0]\n",
    "\n",
    "# Find the highest score across all sequences\n",
    "highest_score_with_seq = max(best_score_with_seq, key=lambda x: x[1][2])\n",
    "\n",
    "print(f\"Highest scoring protein: {highest_score_with_seq[0].name}\")\n",
    "\n",
    "print(\"Sequences identified by the motif\")    \n",
    "for seq,score in best_score_with_seq:\n",
    "    print(seq.info)\n",
    "   \n",
    "print() \n",
    "    \n",
    "print(\"Sequences not identified by the motif\")    \n",
    "not_found = set(sequence.name for sequence in active_seqs) - set(seq.name for seq, score in best_score_with_seq)\n",
    "for sequence in active_seqs:\n",
    "    if sequence.name in not_found:\n",
    "        print(sequence.info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A17\n",
    "# Answer: 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A18\n",
    "# Answer: IPR001279"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A20\n",
    "# Answer: This is because the proteins that don't show up belong to class C or D proteins within the superfamily. Class C and D use serine amino acids\n",
    "# whilst class B use metal ions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A21\n",
    "# sp|A8ABX8|Y1254_IGNH4 UPF0173 metal-dependent hydrolase Igni_1254 OS=Ignicoccus hospitalis (strain KIN4/I / DSM 18386 / JCM 14125) OX=453591 GN=Igni_1254 PE=1 SV=1\n",
    "# This protein has promiscuous β-lactamase activity against Penicillins, Cephalosporins, Carbapenems but is it unknown for Chromogenic\n",
    "# The reason why my motif exists within this protein is because of the likely presence of the histidine residues that account for \n",
    "# the metal binding activity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem A22\n",
    "# Currently, the data that the gibbs sampling model is being trained on is all of the MBL domains. We would need to sanitize/filter\n",
    "# the the mbl_seqs.fa such that it only contains MBL domains from class B/true beta lactamases\n",
    "# We could go a step futher and develop multiple gibbs sampling models that are trained to identify other MBLS within the superfamily which are each \n",
    "# trained on sequence data from their own MBL domain so the motifs from each domain can be accurately identified. \n",
    "# We can cross-validate the models on each data on each of their own positive and negative data sets set to ensure the models are performing accurately. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question B1\n",
    "Answer: $P(\\alpha_A,\\alpha_B,\\alpha_C,\\alpha_D,\\alpha_X,\\alpha_Y,\\alpha_Z) = P(\\alpha_Z) \\times P(\\alpha_X|\\alpha_Z) \\times P(\\alpha_Y|\\alpha_Z) \\times P(\\alpha_A|\\alpha_X) \\times P(\\alpha_B|\\alpha_X) \\times P(\\alpha_C|\\alpha_Y) \\times P(\\alpha_D|\\alpha_Y)$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question B2\n",
    "Answer: We can assume conditional independence for ancestral sequences because we assume that the evolution of nucleotide/amino acids that diverged from a common ancestor are not dependent on each others evolution after their divergence. Their probability for what they are is only dependent on the common ancestor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question B3\n",
    "Answer: Although the statistical simplification of column independence reduces the number of free parameters to handle when calculating ancestor states, from a biological perspective, nucleotides/amino acids occur together within a physical space and the closer they are together, they are more likely to be subjected to the same physiochemical alterations. A mutation at one site might be compensated by a mutation at another site to to maintain the protein's function or structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question B4\n",
    "Answer: R(N) = (N * (N-1)) / 2\n",
    "\n",
    "R(20) = (20*19)/2 = 190"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question B5\n",
    "Answer: the biological reason for the differences in transition probabilities comes down to the physiochemical similarity between the amino acids. Amino acids are more likely to make transitions which are more tolerable and thus more probable. We also need to take into account codon redundancy, and within this context it doesn't take much change to the codonn to transition from An alanine to a valine, where as it takes two nucleotide changes to change from an alanine to a leucine. Its even less likely to transition to a tyrosine as there are only two possible codon combinations that result in a tyrosine as the others result in a stop codon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question B6\n",
    "The transition probabilities are decreasing when the branch length increases for the same amino acid, indicating an inversely proportional relationship between the inputs and the outputs of the function. This is because over long evolutionary timescales, the liklihood of mutations occuring increases. The probabilities tend towards the equilibrium frequencies of the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question B7\n",
    "Answer: The getMLJoint method performs joint reconstruction on the tree. In 'reduced' mode, the algorithm takes a 'shortcut' to reduce the number of joint assignments tested by skipping over certain joint assignments that are deemed suboptimal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question B8\n",
    "Answer: We need to consider the size of the amino acid alphabet and the number of internal nodes.\n",
    "\n",
    "Non-reduced = (20x20)x2 = 800\n",
    "\n",
    "Reduced = (3x3)x2 = 18"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question B9\n",
    "Answer: To infer the ancestral sequence of the last common ancestor, we'll utilize the maximum likelihood ancestral reconstruction approach. This will involve calculating the likelihood of each possible ancestral sequence given the observed sequences at the leaves and the phylogenetic tree.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "column_data: ['T', 'S', 'T', 'T', 'N']\n",
      "ml_joint: ({'root': 'T', 'X': 'T', 'Y': 'T', 'A': 'T'}, -13.05197744602148)\n",
      "column_data: ['H', 'H', 'H', 'D', 'H']\n",
      "ml_joint: ({'root': 'H', 'X': 'H', 'Y': 'H', 'A': 'H'}, -10.014371623910987)\n",
      "column_data: ['H', 'E', 'L', 'G', 'G']\n",
      "ml_joint: ({'root': 'G', 'X': 'G', 'Y': 'G', 'A': 'G'}, -25.77353217273115)\n",
      "column_data: ['H', 'H', 'H', 'H', 'H']\n",
      "ml_joint: ({'root': 'H', 'X': 'H', 'Y': 'H', 'A': 'H'}, -4.228712363558731)\n",
      "column_data: ['Y', 'F', 'P', 'Y', 'Y']\n",
      "ml_joint: ({'root': 'Y', 'X': 'Y', 'Y': 'Y', 'A': 'Y'}, -15.752671367560556)\n",
      "column_data: ['D', 'D', 'D', 'D', 'D']\n",
      "ml_joint: ({'root': 'D', 'X': 'D', 'Y': 'D', 'A': 'D'}, -3.1897325692530076)\n",
      "column_data: ['H', 'H', 'H', 'V', 'I']\n",
      "ml_joint: ({'root': 'H', 'X': 'V', 'Y': 'H', 'A': 'H'}, -18.223339738237687)\n",
      "column_data: ['T', 'A', 'V', 'P', 'F']\n",
      "ml_joint: ({'root': 'A', 'X': 'A', 'Y': 'A', 'A': 'A'}, -26.135781865470523)\n",
      "column_data: ['G', 'G', 'C', 'G', 'G']\n",
      "ml_joint: ({'root': 'G', 'X': 'G', 'Y': 'G', 'A': 'G'}, -10.457783546215055)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('THGHYDHAG', -126.82790269295917)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Question B10\n",
    "\n",
    "import asr\n",
    "\n",
    "# Function to perform joint maximum likelihood ancestral reconstruction for a single column\n",
    "def infer_ancestral_for_column(column_data):\n",
    "    # Starting with root node\n",
    "    root = asr.PhyloBNode(JTT, label='root')\n",
    "\n",
    "    # The way column_data is index is according to the order in which the protein type occurs \n",
    "    internal_node1 = asr.PhyloBNode(JTT, parent=root, distance=0.1022, label='X') \n",
    "    leaf1 = asr.PhyloBNode(JTT, parent=internal_node1, distance=0.5378, label='RNazeZ', annot=column_data[3])\n",
    "    leaf2 = asr.PhyloBNode(JTT, parent=internal_node1, distance=1.3713, label='FDP', annot=column_data[4])\n",
    "\n",
    "    internal_node2 = asr.PhyloBNode(JTT, parent=root, distance=0.1374, label='Y')\n",
    "    leaf3 = asr.PhyloBNode(JTT, parent=internal_node2, distance=0.6324, label='Glyoxalase', annot=column_data[0])\n",
    "\n",
    "    internal_node3 = asr.PhyloBNode(JTT, parent=root, distance=0.178, label='A')\n",
    "    leaf4 = asr.PhyloBNode(JTT, parent=internal_node3, distance=0.6106, label='MBL_Class3', annot=column_data[1])\n",
    "    leaf5 = asr.PhyloBNode(JTT, parent=internal_node3, distance=1.1581, label='Lactonase', annot=column_data[2])\n",
    "\n",
    "    # Construct the Bayesian network\n",
    "    bn = asr.PhyloBNet(root=root)\n",
    "    bn.addNodes([internal_node1, internal_node2,internal_node3, leaf1, leaf2, leaf3, leaf4, leaf5])\n",
    "\n",
    "    # Perform joint maximum likelihood ancestral reconstruction\n",
    "    return bn.getMLJoint()\n",
    "\n",
    "# Iterating over each column of the alignment\n",
    "ancestral_sequence = []\n",
    "log_likelihoods = []\n",
    "sequences = sequence.readFastaFile('mbl_subseqs.fa')\n",
    "JTT = asr.MODELS['JTT']\n",
    "for col_idx in range(sequences[0].length):\n",
    "    column_data = [seq.sequence[col_idx] for seq in sequences]\n",
    "    print(f\"column_data: {column_data}\")\n",
    "    ml_joint = infer_ancestral_for_column(column_data)\n",
    "    print(f\"ml_joint: {ml_joint}\")\n",
    "    ancestral_sequence.append(ml_joint[0]['root'])  # Appending the inferred character for the root\n",
    "    log_likelihoods.append(ml_joint[1])\n",
    "\n",
    "total_log_likelihood = sum(log_likelihoods)\n",
    "''.join(ancestral_sequence), total_log_likelihood\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem B11\n",
    "Answer: THGHYDHAG, −126.69"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question B12\n",
    "Answer: Insufficient sampling of sequences in the superfamily will lead to a skewed representation. If certain subfamilies are overrepresented or underrepresented in the sample, it can bias the probability estimates, leading to potentially inaccurate ancestral sequence reconstructions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question B13:\n",
    "Answer: Node marginalization in the context of a joint probability table involves summing over the probabilities of all possible states of a certain variable, effectively 'eliminating' or 'averaging out' that variable from consideration. This allows us to obtain a probability distribution for the remaining variables without having to specify a particular state for the marginalized variable.\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}