3_blast_msa.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.7.6"
    },
    "colab": {
      "name": "3_blast_msa.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Intertangler/bioinformatics_stockholm/blob/master/3_blast_msa.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Opgj2j-E9E7_"
      },
      "source": [
        "# BLAST"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0vJOFI11GSZH"
      },
      "source": [
        "In this exercise, you will get to try some BLAST and multiple sequence alignment from the context of a Python program. The browser-based web tools are great resources, and we will provide links to these in parallel with the programming exercises. However, the goal here is to experience some of the flexibility that comes with being able to put together our own tools to do  specific tasks and answer particular questions that might not perfectly overlap with the tools provided on websites.\n",
        "\n",
        "We will use a database called Swiss-Prot, that was created in 1986 by Amos Bairoch during his PhD and developed by the Swiss Institute of Bioinformatics. Its aim was to provide reliable protein sequences and a minimal level of redundancy. It is a good source for retrieving high quality records. Today it is integrated with Uni-Prot.\n",
        "\n",
        "First run the two cells below to install biopython and download the entire Swiss-Prot database onto Google's server. Other databases can be downloaded in this same way through NCBI's file transfer page: https://ftp.ncbi.nlm.nih.gov/blast/db/ Including the entire genbank nucleotide database for example as well as specialized databases for certain genomes, viruses, etc..\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PSH-pmrswRwp"
      },
      "source": [
        "pip install biopython"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vhbtVG0effHV"
      },
      "source": [
        "!wget https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/swissprot.gz"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pMN6TJSZLA7x"
      },
      "source": [
        "Run the line below to unzip the compressed database file. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ewax8T8pgFT7"
      },
      "source": [
        "!gunzip swissprot.gz"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IBdgBycmLTv1"
      },
      "source": [
        "Next run the cell below to download two home-made modules that we will use in this exercise; they are based mostly on code from Migeul Rocha https://github.com/miguelfrocha/BioinformaticsAlgorithmsBook but we have added a few things for the purposes of visualization and the goals of this exercise. The first module will help us perform multiple sequence alignments and visualize them, and the second is a simplified BLAST program (which does not utilize Blosum or a substitution matrix - only a fixed match/mismatch score) that we will be able to run locally."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Nx_8jqHXe_kP"
      },
      "source": [
        "import sys\n",
        "from os import path\n",
        "\n",
        "if path.exists(\"/content/bioinformatics_stockholm/\"):\n",
        "  pass\n",
        "else:\n",
        "  !git clone https://github.com/Intertangler/bioinformatics_stockholm\n",
        "sys.path.insert(0,'/content')\n",
        "\n",
        "import importlib.util\n",
        "spec = importlib.util.spec_from_file_location(\"msa\", \"/content/bioinformatics_stockholm/msa/msa.py\")\n",
        "msa = importlib.util.module_from_spec(spec)\n",
        "spec.loader.exec_module(msa)\n",
        "\n",
        "spec = importlib.util.spec_from_file_location(\"miniblast\", \"/content/bioinformatics_stockholm/Miniblast/miniblast.py\")\n",
        "miniblast = importlib.util.module_from_spec(spec)\n",
        "spec.loader.exec_module(miniblast)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7W_83x4PNkN8"
      },
      "source": [
        "## web BLAST\n",
        "First, if you have never used it before, check out the web-based BLAST hosted by the NIH. https://blast.ncbi.nlm.nih.gov/Blast.cgi Note the multiple types of BLAST that you can choose - for example nucleotide or protein-based depending on the type of sequence you will use for your query. If you go to protein blast, you will see an input box that accepts sequences or accession numbers. Note the various options such as choice of database (the standard one for protein searches is called \"nr\"). You can try searching for this sequence for example: \"MKWVTFISLLFLFSSAYSRGVFRRDAHKSEVAHRFKDLGEENFKALVLIAFAQYLQQCPFEDHVKLVNEVTEFAKTCVADESAENCDKSLHTLFGDKLCTVATLRETYGEMADCCAKQEPERNECFLQHKDDNPNLPRLVRPEVDVMCTAFHDNEETFLKKYLYEIA\". You should get several hits related to the protein albumin. Their significance level should be high, indicated by correspondingly low E-values (expected number of random matches to your input)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r39PgrBdccnA"
      },
      "source": [
        "## local BLAST\n",
        "Now that we've tried the internet-based version, let's do some local BLASTing on the Swiss-Prot database that we have downloaded. First run the cell below to define our Blosum62 matrix that will be used for our alignments."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3iZgXZVOomxe"
      },
      "source": [
        "blosum62 = {\n",
        "    ('W', 'F'): 1, ('L', 'R'): -2, ('S', 'P'): -1, ('V', 'T'): 0,\n",
        "    ('Q', 'Q'): 5, ('N', 'A'): -2, ('Z', 'Y'): -2, ('W', 'R'): -3,\n",
        "    ('Q', 'A'): -1, ('S', 'D'): 0, ('H', 'H'): 8, ('S', 'H'): -1,\n",
        "    ('H', 'D'): -1, ('L', 'N'): -3, ('W', 'A'): -3, ('Y', 'M'): -1,\n",
        "    ('G', 'R'): -2, ('Y', 'I'): -1, ('Y', 'E'): -2, ('B', 'Y'): -3,\n",
        "    ('Y', 'A'): -2, ('V', 'D'): -3, ('B', 'S'): 0, ('Y', 'Y'): 7,\n",
        "    ('G', 'N'): 0, ('E', 'C'): -4, ('Y', 'Q'): -1, ('Z', 'Z'): 4,\n",
        "    ('V', 'A'): 0, ('C', 'C'): 9, ('M', 'R'): -1, ('V', 'E'): -2,\n",
        "    ('T', 'N'): 0, ('P', 'P'): 7, ('V', 'I'): 3, ('V', 'S'): -2,\n",
        "    ('Z', 'P'): -1, ('V', 'M'): 1, ('T', 'F'): -2, ('V', 'Q'): -2,\n",
        "    ('K', 'K'): 5, ('P', 'D'): -1, ('I', 'H'): -3, ('I', 'D'): -3,\n",
        "    ('T', 'R'): -1, ('P', 'L'): -3, ('K', 'G'): -2, ('M', 'N'): -2,\n",
        "    ('P', 'H'): -2, ('F', 'Q'): -3, ('Z', 'G'): -2, ('X', 'L'): -1,\n",
        "    ('T', 'M'): -1, ('Z', 'C'): -3, ('X', 'H'): -1, ('D', 'R'): -2,\n",
        "    ('B', 'W'): -4, ('X', 'D'): -1, ('Z', 'K'): 1, ('F', 'A'): -2,\n",
        "    ('Z', 'W'): -3, ('F', 'E'): -3, ('D', 'N'): 1, ('B', 'K'): 0,\n",
        "    ('X', 'X'): -1, ('F', 'I'): 0, ('B', 'G'): -1, ('X', 'T'): 0,\n",
        "    ('F', 'M'): 0, ('B', 'C'): -3, ('Z', 'I'): -3, ('Z', 'V'): -2,\n",
        "    ('S', 'S'): 4, ('L', 'Q'): -2, ('W', 'E'): -3, ('Q', 'R'): 1,\n",
        "    ('N', 'N'): 6, ('W', 'M'): -1, ('Q', 'C'): -3, ('W', 'I'): -3,\n",
        "    ('S', 'C'): -1, ('L', 'A'): -1, ('S', 'G'): 0, ('L', 'E'): -3,\n",
        "    ('W', 'Q'): -2, ('H', 'G'): -2, ('S', 'K'): 0, ('Q', 'N'): 0,\n",
        "    ('N', 'R'): 0, ('H', 'C'): -3, ('Y', 'N'): -2, ('G', 'Q'): -2,\n",
        "    ('Y', 'F'): 3, ('C', 'A'): 0, ('V', 'L'): 1, ('G', 'E'): -2,\n",
        "    ('G', 'A'): 0, ('K', 'R'): 2, ('E', 'D'): 2, ('Y', 'R'): -2,\n",
        "    ('M', 'Q'): 0, ('T', 'I'): -1, ('C', 'D'): -3, ('V', 'F'): -1,\n",
        "    ('T', 'A'): 0, ('T', 'P'): -1, ('B', 'P'): -2, ('T', 'E'): -1,\n",
        "    ('V', 'N'): -3, ('P', 'G'): -2, ('M', 'A'): -1, ('K', 'H'): -1,\n",
        "    ('V', 'R'): -3, ('P', 'C'): -3, ('M', 'E'): -2, ('K', 'L'): -2,\n",
        "    ('V', 'V'): 4, ('M', 'I'): 1, ('T', 'Q'): -1, ('I', 'G'): -4,\n",
        "    ('P', 'K'): -1, ('M', 'M'): 5, ('K', 'D'): -1, ('I', 'C'): -1,\n",
        "    ('Z', 'D'): 1, ('F', 'R'): -3, ('X', 'K'): -1, ('Q', 'D'): 0,\n",
        "    ('X', 'G'): -1, ('Z', 'L'): -3, ('X', 'C'): -2, ('Z', 'H'): 0,\n",
        "    ('B', 'L'): -4, ('B', 'H'): 0, ('F', 'F'): 6, ('X', 'W'): -2,\n",
        "    ('B', 'D'): 4, ('D', 'A'): -2, ('S', 'L'): -2, ('X', 'S'): 0,\n",
        "    ('F', 'N'): -3, ('S', 'R'): -1, ('W', 'D'): -4, ('V', 'Y'): -1,\n",
        "    ('W', 'L'): -2, ('H', 'R'): 0, ('W', 'H'): -2, ('H', 'N'): 1,\n",
        "    ('W', 'T'): -2, ('T', 'T'): 5, ('S', 'F'): -2, ('W', 'P'): -4,\n",
        "    ('L', 'D'): -4, ('B', 'I'): -3, ('L', 'H'): -3, ('S', 'N'): 1,\n",
        "    ('B', 'T'): -1, ('L', 'L'): 4, ('Y', 'K'): -2, ('E', 'Q'): 2,\n",
        "    ('Y', 'G'): -3, ('Z', 'S'): 0, ('Y', 'C'): -2, ('G', 'D'): -1,\n",
        "    ('B', 'V'): -3, ('E', 'A'): -1, ('Y', 'W'): 2, ('E', 'E'): 5,\n",
        "    ('Y', 'S'): -2, ('C', 'N'): -3, ('V', 'C'): -1, ('T', 'H'): -2,\n",
        "    ('P', 'R'): -2, ('V', 'G'): -3, ('T', 'L'): -1, ('V', 'K'): -2,\n",
        "    ('K', 'Q'): 1, ('R', 'A'): -1, ('I', 'R'): -3, ('T', 'D'): -1,\n",
        "    ('P', 'F'): -4, ('I', 'N'): -3, ('K', 'I'): -3, ('M', 'D'): -3,\n",
        "    ('V', 'W'): -3, ('W', 'W'): 11, ('M', 'H'): -2, ('P', 'N'): -2,\n",
        "    ('K', 'A'): -1, ('M', 'L'): 2, ('K', 'E'): 1, ('Z', 'E'): 4,\n",
        "    ('X', 'N'): -1, ('Z', 'A'): -1, ('Z', 'M'): -1, ('X', 'F'): -1,\n",
        "    ('K', 'C'): -3, ('B', 'Q'): 0, ('X', 'B'): -1, ('B', 'M'): -3,\n",
        "    ('F', 'C'): -2, ('Z', 'Q'): 3, ('X', 'Z'): -1, ('F', 'G'): -3,\n",
        "    ('B', 'E'): 1, ('X', 'V'): -1, ('F', 'K'): -3, ('B', 'A'): -2,\n",
        "    ('X', 'R'): -1, ('D', 'D'): 6, ('W', 'G'): -2, ('Z', 'F'): -3,\n",
        "    ('S', 'Q'): 0, ('W', 'C'): -2, ('W', 'K'): -3, ('H', 'Q'): 0,\n",
        "    ('L', 'C'): -1, ('W', 'N'): -4, ('S', 'A'): 1, ('L', 'G'): -4,\n",
        "    ('W', 'S'): -3, ('S', 'E'): 0, ('H', 'E'): 0, ('S', 'I'): -2,\n",
        "    ('H', 'A'): -2, ('S', 'M'): -1, ('Y', 'L'): -1, ('Y', 'H'): 2,\n",
        "    ('Y', 'D'): -3, ('E', 'R'): 0, ('X', 'P'): -2, ('G', 'G'): 6,\n",
        "    ('G', 'C'): -3, ('E', 'N'): 0, ('Y', 'T'): -2, ('Y', 'P'): -3,\n",
        "    ('T', 'K'): -1, ('A', 'A'): 4, ('P', 'Q'): -1, ('T', 'C'): -1,\n",
        "    ('V', 'H'): -3, ('T', 'G'): -2, ('I', 'Q'): -3, ('Z', 'T'): -1,\n",
        "    ('C', 'R'): -3, ('V', 'P'): -2, ('P', 'E'): -1, ('M', 'C'): -1,\n",
        "    ('K', 'N'): 0, ('I', 'I'): 4, ('P', 'A'): -1, ('M', 'G'): -3,\n",
        "    ('T', 'S'): 1, ('I', 'E'): -3, ('P', 'M'): -2, ('M', 'K'): -1,\n",
        "    ('I', 'A'): -1, ('P', 'I'): -3, ('R', 'R'): 5, ('X', 'M'): -1,\n",
        "    ('L', 'I'): 2, ('X', 'I'): -1, ('Z', 'B'): 1, ('X', 'E'): -1,\n",
        "    ('Z', 'N'): 0, ('X', 'A'): 0, ('B', 'R'): -1, ('B', 'N'): 3,\n",
        "    ('F', 'D'): -3, ('X', 'Y'): -1, ('Z', 'R'): 0, ('F', 'H'): -1,\n",
        "    ('B', 'F'): -3, ('F', 'L'): 0, ('X', 'Q'): -1, ('B', 'B'): 4\n",
        "}\n",
        "Blosum62 = {}\n",
        "for i in blosum62:\n",
        "  Blosum62[i] = blosum62[i]\n",
        "  Blosum62[i[::-1]] = blosum62[i] #adds the reverse pairs so that we can lookup ('F', 'H') or ('H', 'F') for example"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CIJj3xbyd3EK"
      },
      "source": [
        "Next we will try BLASTing for the same albumin sequence we used in the web-BLAST. \"MKWVTFISLLFLFSSAYSRGVFRRDAHKSEVAHRFKDLGEENFKALVLIAFAQYLQQCPFEDHVKLVNEVTEFAKTCVADESAENCDKSLHTLFGDKLCTVATLRETYGEMADCCAKQEPERNECFLQHKDDNPNLPRLVRPEVDVMCTAFHDNEETFLKKYLYEIA\". \n",
        "The cell below will first run a function miniblast.read_database_fasta() which will access the lines of the Swiss-Prot database residing on the server. Then the miniblast.best_alignment_fasta() function will execute the search, with a query sequence, the database, and the word- or k-mer size as input parameters. The resulting output will be attached to a variable we call result, and in a subsequent loop we can print out the contents of the result."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WVUgYJPDd9Ul"
      },
      "source": [
        "from Bio import SeqIO\n",
        "from Bio.Seq import Seq\n",
        "from Bio import pairwise2\n",
        "import numpy as np\n",
        "db = miniblast.read_database_fasta('/content/swissprot')\n",
        "query = \"MKWVTFISLLFLFSSAYSRGVFRRDAHKSEVAHRFKDLGEENFKALVLIAFAQYLQQCPFEDHVKLVNEVTEFAKTCVADESAENCDKSLHTLFGDKLCTVATLRETYGEMADCCAKQEPERNECFLQHKDDNPNLPRLVRPEVDVMCTAFHDNEETFLKKYLYEIA\"\n",
        "result = miniblast.best_alignment_fasta(db,query,w=3)\n",
        "print(result)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WnULxicc1CAe"
      },
      "source": [
        "The database object db will return the result of an entry if the position of the entry i is used as an index."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uUVe9kGMht1z"
      },
      "source": [
        "for i in result:\n",
        "  print(i)\n",
        "  print(db[i])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t7z4_qhLgdMs"
      },
      "source": [
        "Note how we get some additional hits apparently unrelated to our query (which we already know to be albumin). To better understand which hits are genuine, it is always important to examine the alignments directly. This miniblast module does not perform the final Smith-Waterman local alignment typical in the BLAST workflow, so we can again make use of biopython's pairwise2.align() function in order to perform the final local alignment on the winning results with our query. We will simply set up a loop that goes through the database, looking up the indices produced by the BLAST function above, and which performs alignment with our query and formatting to visually view the alignments like we did with the NGS data last time.\n",
        "\n",
        "If we check the relative scores and visually inspect each alignment now, we see that albumin (indicated where it says RecName: Full=Albumin) is indeed our winner, as the other candidates have only sparse agreement with a lot of gaps. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aRBh6TvSgcus"
      },
      "source": [
        "for i in result:\n",
        "  print(db[i])\n",
        "  target_sequence = db[i].seq\n",
        "  aln = pairwise2.align.localds(query,target_sequence, match_dict = Blosum62, open = -3, extend = -1)[0]\n",
        "  print(pairwise2.format_alignment(align1=aln[0], align2=aln[1], score=aln[2], begin=aln[3], end=aln[4], full_sequences=True))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pgkkgpTPjIEA"
      },
      "source": [
        "## Problem 3.1\n",
        "\n",
        "Consider a hypothetical scenario: the following sequences are fragmented data that have been isolated from a pathogen whose current identity is unknown. It is believed that the sequences are all derived from a single protein related to the pathogen. To make things interesting, we've tossed in a few totally random sequences to simulate a bit of experimental noise. Thus, you could go through all the sequences and web-BLAST them individually, but a more convenient approach will be to write a program that systematically performs miniblast on our Swiss-Prot database, and which at the same time generates alignments for you to check as we did above. After running the program, visually inspect the results or locate the high-scoring alignments and determine the protein of interest. **What is the protein's name? For your submission, report the answer as the name of the entry, i.e. where it says RecName: Full=.... in the database record** Did you discover the identity of the pathogen? Note that the BLAST function can take a few minutes to run with multiple queries.\n",
        "\n",
        "fragments = [\"EHFTGQLPMKQVTATIMIRNCEMEIRHIWCMREIVFEETQGEM\",\"VSWEDPLVISTDANFCGNHWYPSTQDHEACNKF\",\"FWMLMIYNVAGNKWVTVYYGVPVWKEAKTTLFCASDAKGY\",\"TTEYFGKGIPLYRKFEWSDPPSPFSAKQTYHWGKMTTSYYFDH\",\"EREVHNIWATHACVPTDPSPQEMFLVNVTENFNMWKNDMVDQMHEDIISL\",\"THGIKPVVSTQLLLNGSLAEEEIIIRSENLTDNAKTIIVHLNESIEIECVRPNNNTRKSVRIGPGQTFYATGDIIGDIRAAYCDINKDKWNRALQRVSKKLRELYPNTTTIH\",\"DKWNRALQRVSKKLRELYPNTTTIHFNSSSGGDPEITMHSFNCRGEFFYCDTSRLFNGTYNGTDSNPNNGTDSDSNSNSNITLPCRIKQIINMWQEVGRAMYAPPIAGSITCKSNITGLLLTRDGGKEENSTTNRIEIFRPGGGNMKDN\",\"LSQNQQEQNEKDLLSLDKWQTLWNWFDITNWLWYIKIFIMIVGGLIGLRIIFAVLSIINRVRQGYSPLSLQTLIPNPREPDRLGRIEEEGGEQDRNRSIRLVNGFLALAWDDLRNLCLFSYHQLRNFILIVARAVELLGRSSLRGLQRGWEAL\", \"FLALAWDDLRNLCLFSYHQLRNFILIVARAVELLGRSSLRGLQRGWEALKYLGSIVQYWGLELKKSAVSLFDAIAIAVAEGTDRII\", \"FKGRWKESSMPYFWMFQHADTSIIIWFCKMWCPTRTTSFKT\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B-EOYVqH86oS"
      },
      "source": [
        "# multiple sequence alignment\n",
        "\n",
        "The code below will performa multiple sequence alignment on a dictionary of input sequences (with their names or ID's as keys, and their sequences as values). Try running it for the test sequences below and have a look at the output alignment object that is produced, as well as the color interactive visualization we have made below that. The next problem will simply be a matter of replacing the input sequences and rerunning this program.\n",
        "\n",
        "In case you are curious, these sequences are from the zinc finger family that contains members involved in the regulation of transcription processes. Zinc finger domains play a role in the mediation of protein-DNA interactions. The multiple sequence alignment reveals to us some of the common/conserved structure necessary for this function."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "K50uKU4IVGBQ"
      },
      "source": [
        "from Bio import AlignIO\n",
        "from Bio.Align import MultipleSeqAlignment\n",
        "\n",
        "example_seqs = {\n",
        "\"ATRX_HUMAN/159-296\" :     \"KRGEDGLHGIVSCTACGQQVNHFQKDSIYRHPSLQVLICKNCFKYYMSDDISRDSDGMDEQCRWCAEGGNLICCDFCHNAFCKKCILRNLGRKELSTIMDENNQWYCYICHPEPLLDLVTACNSVFENLEQLLQQNKK\",\n",
        "\"ATRX_MACEU/21-158\"  :     \"KKRGEGLHGIVSCTACGQQVNHFQKDSIYRHPTLKVLICKNCYKYYMSDDISRDADGMDEQCRWCAEGGNLICCDFCHNAFCKKCILRNLGRKELSAIMDENSQWYCYICRPEPLLDLVTACHSVFKNLEQLLQQNKK\",\n",
        "\"ATRX_MOUSE/158-295\" :     \"KRGEDGLHGIVSCTACGQQVNHFQKDSIYRHPSLKVLICKNCFKYYMSDDISRDSDGMDEQCRWCAEGGNLICCDFCHNAFCKKCILRNLGRKELSTIMDENNQWYCYICQPEPLLDLVTACNSVFENLEQLLQQNKK\",\n",
        "\"ATRX_PANTR/159-296\" :     \"KRGEDGLHGIVSCTACGQQVNHFQKDSIYRHPSLQVLICKNCFKYYMSDDISRDSDGMDEQCRWCAEGGNLICCDFCHNAFCKKCILRNLGRKELSTIMDENNQWYCYICHPEPLLDLVTACNSVFENLEQLLQQNKK\",\n",
        "\"ATRX_PONPY/159-296\" :     \"KRGEDGLHGIVSCTACGQQVNHFQKDSIYRHPSLQVLICKNCFKYYMSDDISRDSDGMDEQCRWCAEGGNLICCDFCHNAFCKKCILRNLGRKELSTIMDENNQWYCYICHPEPLLDLVTACNSVFDNLEQLLQQNKK\",\n",
        "\"ATRY_MACEU/1-96\"    :     \"VICTACGQQVNQFQKDSIYRHPTLNVLICKRWCAEGGNLICCDSCHNAFCKKCIWRNLGRKEISKIMNEKNEWHCYICCPEPLLDLIAVCDSVLEN\",\n",
        "\"CHR20_ARATH/472-601\":     \"RDDSQNPANNFRCTACNKVAVEVHSHPLLEVIVCMDCKRSIEDRVSKVDDSLERHCEWCGHIADLIDCRTCEKLFCASCIKRNIGEEYMSEAQSSGWDCCCCSPIPLQRLTLELEKAMRDKKSIELSSDS\",\n",
        "\"DNM3A_CHICK/447-579\":     \"EVRQKCRNIEDICISCGSLNVTLEHPLFIGGMCQNCKNCFLECAYQYDDDGYQSYCTICCGGREVLMCGNNNCCRCFCVECVDLLVGPGAAQAAIKEDPWNCYMCGHKGVYGLLRRREDWPSRLQMFFANNHD\",\n",
        "\"DNM3A_HUMAN/482-614\":     \"EVRQKCRNIEDICISCGSLNVTLEHPLFVGGMCQNCKNCFLECAYQYDDDGYQSYCTICCGGREVLMCGNNNCCRCFCVECVDLLVGPGAAQAAIKEDPWNCYMCGHKGTYGLLRRREDWPSRLQMFFANNHD\",\n",
        "\"DNM3A_MOUSE/478-610\":     \"EVRQKCRNIEDICISCGSLNVTLEHPLFIGGMCQNCKNCFLECAYQYDDDGYQSYCTICCGGREVLMCGNNNCCRCFCVECVDLLVGPGAAQAAIKEDPWNCYMCGHKGTYGLLRRREDWPSRLQMFFANNHD\",\n",
        "\"DNM3A_RAT/478-610\"  :     \"EVRQKCRNIEDICISCGSLNVTLEHPLFIGGMCQNCKNCFLECAYQYDDDGYQSYCTVCCEGRELLLCSNTSCCRCFCVECLEVLVGTGTAAEAKLQEPWSCYMCLPQRCHGVLRRREDWPSRLQMFFANNHD\",\n",
        "\"DNM3B_HUMAN/423-555\":     \"DVANNKSSLEDGCLSCGRKNPVSFHPLFEGGLCQTCRDRFLELFYMYDDDGYQSYCTVCCEGRELLLCSNTSCCRCFCVECLEVLVGAGTAEDAKLQEPWSCYMCLPQRCHGVLRRRKDWNVRLQAFFTSDTG\",\n",
        "\"DNM3B_MOUSE/428-560\":     \"EVTNNKGNLEDRCLSCGKKNPVSFHPLFEGGLCQSCRDRFLELFYMYDEDGYQSYCTVCCEGRELLLCSNTSCCRCFCVECLEVLVGAGTAEDVKLQEPWSCYMCLPQRCHGVLRRRKDWNMRLQDFFTTDPD\",\n",
        "\"DNM3C_MOUSE/309-441\":     \"DVTNNKGNLEDHCLSCGRKDPVSFHPLFEGGLCQSCRDRFLELFYMYDEDGYQSYCSICCSGETLLICGNPDCTRCYCFECVDSLVGPGTSGKVHAMSNWVCYLCLPSSRSGLLQRRKDWNMRLQDFFTTDPD\",\n",
        "\"DNM3L_HUMAN/41-173\" :     \"EVKANQRNIEDICICCGSLQVHTQHPLFEGGICAPCKDKFLDALFLYDDDGYQSYCTICCSGGTLFICESPDCTRCYCFECVDILVGPGTSERINAMACWVCFLCLPFSRSGLLQRRRKWRSQLKAFYDRESE\",\n",
        "\"DNM3L_MOUSE/75-207\" :     \"EVKVNRRSIEDICLCCGTLQVYTRHPLFEGGLCAPCKDKFLESLFLYDDDGHQSYCTICCSGHTLFICESPDCTRCYCFECVDILVGPGTSERINAMACWVCFLCLPFSRSGLLQRRKRWRHQLKAFHDQEGA\",\n",
        "\"DNM3L_RAT/76-208\"   :     \"EVNVNQRNIEDICLCCGSLQVYAQHPLFEGGICAPCKDKFLETLFLYDEDGHQSYCTICCGGREVLMCGNNNCCRCFCVECVDLLVGPGAAQAAIKEDPWNCYMCGHKGTYGLLRRRKKWRHQLKAFHDREGA\",\n",
        "}\n",
        "\n",
        "\n",
        "\n",
        "Example_seqs = [msa.MySeq(val, \"PROTEIN\") for val in example_seqs.values()]\n",
        "\n",
        "sm = msa.SubstMatrix()\n",
        "# sm.create_submat(1, -1, \"ACGT\")\n",
        "sm.read_submat_file(\"/content/bioinformatics_stockholm/msa/blosum62.mat\", \"\\t\")\n",
        "aseq = msa.PairwiseAlignment(sm, -1)\n",
        "ma = msa.MultipleAlignment(Example_seqs, aseq)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SHx54EqI-btV"
      },
      "source": [
        "The following code will take the alignment object produced in the last cell and save it to the server as a FASTA file. Note the filename \"multiple_alignment_example.fsa\" as the destination for the data. Each fasta entry is a sequence \"record\" that gets made here by supplying a sequence and an id.\n",
        "\n",
        "The FASTA format is simply a greater-than sign followed by identifying information, and then sequence information below that (for example):\n",
        "\n",
        "$> SEQUENCE 1$\n",
        "$MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG$\n",
        "$LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK$\n",
        "$IPQFASRKQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL$\n",
        "$MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL$\n",
        "\n",
        "And the file output will basically be a text file with entries that look like this. It is the same format we used for our NGS data in the previous exercises."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lIS5EryOVP7Y"
      },
      "source": [
        "from Bio import SeqIO\n",
        "from Bio.SeqRecord import SeqRecord\n",
        "from Bio.Seq import Seq\n",
        "import numpy as np\n",
        "i = 0\n",
        "with open(\"/content/bioinformatics_stockholm/multiple_alignment_example.fsa\", 'w') as filetowrite:\n",
        "    filetowrite.write('')\n",
        "for seq in ma.align_consensus():\n",
        "  # print(seq)\n",
        "  sorted_rec = SeqRecord(Seq(seq),\n",
        "                          id=str(list(example_seqs.keys())[i]),\n",
        "                          name=str(list(example_seqs.keys())[i]),\n",
        "                          description=f\"sequence: {i}\")\n",
        "  print(sorted_rec.seq)\n",
        "  with open(\"/content/bioinformatics_stockholm/multiple_alignment_example.fsa\", \"a\", ) as destination:\n",
        "    SeqIO.write(sorted_rec, destination, \"fasta\")  \n",
        "  i += 1\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gaarr6n8Yxz5"
      },
      "source": [
        "The code below will take the FASTA file produced in the previous step and produce an interactive plot that can be used to visualize the multiple alignment. You may click and drag the window to scroll through the positions."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eeqcFgIyQTaQ"
      },
      "source": [
        "import panel as pn\n",
        "pn.extension()\n",
        "aln = AlignIO.read(\"/content/bioinformatics_stockholm/multiple_alignment_example.fsa\", \"fasta\")\n",
        "p = msa.view_alignment(aln, plot_width=900)\n",
        "pn.pane.Bokeh(p)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "txWXdXA_AReG"
      },
      "source": [
        "## Problem 3.2\n",
        "Cas9 proteins are part of the CRISPR bacterial defense system against bacteria-infecting viruses, and Cas9 proteins are common throughout the bacterial kingdom. The proteins vary widely in size and sequence, however all known Cas9 enzymes contain a domain that cleaves the DNA strand complementary to the guide RNA. The following sequences were published in 2014 by Jennifer Doudna and colleagues, and their multiple sequence analysis of this collection led to their discovery of the protein's function, a discovery which was awarded the Nobel prize in 2020. **Use the scripts demonstrated above to perform a multiple sequence alignment on these sequences. Which amino acid (single letter code) is most represented in column/position 99 in the resulting multiple sequence alignment you get?** Note, you should not change the order of the sequences from how they are represented now. We will explore why this is the case later - but essentially the result may change if you alter the order. \n",
        "\n",
        "\"CAS9_ACTNH/513-675\": \"GTTIGYHTCQLDHIVPQAGPGSNNRRGNLVAVCERCNRSKSNTPFAVWAQKCGIPGTTIGYHTCQLDHIVPQAGPGSNNRRGNLVAVCERCNRSKSNTPFAVWAQKCGIPHVGVKEAIGRVRGWRKQTPNTSSEDLTRLKKEVIARLRRTQEDPEIDERSM\",\n",
        "\"CAS9A_STRTD/516-684\": \"CLYTGKTISIHDLINNSNQFEVDHILPLSITFDDSLANKVLVYATANQEKGQRTPYQALDCLYTGKTISIHDLINNSNQFEVDHILPLSITFDDSLANKVLVYATANQEKGQRTPYQALDSMDDAWSFRELKAFVRESKTLSNKKKEYLLTEEDISKFDVRKKFIERNL\",\n",
        "\"CAS9B_STRTD/771-928\": \"LQNGKDMYTGDDLDIDRLSNYDIDHIIPQAFLKDNSIDNKVLVSSASNRGKSDDVPSLLQNGKDMYTGDDLDIDRLSNYDIDHIIPQAFLKDNSIDNKVLVSSASNRGKSDDVPSLEVVKKRKTFWYQLLKSKLISQRKFDNLTKAERGGLSPEDKAGFIQRQL\",\n",
        "\"CAS9_CAMJE/487-637\": \"KIKISDLQDEKMLEIDHIYPYSRSFDDSYMNKVLVFTKQNQEKLNQTPFEAFGNDSAKWKIKISDLQDEKMLEIDHIYPYSRSFDDSYMNKVLVFTKQNQEKLNQTPFEAFGNDSAKWQKIEVLAKNLPTKKQKRILDKNYKDKEQKNFKDRNL\",\n",
        "\"CAS9_CORDI/504-665\": \"GSPITFSNSEMDHIVPRAGQGSTNTRENLVAVCHRCNQSKGNTPFAIWAKNTSIEGSPITFSNSEMDHIVPRAGQGSTNTRENLVAVCHRCNQSKGNTPFAIWAKNTSIEGVSVKEAVERTRHWVTDTGMRSTDFKKFTKAVVERFQRATMDEEIDARSM\",\n",
        "\"CAS9_FRATN/897-1046\": \"NLTDGDFDGAKEELDHIIPRSHKKYGTLNDEANLICVTRGDNKNKGNRIFCLRDLADNYNLTDGDFDGAKEELDHIIPRSHKKYGTLNDEANLICVTRGDNKNKGNRIFCLRDLADNYKLKQFETTDDLEIEKKIADTIWDANKKDFKFGNY\",\n",
        "\"CAS9_LISIN/773-924\": \"MYTGQDLDIHNLSNYDIDHIVPQSFITDNSIDNLVLTSSAGNREKGDDVPPLEIVRKRMYTGQDLDIHNLSNYDIDHIVPQSFITDNSIDNLVLTSSAGNREKGDDVPPLEIVRKRKVFWEKLYQGNLMSKRKFDYLTKAERGGLTEADKARFIHRQL\",\n",
        "\"CAS9_NEIM8/512-667\": \"SGKEINLGRLNEKGYVEIDHALPFSRTWDDSFNNKVLVLGSENQNKGNQTPYEYFNGKDNSGKEINLGRLNEKGYVEIDHALPFSRTWDDSFNNKVLVLGSENQNKGNQTPYEYFNGKDNSREWQEFKARVETSRFPRSKKQRILLQKFDEDGFKERNL\",\n",
        "\"CAS9_NEIMA/512-667\": \"SGKEINLGRLNEKGYVEIDHALPFSRTWDDSFNNKVLVLGSENQNKGNQTPYEYFNGKDNSGKEINLGRLNEKGYVEIDHALPFSRTWDDSFNNKVLVLGSENQNKGNQTPYEYFNGKDNSREWQEFKARVETSRFPRSKKQRILLQKFDEDGFKERNL\",\n",
        "\"CAS9_PASMU/504-659\": \"SGKEINIHRLNEKGYVEIDHALPFSRTWDDSFNNKVLVLASENQNKGNQTPYEWLQGKINSGKEINIHRLNEKGYVEIDHALPFSRTWDDSFNNKVLVLASENQNKGNQTPYEWLQGKINSERWKNFVALVLGSQCSAAKKQRLLTQVIDDNKFIDRNL\",\n",
        "\"CAS9_STAAU/480-646\": \"LEAIPLEDLLNNPFNYEVDHIIPRSVSFDNSFNNKVLVKQEENSKKGNRTPFQYLSSSDSLEAIPLEDLLNNPFNYEVDHIIPRSVSFDNSFNNKVLVKQEENSKKGNRTPFQYLSSSDSKISYETFKKHILNLAKGKGRISKTKKEYLLEERDINRFSVQKDFINRNL\",\n",
        "\"CAS9_STRMU/770-921\": \"MYTGEELDIDYLSQYDIDHIIPQAFIKDNSIDNRVLTSSKENRGKSDDVPSKDVVRKMMYTGEELDIDYLSQYDIDHIIPQAFIKDNSIDNRVLTSSKENRGKSDDVPSKDVVRKMKSYWSKLLSAKLITQRKFDNLTKAERGGLTDDDKAGFIKRQL\",\n",
        "\"CAS9_STRP1/770-921\": \"MYVDQELDINRLSDYDVDHIVPQSFLKDDSIDNKVLTRSDKNRGKSDNVPSEEVVKKMMYVDQELDINRLSDYDVDHIVPQSFLKDDSIDNKVLTRSDKNRGKSDNVPSEEVVKKMKNYWRQLLNAKLITQRKFDNLTKAERGGLSELDKAGFIKRQL\",\n",
        "\"CAS9_STRTR/792-949\": \"LQNGKDMYTGDDLDIDRLSNYDIDHIIPQAFLKDNSIDNKVLVSSASNRGKSDDFPSLLQNGKDMYTGDDLDIDRLSNYDIDHIIPQAFLKDNSIDNKVLVSSASNRGKSDDFPSLEVVKKRKTFWYQLLKSKLISQRKFDNLTKAERGGLLPEDKAGFIQRQL\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d3QsvTPsZOPI"
      },
      "source": [
        "## Problem 3.3\n",
        "The following set of sequences are from the Nsp3e nucleic acid-binding domain of betacoronaviruses - the family includes members like SARS, SARS-CoV-2, BatCoV RaTG13 and Bat-SARS-like coronavirus (BATSL-CoVZXC21 and BAT-SL-CoVZC45). \n",
        "\n",
        "**Generate a \"dumb\" consensus for the Nsp3e sequences.**  First run the code below, and it will create a list the_sequences = ma.align_consensus() which is a list of each of the aligned (gaps inserted) sequences in the multiple sequence alignment. Your task is to figure out how to determine which base is the most prevalent in each column, and print out a final consensus sequence. **For canvas: what is the consensus sequence? Include any gaps in this answer.** To help you out, we have provided a small function called most_frequent() which will check which is the most-prevalent character of a *list* or *array* of characters. For example if you give it [\"A\", \"A\", \"C\"] it will return the string \"A\". \n",
        "\n",
        "(note again - do not change the order of sequences from how they are currently listed as this will affect the result of the alignment - we will discuss why later). \n",
        "\n",
        "\n",
        "Additional guidance: Wondering where to start? This problem is a matter of knowing how to index/slice arrays, lists, and strings in an intermediate/advanced way. The goal is essentially to fetch the columns of this multiple sequence alignment, however the problem is that the sequences are provided as a list of linear strings. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eU5ENah4zj3J"
      },
      "source": [
        "def most_frequent(List):\n",
        "    List = list(List)\n",
        "    return max(set(List), key = List.count) "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DSMdWA2bZcjF"
      },
      "source": [
        "from Bio import AlignIO\n",
        "from Bio.Align import MultipleSeqAlignment\n",
        "\n",
        "Nsp3e_seqs = {\n",
        "\"R1AB_BC133/1885-2002\":\"TIKYSPATILPGSVYSNSCLVGVDGTPGSDTISKFFNDLLGFDETKPISKKLTYSLLPNEDGDVLLSEFNNYNPVYKKGVMLKGKPILWVNNGVCDSALNKPNRASLRQLYDVAPIVL\",\n",
        "\"R1AB_BC279/1894-2004\":\"PIDLVPTQPLPNASFDNFKLTCSNTKFADDLNQMTGFKKPASRELSVTFFPDLNGDVVAIDYRHYSASFKKGAKLLHKPIIWHINQTTNKTTYKPNTWCLRCLWSTKPVET\",\n",
        "\"R1AB_BCHK3/1882-1992\":\"PIDLVPTQPMPNASFDNFKLTCSNTKFADDLNQMTGFKKPASRELTVTFFPDLNGDVVAIDYRHYSTSFKKGAKLVHKPILWHINQTTNKTTYKPNIWCLRCLWSTKPVDT\",\n",
        "\"R1AB_BCHK4/1878-1995\":\"TIKYSPATILPGSVYSNSCLVGVDGTPGSDTISKFFNDLLGFDETKPISKKLTYSLLPNEDGDVLLSEFSNYNPVYKKGVMLKGKPILWVNNGVCDSALNKPNRASLRQLYDVAPIVL\",\n",
        "\"R1AB_BCHK5/1916-2033\":\"VIEYSPATILSGSVYTNSCLVGHDGTIGSDAISSSFNNLLGFDNSKPVSKKLTYSFFPDFEGDVILTEYSTYDPIYKNGAMLHGKPILWVNNSKFDSALNKFNRATLRQVYDIAPVTL\",\n",
        "\"R1AB_BCHK9/1770-1870\":\"PIEVVAAPKLVTSYDGFYLSSCQNPQLAESFNKAINATKTGPMKLLTMYPNVAGDVVAISDDNVVAHPYGSLHMGKPVLFVTRPNTWKKLVPLLSTVVVNT\",\n",
        "\"R1AB_BCRP3/1886-1996\":\"PIDLIPTQPLPNASFDNFKLTCSNTKFADDLNQMTGFTKPASRELSVTFFPDLNGDVVAIDYRHYSASFKKGAKLLHKPIVWHINQATTKTTFKPNTWCLRCLWSTKPVDT\",\n",
        "\"R1A_BC133/1885-2002\":\"TIKYSPATILPGSVYSNSCLVGVDGTPGSDTISKFFNDLLGFDETKPISKKLTYSLLPNEDGDVLLSEFNNYNPVYKKGVMLKGKPILWVNNGVCDSALNKPNRASLRQLYDVAPIVL\",\n",
        "\"R1A_BC279/1894-2004\":\"PIDLVPTQPLPNASFDNFKLTCSNTKFADDLNQMTGFKKPASRELSVTFFPDLNGDVVAIDYRHYSASFKKGAKLLHKPIIWHINQTTNKTTYKPNTWCLRCLWSTKPVET\",\n",
        "\"R1A_BCHK3/1882-1992\":\"PIDLVPTQPMPNASFDNFKLTCSNTKFADDLNQMTGFKKPASRELTVTFFPDLNGDVVAIDYRHYSTSFKKGAKLVHKPILWHINQTTNKTTYKPNIWCLRCLWSTKPVDT\",\n",
        "\"R1A_BCHK4/1878-1995\":\"TIKYSPATILPGSVYSNSCLVGVDGTPGSDTISKFFNDLLGFDETKPISKKLTYSLLPNEDGDVLLSEFSNYNPVYKKGVMLKGKPILWVNNGVCDSALNKPNRASLRQLYDVAPIVL\",\n",
        "\"R1A_BCHK5/1916-2033\":\"VIEYSPATILSGSVYTNSCLVGHDGTIGSDAISSSFNNLLGFDNSKPVSKKLTYSFFPDFEGDVILTEYSTYDPIYKNGAMLHGKPILWVNNSKFDSALNKFNRATLRQVYDIAPVTL\",\n",
        "\"R1A_BCHK9/1770-1870\":\"PIEVVAAPKLVTSYDGFYLSSCQNPQLAESFNKAINATKTGPMKLLTMYPNVAGDVVAISDDNVVAHPYGSLHMGKPVLFVTRPNTWKKLVPLLSTVVVNT\",\n",
        "\"R1A_BCRP3/1886-1996\":\"PIDLIPTQPLPNASFDNFKLTCSNTKFADDLNQMTGFTKPASRELSVTFFPDLNGDVVAIDYRHYSASFKKGAKLLHKPIVWHINQATTKTTFKPNTWCLRCLWSTKPVDT\",\n",
        "\"R1AB_CVBEN/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTICDILNAKLGFDSSKEFVEYKVTEWPTATGDVVLATDDLYVKRYERGCITFGKPVIWLSHEQASLNSLTYFNRPLLVD\",\n",
        "\"R1AB_CVBLU/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTICDILNAKLGFDSSKEFVEYKVTEWPTATGDVVLATDDLYVKRYERGCITFGKPVIWLSHEQASLNSLTYFNRPLLVD\",\n",
        "\"R1AB_CVBM/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTVCDILNAKLGFDSSKEFVEYKVTEWPTATGDVVLATDDLYVKRYERGCITFGKPVIWLSHEQASLNSLTYFNRPLLVD\",\n",
        "\"R1AB_CVBQ/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTVCDILNAKLGFDSSKEFVEYKVTEWPTATGDVVLATDDLYVKRYERGCITFGKPVIWLSHEQASLNSLTYFNRPLLVD\",\n",
        "\"R1AB_CVHN1/1992-2093\":\"IKAQFKPFAKVDGVYTNFKLVGHDICAQLNDKLGFNVDLPFVEYKVTVWPVATGDVVLASDDLYVKRYFKGCETFGKPVIWFCHDEASLNSLTYFNKPSFKS\",\n",
        "\"R1AB_CVHN2/1962-2063\":\"IKAQFKPFAKVDGVYTNFKLVGHDICAQLNDKLGFNVDLPFVEYKVTVWPVATGDVVLASDDLYVKRYFKGCETFGKPVIWLCHDEASLNSLTYFNKPSFKS\",\n",
        "\"R1AB_CVHN5/1942-2043\":\"IKAQFKPFAKVDGVYTNFKLVGHDICAQLNDKLGFNVDLPFVEYKVTVWPVATGDVVLASDDLYVKRYFKGCETFGKPVIWFCHDEASLNSLTYFNKPSFKS\",\n",
        "\"R1AB_CVHOC/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTVCDSLNSKLGFDSSKEFVEYKITEWPTATGDVVLANDDLYVKRYERGCITFGKPVIWLSHEKASLNSLTYFNRPLLVD\",\n",
        "\"R1AB_CVM2/1898-1999\":\"IKAQFRTFEKVEGVYTNFKLVGHSIAEKFNAKLGFDCNSPFTEYKITEWPTATGDVVLASDDLYVSRYSGGCVTFGKPVIWLGHEEASLKSLTYFNRPSVVC\",\n",
        "\"R1AB_CVMA5/1951-2052\":\"IKAQFRTFEKVDGVYTNFKLVGHSIAEKLNAKLGFDCNSPFVEYKITEWPTATGDVVLASDDLYVSRYSSGCITFGKPVVWLGHEEASLKSLTYFNRPSVVC\",\n",
        "\"R1AB_CVMJH/1950-2051\":\"IKAQFRTFEKVEGVYTNFKLVGHDIAEKLNAKLGFDCNSPFMEYKITEWPTATGDVVLASDDLYVSRYSGGCVTFGKPVIWRGHEEASLKSLTYFNRPSVVC\",\n",
        "\"R1AB_MERS1/1837-1954\":\"PVTYSPATILAGSVYTNSCLVSSDGQPGGDAISLSFNNLLGFDSSKPVTKKYTYSFLPKEDGDVLLAEFDTYDPIYKNGAMYKGKPILWVNKASYDTNLNKFNRASLRQIFDVAPIEL\",\n",
        "\"R1AB_SARS/1888-1998\":\"PIDLVPTQPLPNASFDNFKLTCSNTKFADDLNQMTGFTKPASRELSVTFFPDLNGDVVAIDYRHYSASFKKGAKLLHKPIVWHINQATTKTTFKPNTWCLRCLWSTKPVDT\",\n",
        "\"R1AB_SARS2/1911-2021\":\"PIDLVPNQPYPNASFDNFKFVCDNIKFADDLNQLTGYKKPASRELKVTFFPDLNGDVVAIDYKHYTPSFKKGAKLLHKPIVWHVNNATNKATYKPNTWCIRCLWSTKPVET\",\n",
        "\"R1A_CVBEN/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTICDILNAKLGFDSSKEFVEYKVTEWPTATGDVVLATDDLYVKRYERGCITFGKPVIWLSHEQASLNSLTYFNRPLLVD\",\n",
        "\"R1A_CVBLU/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTICDILNAKLGFDSSKEFVEYKVTEWPTATGDVVLATDDLYVKRYERGCITFGKPVIWLSHEQASLNSLTYFNRPLLVD\",\n",
        "\"R1A_CVBM/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTVCDILNAKLGFDSSKEFVEYKVTEWPTATGDVVLATDDLYVKRYERGCITFGKPVIWLSHEQASLNSLTYFNRPLLVD\",\n",
        "\"R1A_CVBQ/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTVCDILNAKLGFDSSKEFVEYKVTEWPTATGDVVLATDDLYVKRYERGCITFGKPVIWLSHEQASLNSLTYFNRPLLVD\",\n",
        "\"R1A_CVHN1/1992-2093\":\"IKAQFKPFAKVDGVYTNFKLVGHDICAQLNDKLGFNVDLPFVEYKVTVWPVATGDVVLASDDLYVKRYFKGCETFGKPVIWFCHDEASLNSLTYFNKPSFKS\",\n",
        "\"R1A_CVHN2/1962-2063\":\"IKAQFKPFAKVDGVYTNFKLVGHDICAQLNDKLGFNVDLPFVEYKVTVWPVATGDVVLASDDLYVKRYFKGCETFGKPVIWLCHDEASLNSLTYFNKPSFKS\",\n",
        "\"R1A_CVHN5/1942-2043\":\"IKAQFKPFAKVDGVYTNFKLVGHDICAQLNDKLGFNVDLPFVEYKVTVWPVATGDVVLASDDLYVKRYFKGCETFGKPVIWFCHDEASLNSLTYFNKPSFKS\",\n",
        "\"R1A_CVHOC/1906-2007\":\"IKAQFKTFEKVDGVYTNFKLIGHTVCDSLNSKLGFDSSKEFVEYKITEWPTATGDVVLANDDLYVKRYERGCITFGKPVIWLSHEKASLNSLTYFNRPLLVD\",\n",
        "\"R1A_CVM2/1898-1999\":\"IKAQFRTFEKVEGVYTNFKLVGHSIAEKFNAKLGFDCNSPFTEYKITEWPTATGDVVLASDDLYVSRYSGGCVTFGKPVIWLGHEEASLKSLTYFNRPSVVC\",\n",
        "\"R1A_CVMA5/1951-2052\":\"IKAQFRTFEKVDGVYTNFKLVGHSIAEKLNAKLGFDCNSPFVEYKITEWPTATGDVVLASDDLYVSRYSSGCITFGKPVVWLGHEEASLKSLTYFNRPSVVC\",\n",
        "\"R1A_CVMJH/1950-2051\":\"IKAQFRTFEKVEGVYTNFKLVGHDIAEKLNAKLGFDCNSPFMEYKITEWPTATGDVVLASDDLYVSRYSGGCVTFGKPVIWRGHEEASLKSLTYFNRPSVVC\",\n",
        "\"R1A_MERS1/1837-1954\":\"PVTYSPATILAGSVYTNSCLVSSDGQPGGDAISLSFNNLLGFDSSKPVTKKYTYSFLPKEDGDVLLAEFDTYDPIYKNGAMYKGKPILWVNKASYDTNLNKFNRASLRQIFDVAPIEL\",\n",
        "\"R1A_SARS/1888-1998\":\"PIDLVPTQPLPNASFDNFKLTCSNTKFADDLNQMTGFTKPASRELSVTFFPDLNGDVVAIDYRHYSASFKKGAKLLHKPIVWHINQATTKTTFKPNTWCLRCLWSTKPVDT\",\n",
        "\"R1A_SARS2/1911-2021\":\"PIDLVPNQPYPNASFDNFKFVCDNIKFADDLNQLTGYKKPASRELKVTFFPDLNGDVVAIDYKHYTPSFKKGAKLLHKPIVWHVNNATNKATYKPNTWCIRCLWSTKPVET\"}\n",
        "\n",
        "Nsp3e_Seqs = [msa.MySeq(val, \"PROTEIN\") for val in Nsp3e_seqs.values()]\n",
        "\n",
        "sm = msa.SubstMatrix()\n",
        "# sm.create_submat(1, -1, \"ACGT\")\n",
        "sm.read_submat_file(\"/content/bioinformatics_stockholm/msa/blosum62.mat\", \"\\t\")\n",
        "aseq = msa.PairwiseAlignment(sm, -1)\n",
        "ma_Nsp3e = msa.MultipleAlignment(Example_seqs, aseq)\n",
        "the_sequences = ma_Nsp3e.align_consensus()\n",
        "print(the_sequences)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mbgCEI7c3TDX"
      },
      "source": [
        "## optional challenge problem\n",
        "**Generate an amino acid frequency table like those discussed in the lecture. Use the multiple alignment data from the zinc finger example earlier in this section. Your matrix should have a row for every amino acid type and a column corresponding to each column in the multiple alignment. And elements of the matrix should have values corresponding to the amino acid frequency in that column."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tKFnhSYk48CD"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}