2_pairwise_sequence_comparison.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": "2_pairwise_sequence_comparison.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/2_pairwise_sequence_comparison.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yOPSrMNgZIyO"
      },
      "source": [
        "# Dot plots \n",
        "In today's exercise we will use a homemade module to produce dot plots and visualize dynamic programming alignments.\n",
        "\n",
        "Run the code below to download biopython on the server and install our dot_plotting module as well as the dynamic_programming module."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "thgytQXvpjeC"
      },
      "source": [
        ""
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "91UtzNQS4Dn_"
      },
      "source": [
        "pip install biopython"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "15VmTe7Epjtm"
      },
      "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(\"dot_plotting\", \"/content/bioinformatics_stockholm/dot_plotting/dot_plotting.py\")\n",
        "dot_plotting = importlib.util.module_from_spec(spec)\n",
        "spec.loader.exec_module(dot_plotting)\n",
        "\n",
        "spec = importlib.util.spec_from_file_location(\"dynamic_programming\", \"/content/bioinformatics_stockholm/dynamic_programming/dynamic_programming.py\")\n",
        "dynamic_programming = importlib.util.module_from_spec(spec)\n",
        "spec.loader.exec_module(dynamic_programming)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ldxjHpTf5uEv"
      },
      "source": [
        "## dot plotting related genes from NCBI databases\n",
        "\n",
        "Recall that we can use the following block of code to fetch a full record from NCBI's nucleotide database. And that we can use the code below that to fetch sequences from specific records if we supply their accession numbers/IDs. Note here how we are translating the sequences to protein after fetching them from the nucleotide records."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4sGCsW-dU0Yr"
      },
      "source": [
        "from Bio import Entrez\n",
        "Entrez.email = \"ithof@kth.se\"  # Always tell NCBI who you are\n",
        "handle = Entrez.efetch(db=\"nucleotide\", id=\"AJ457860.1\", rettype=\"gb\", retmode=\"text\")\n",
        "print(handle.read())\n",
        "handle.close()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lFd44u0o7EhI"
      },
      "source": [
        "from Bio import SeqIO\n",
        "from Bio import Entrez\n",
        "Entrez.email = \"ithof@kth.se\" \n",
        "handle = Entrez.efetch(db=\"nucleotide\", id=\"AJ457860.1\", rettype=\"gb\", retmode=\"text\")\n",
        "mouse_albumin = SeqIO.read(handle, \"genbank\").seq.translate()  # translate!\n",
        "handle.close()\n",
        "handle = Entrez.efetch(db=\"nucleotide\", id=\"X74045.1\", rettype=\"gb\", retmode=\"text\")\n",
        "hourse_albumin = SeqIO.read(handle, \"genbank\").seq.translate() # translate!\n",
        "handle.close()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rfmMEiO5VqhW"
      },
      "source": [
        "Try running our dot-plotting function make_quick_plot() which takes as its input the size of a sliding window and two sequences."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N477GUtx8fuQ"
      },
      "source": [
        "dot_plotting.make_quick_plot(window = 3, sequence_1 = mouse_albumin, sequence_2 = hourse_albumin) # this is our homemade dotplot function - it takes as parameters the window size, and the two sequences"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nInhVi3lV_jJ"
      },
      "source": [
        "## Problem 2.1\n",
        "Go through the following accession numbers for different versions (i.e. derived from different species) of the serum albumin protein and, using dot plots, visually inspect them for their similarity with the human one: AJ457860.1. You may find it convenient to make a function for importing the sequences and generating the dotplots, as this could reduce clutter. What is the latin name (with the form: Genus species) of the animal with serum albumin most similar to that of the human? Do you find the result surprising?   \n",
        "\n",
        "M90463.1 \n",
        "\n",
        "DQ145791.1  \n",
        "\n",
        "XM_014546593.1 \n",
        "\n",
        "NM_134326.2 "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tjWZwGNStwQv"
      },
      "source": [
        "# Substitution scoring of aligned sequences\n",
        "\n",
        "The following python dictionary will be used to store the BLOSUM62 matrix. We can use it just like other dictionaries by supplying the correct key, which will return the corresponding substitution score. The key must have the following format (for example for amino acids W and F):\n",
        "\n",
        "(\"W\",\"F\")\n",
        "\n",
        "Try accessing a few scores this way for different pairs. Remember that to access the value of a dictionary, we supply the key in brackets: dict[key]\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IE1A0eFPyEnX"
      },
      "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": "code",
      "metadata": {
        "id": "BBe6Qy1CynAM"
      },
      "source": [
        "example_1 = Blosum62[(\"W\",\"F\")]\n",
        "example_2 = Blosum62[('N','B')]\n",
        "print(example_1, example_2)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "POINRp7by_BH"
      },
      "source": [
        "## problem 2.2\n",
        "**Make a function that calculates the scores of following two alignments of these sequences using our Blosum62 matrix. Treat all gaps with a fixed penalty of -5. The sequences are globins from two different species of worm.**\n",
        "\n",
        "**Alignment 1:** \n",
        "\n",
        "*Caenorhabdatis elegans:* \n",
        "\n",
        "\"--------------MSMNRQEISDLCVKSLEGRMVGTEAQNIENGNAFYRYFFTNFPDLRVYFKGAEKYTADDVKKSERFDKQGQRILLACHLLANVYTNEEVFKGYVRETINRHRIYKMDPALWMAFF-TVFTGYLESVGC-LNDQQKAAWMALGKEFNAESQTHLKNSNLPHV\"\n",
        "\n",
        "*Nippostrongylus brasiliensis:* \n",
        "\n",
        "\"MLWFVAVCFAIASVSAMSPADVKKHTVESMKAVPVGRD--KAQNGIDFYKFFFTHHKDLRKFFKGAENFGADDVQKSKRFEKQGTALLLAVHVLANVYDNQAVFHGFVRELMNRHEKRGVDPKLWKIFFDDVWVPFLESKGAKLSGDAKAAWKELNKNFNSEAQHQLEKLGLPHA\"\n",
        "\n",
        "\n",
        "**Alignment 2:** \n",
        "\n",
        "*Caenorhabdatis elegans:* \n",
        "\n",
        "\"-------------MS-MNRQEIS--DLCVK--SLEGRM----VGTE-AQNIENG-NAFYRYFFTNFPDLRVYFKGAEKYTADDVKKSERFDKQGQRILLACHLLANVYTNEEVFKGYVRETINRH--R-IYKMDPALW-MAFF--T-V-FTGYLESVGC-LN-DQQKAAWMALGKEFNAESQTH-L-KNS-NLPHV\"\n",
        "\n",
        "*Nippostrongylus brasiliensis:* \n",
        "\n",
        "\"MLWFVAVCFAIASVSAM-----SPAD--VKKHTVES-MKAVPVGRDKAQ---NGID-FYKFFFTHHKDLRKFFKGAENFGADDVQKSKRFEKQGTALLLAVHVLANVYDNQAVFHGFVRELMNRHEKRGV---DPKLWKI-FFDDVWVPF---LESKGAKLSGD-AKAAWKELNKNFNSEAQ-HQLEK--LGLPHA\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true,
        "id": "wh79XGjHZIyQ"
      },
      "source": [
        "# Dynamic programming. \n",
        "## visualizing alignment algorithms\n",
        "Try the examples below which run a homemade Needleman-Wunsch visualization, and below that a corresponding Smith-Waterman visualization. Note the dependence of the result on gap penalty value - try changing it (the last argument supplied to the function). Pretty neat right?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0WogJ9AzZIyQ"
      },
      "source": [
        "# import dynamic_programming.dynamic_programming as dp # another homemade module - download and unzip the folder and add it to your working directory\n",
        "sequence_1 = 'THISLINE' # this is an ordinary string (a piece of text), but we have been careful to use characters that are also blosum62-compatible amino acid residue codes\n",
        "sequence_2 = 'ISALIGNED'\n",
        "viz_align = dynamic_programming.needlemanwunsch(sequence_1,sequence_2,fulloutput=False, linear_gap_penalty=-4) # this will run the vizualization of needleman wunsch, it takes as its inputs an option for outputting more data "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Fovv-XH_ZIyR"
      },
      "source": [
        "viz_align = dynamic_programming.smithwaterman(sequence_1,sequence_2,fulloutput=False, linear_gap_penalty=-3)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6zHdr0ggJ6kQ"
      },
      "source": [
        "## generating pairwise alignments with biopython\n",
        "\n",
        "Visualizing alignments is good for understanding how they work, but often we are more interested in simply getting the alignments and lots of them. \n",
        "\n",
        "Let's try out biopython's pairwise alignment functions. The functions are a bit complex, with many adjustable parameters, and the outputs (alignments) can be formatted in multiple ways. Thus we will try to show here a bit how it works. The first step is to import the proper modules, Seq for handling sequences and pairwise2 which has our alignment function. Then we will define a score matrix using a dictionary for DNA."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DBLMV67-PDh4",
        "outputId": "3716ccee-7f99-461d-8b7c-2c2f5164b60f"
      },
      "source": [
        "from Bio.Seq import Seq\n",
        "from Bio import pairwise2\n",
        "score_matrix = {\n",
        "    ('A', 'A'): 4, ('A', 'T'): -4, ('A', 'C'): -4, ('A', 'G'): -4,\n",
        "    ('T', 'A'): -4, ('T', 'T'): 4, ('T', 'C'): -4, ('T', 'G'): -4,\n",
        "    ('C', 'A'): -4, ('C', 'T'): -4, ('C', 'C'): 4, ('C', 'G'): -4,\n",
        "    ('G', 'A'): -4, ('G', 'T'): -4, ('G', 'C'): -4, ('G', 'G'): 4 }\n",
        "\n",
        "score_matrix[('A', 'A')]  # example - if we want to fetch a particular substitution event"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "4"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "POQS4HRFPnwy"
      },
      "source": [
        "Next, we want to use pairwise2's local alignment function which utilizes the Smith-Waterman procedure we discussed in lecture. Let's check what the function needs using help()"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nxdvMYyrPctU"
      },
      "source": [
        "help(pairwise2.align.localds)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ejr8Zw9sPzRR"
      },
      "source": [
        "It's a lot of text, but the important thing we get from this is what the function needs as input parameters:\n",
        "\n",
        "**localds(sequenceA, sequenceB, match_dict, open, extend)**\n",
        "\n",
        "From this we understand that the first two inputs are sequences that will be aligned, the third parameter is our score matrix dictionary, and the last two are opening and extension gap penalties.\n",
        "\n",
        "Let's try this out then for two example sequences:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aHwLupJ2OPS-"
      },
      "source": [
        "target_sequence = Seq(\"GACGATCAGCATCGAGCATCGAATCGAGCGACTAGCGACGATCGACTACGAGCTACGAGCTACGAC\")\n",
        "query_sequence = Seq(\"ATCAGCGATAGCTCGA\")\n",
        "\n",
        "testaln = pairwise2.align.localds(query_sequence,target_sequence, match_dict = score_matrix, open = -3, extend = -1)[0] # the [0] is because aligments sometimes return multiple results and always spit out a list, we just want the first result \n",
        "print(testaln)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3CECefbRQSZ3"
      },
      "source": [
        "The function has succeeded, but the format of the alignment is a bit hard to read. It is in fact a list of multiple components: the first sequence with gaps, the second sequence with gaps, the score, the beginning of the aligned region, and the end of the aligned region. To access the different components we would use index rules, e.g. testaln[0] to fetch the first sequence.\n",
        "\n",
        "Let's fetch the score of the alignment for example:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gZkrhHthY4_R"
      },
      "source": [
        "the_score = testaln[2]\n",
        "print(the_score)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xkoFfHsIZImw"
      },
      "source": [
        "We can use biopython's formatting module to help us better view the alignment though. Again we check help() to find out what information the function needs to work."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6V_rVDapSC6a"
      },
      "source": [
        "help(pairwise2.format_alignment)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LJcdXeauTS4G"
      },
      "source": [
        "We see here that the input parameters for the format_alignment() function are the different elements of the alignment object that was spit out earlier. We plug them in then to get a nicer looking view:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KJNoNABiQaFy"
      },
      "source": [
        "print(pairwise2.format_alignment(align1=testaln[0], align2=testaln[1], score=testaln[2], begin=testaln[3], end=testaln[4], full_sequences=True))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pN-A1_vvN7Uv"
      },
      "source": [
        "## problem 2.3\n",
        "Now to practice a bit of brute force. We will use again our sample NGS dataset, an experimental result with multiple \"reads\" or individual short sequences comprising the dataset. Let’s search for the prevalence of a known reference sequence which we expect to find in our data “CTGATTTCTTTCCCTTCTTCCTCAGTTTGGATCAGATCGG”. \n",
        "\n",
        "**Use local alignment with biopython to systematically perform pairwise alignments between the reference sequence and all sequences in the NGS dataset using a gap opening penalty of -5 and an extension penalty of -1. How many alignments have a score greater than 100?**\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sLDbccboVFTx"
      },
      "source": [
        "### the code here will fetch the NGS fasta data and store all sequences into a list called all_sequences\n",
        "import sys\n",
        "from os import path\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",
        "\n",
        "all_sequences = [] \n",
        "for seq_record in SeqIO.parse(\"/content/bioinformatics_stockholm/ngs_data.fsa\", \"fasta\"): # this goes through each line of the file\n",
        "    all_sequences.append(seq_record.seq) #this adds the line from the file to our growing list"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "x3hmoUnsakKl"
      },
      "source": [
        "# optional challenge problem\n",
        "Say we have a set of error-free mini-reads from an NGS experiment:\n",
        "segs = [\"AAAAATGCTAA\", \"TGCTAAACTCT\", \"ACTCTTGTCCT\", \"GTCCTTACGACT\",\"CGACTTACGATAGC\", \"ATAGCGCATAGGGT\"]\n",
        "\n",
        "Each sequence overlaps partially with one or two of the other sequences. Write an algorithm that assembles a single sequence from these overlapping pieces. This is similar to problems in genome assembly."
      ]
    }
  ]
}