genome_assembly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7247fb34",
   "metadata": {},
   "source": [
    "# Gene structure\n",
    "\n",
    "A **gene** is a piece of DNA which is **transcribed into RNA** and then **translated into amino acids**. \n",
    "\n",
    "However, not the entire DNA sequence is translated into the final protein sequence, e.g. **introns are excluded** and **exons can be combined** in different ways to give different products known as **protein isoforms**.\n",
    "\n",
    "Below a figure of the **P53 gene** highlighting the gene structure as provided by the crresponding [transcript page](https://www.ensembl.org/Homo_sapiens/Transcript/Summary?g=ENSG00000141510;r=17:7668402-7687538;t=ENST00000269305) in the **Ensembl database**.\n",
    "\n",
    "\n",
    "![zscore](figures/ensemble_p53.png)\n",
    "\n",
    "\n",
    "- Exons --> boxes\n",
    "    - Coding --> filled boxes\n",
    "    - Non coding (UnTranslated Regions - UTR) --> white boxes\n",
    "- Introns --> lines\n",
    "- Antisense (reverse) direction. From 5' (right) --> to 3' (left)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18e2e31d",
   "metadata": {},
   "source": [
    "# Genome assembly\n",
    "\n",
    "In this Notebook you will learn how to map fragments to a known reference sequence. This is a common task in genome assembly and population studies. \n",
    "\n",
    "- Genome assembly --> identify overlapping regions of DNA fragments\n",
    "- Population studies --> map seqeunced fragments to a reference sequence and identify mutations / variations\n",
    "\n",
    "Specifically, you will learn how to use pairwise alignments to map exons sequences to a reference DNA sequence.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Glossary\n",
    "- DNA → Full gene DNA sequence (exons + introns) - 19,170 bps\n",
    "- cDNA → Complementary DNA / Transcript (exons, red + white blocks) - 2,512 bps\n",
    "    - The same gene can have multiple transcripts (different combinations of exons)\n",
    "    - Only one transcript is the canonical transcript, the others are called isoforms\n",
    "- CDS → Coding DNA sequence (exons, red blocks) - 1,182 bps\n",
    "- Protein → Translated sequence - 393 residues (1/3 CDS)\n",
    "\n",
    "\n",
    "### Dataset\n",
    "The dataset fasta file (ensemble_p53.fasta) has been created from the Ensembl transcript 1 page (see above) of P53 and clicking on:\n",
    "\n",
    ">\"Export data\" -> \"Genomics None\" -> \"Next\" -> \"Text\"\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1d745bd0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/damiano/miniconda3/envs/biological_data/lib/python3.9/site-packages/Bio/SubsMat/__init__.py:126: BiopythonDeprecationWarning: Bio.SubsMat has been deprecated, and we intend to remove it in a future release of Biopython. As an alternative, please consider using Bio.Align.substitution_matrices as a replacement, and contact the Biopython developers if you still need the Bio.SubsMat module.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "from Bio.Seq import Seq\n",
    "#from Bio.Alphabet import generic_dna, generic_protein (only for old Python)\n",
    "from Bio import SeqIO\n",
    "\n",
    "from Bio import pairwise2\n",
    "from Bio.SubsMat import MatrixInfo\n",
    "\n",
    "import re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b1492f48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['cdna', 'cds', 'ENSE00003753508_exon', 'ENSE00002667911_exon', 'ENSE00002419584_exon', 'ENSE00003625790_exon', 'ENSE00003518480_exon', 'ENSE00003723991_exon', 'ENSE00003712342_exon', 'ENSE00003725258_exon', 'ENSE00003786593_exon', 'ENSE00003545950_exon', 'ENSE00002037735_exon', 'intron_1', 'intron_2', 'intron_3', 'intron_4', 'intron_5', 'intron_6', 'intron_7', 'intron_8', 'intron_9', 'intron_10', 'peptide', 'utr3', 'utr5'])\n",
      "CDS\n",
      "ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCAGACCTATGGAAACTACTTCCTGAAAACAACGTTCTGTCCCCCTTGCCGTCCCAAGCAATGGATGATTTGATGCTGTCCCCGGACGATATTGAACAATGGTTCACTGAAGACCCAGGTCCAGATGAAGCTCCCAGAATGCCAGAGGCTGCTCCCCCCGTGGCCCCTGCACCAGCAGCTCCTACACCGGCGGCCCCTGCACCAGCCCCCTCCTGGCCCCTGTCATCTTCTGTCCCTTCCCAGAAAACCTACCAGGGCAGCTACGGTTTCCGTCTGGGCTTCTTGCATTCTGGGACAGCCAAGTCTGTGACTTGCACGTACTCCCCTGCCCTCAACAAGATGTTTTGCCAACTGGCCAAGACCTGCCCTGTGCAGCTGTGGGTTGATTCCACACCCCCGCCCGGCACCCGCGTCCGCGCCATGGCCATCTACAAGCAGTCACAGCACATGACGGAGGTTGTGAGGCGCTGCCCCCACCATGAGCGCTGCTCAGATAGCGATGGTCTGGCCCCTCCTCAGCATCTTATCCGAGTGGAAGGAAATTTGCGTGTGGAGTATTTGGATGACAGAAACACTTTTCGACATAGTGTGGTGGTGCCCTATGAGCCGCCTGAGGTTGGCTCTGACTGTACCACCATCCACTACAACTACATGTGTAACAGTTCCTGCATGGGCGGCATGAACCGGAGGCCCATCCTCACCATCATCACACTGGAAGACTCCAGTGGTAATCTACTGGGACGGAACAGCTTTGAGGTGCGTGTTTGTGCCTGTCCTGGGAGAGACCGGCGCACAGAGGAAGAGAATCTCCGCAAGAAAGGGGAGCCTCACCACGAGCTGCCCCCAGGGAGCACTAAGCGAGCACTGCCCAACAACACCAGCTCCTCTCCCCAGCCAAAGAAGAAACCACTGGATGGAGAATATTTCACCCTTCAGATCCGTGGGCGTGAGCGCTTCGAGATGTTCCGAGAGCTGAATGAGGCCTTGGAACTCAAGGATGCCCAGGCTGGGAAGGAGCCAGGGGGGAGCAGGGCTCACTCCAGCCACCTGAAGTCCAAAAAGGGTCAGTCTACCTCCCGCCATAAAAAACTCATGTTCAAGACAGAAGGGCCTGACTCAGACTGA\n",
      "\n",
      "\n",
      "Complement\n",
      "TACCTCCTCGGCGTCAGTCTAGGATCGCAGCTCGGGGGAGACTCAGTCCTTTGTAAAAGTCTGGATACCTTTGATGAAGGACTTTTGTTGCAAGACAGGGGGAACGGCAGGGTTCGTTACCTACTAAACTACGACAGGGGCCTGCTATAACTTGTTACCAAGTGACTTCTGGGTCCAGGTCTACTTCGAGGGTCTTACGGTCTCCGACGAGGGGGGCACCGGGGACGTGGTCGTCGAGGATGTGGCCGCCGGGGACGTGGTCGGGGGAGGACCGGGGACAGTAGAAGACAGGGAAGGGTCTTTTGGATGGTCCCGTCGATGCCAAAGGCAGACCCGAAGAACGTAAGACCCTGTCGGTTCAGACACTGAACGTGCATGAGGGGACGGGAGTTGTTCTACAAAACGGTTGACCGGTTCTGGACGGGACACGTCGACACCCAACTAAGGTGTGGGGGCGGGCCGTGGGCGCAGGCGCGGTACCGGTAGATGTTCGTCAGTGTCGTGTACTGCCTCCAACACTCCGCGACGGGGGTGGTACTCGCGACGAGTCTATCGCTACCAGACCGGGGAGGAGTCGTAGAATAGGCTCACCTTCCTTTAAACGCACACCTCATAAACCTACTGTCTTTGTGAAAAGCTGTATCACACCACCACGGGATACTCGGCGGACTCCAACCGAGACTGACATGGTGGTAGGTGATGTTGATGTACACATTGTCAAGGACGTACCCGCCGTACTTGGCCTCCGGGTAGGAGTGGTAGTAGTGTGACCTTCTGAGGTCACCATTAGATGACCCTGCCTTGTCGAAACTCCACGCACAAACACGGACAGGACCCTCTCTGGCCGCGTGTCTCCTTCTCTTAGAGGCGTTCTTTCCCCTCGGAGTGGTGCTCGACGGGGGTCCCTCGTGATTCGCTCGTGACGGGTTGTTGTGGTCGAGGAGAGGGGTCGGTTTCTTCTTTGGTGACCTACCTCTTATAAAGTGGGAAGTCTAGGCACCCGCACTCGCGAAGCTCTACAAGGCTCTCGACTTACTCCGGAACCTTGAGTTCCTACGGGTCCGACCCTTCCTCGGTCCCCCCTCGTCCCGAGTGAGGTCGGTGGACTTCAGGTTTTTCCCAGTCAGATGGAGGGCGGTATTTTTTGAGTACAAGTTCTGTCTTCCCGGACTGAGTCTGACT\n",
      "\n",
      "\n",
      "Reverse complement\n",
      "TCAGTCTGAGTCAGGCCCTTCTGTCTTGAACATGAGTTTTTTATGGCGGGAGGTAGACTGACCCTTTTTGGACTTCAGGTGGCTGGAGTGAGCCCTGCTCCCCCCTGGCTCCTTCCCAGCCTGGGCATCCTTGAGTTCCAAGGCCTCATTCAGCTCTCGGAACATCTCGAAGCGCTCACGCCCACGGATCTGAAGGGTGAAATATTCTCCATCCAGTGGTTTCTTCTTTGGCTGGGGAGAGGAGCTGGTGTTGTTGGGCAGTGCTCGCTTAGTGCTCCCTGGGGGCAGCTCGTGGTGAGGCTCCCCTTTCTTGCGGAGATTCTCTTCCTCTGTGCGCCGGTCTCTCCCAGGACAGGCACAAACACGCACCTCAAAGCTGTTCCGTCCCAGTAGATTACCACTGGAGTCTTCCAGTGTGATGATGGTGAGGATGGGCCTCCGGTTCATGCCGCCCATGCAGGAACTGTTACACATGTAGTTGTAGTGGATGGTGGTACAGTCAGAGCCAACCTCAGGCGGCTCATAGGGCACCACCACACTATGTCGAAAAGTGTTTCTGTCATCCAAATACTCCACACGCAAATTTCCTTCCACTCGGATAAGATGCTGAGGAGGGGCCAGACCATCGCTATCTGAGCAGCGCTCATGGTGGGGGCAGCGCCTCACAACCTCCGTCATGTGCTGTGACTGCTTGTAGATGGCCATGGCGCGGACGCGGGTGCCGGGCGGGGGTGTGGAATCAACCCACAGCTGCACAGGGCAGGTCTTGGCCAGTTGGCAAAACATCTTGTTGAGGGCAGGGGAGTACGTGCAAGTCACAGACTTGGCTGTCCCAGAATGCAAGAAGCCCAGACGGAAACCGTAGCTGCCCTGGTAGGTTTTCTGGGAAGGGACAGAAGATGACAGGGGCCAGGAGGGGGCTGGTGCAGGGGCCGCCGGTGTAGGAGCTGCTGGTGCAGGGGCCACGGGGGGAGCAGCCTCTGGCATTCTGGGAGCTTCATCTGGACCTGGGTCTTCAGTGAACCATTGTTCAATATCGTCCGGGGACAGCATCAAATCATCCATTGCTTGGGACGGCAAGGGGGACAGAACGTTGTTTTCAGGAAGTAGTTTCCATAGGTCTGAAAATGTTTCCTGACTCAGAGGGGGCTCGACGCTAGGATCTGACTGCGGCTCCTCCAT\n",
      "\n",
      "\n",
      "Transcription (DNA --> RNA)\n",
      "AUGGAGGAGCCGCAGUCAGAUCCUAGCGUCGAGCCCCCUCUGAGUCAGGAAACAUUUUCAGACCUAUGGAAACUACUUCCUGAAAACAACGUUCUGUCCCCCUUGCCGUCCCAAGCAAUGGAUGAUUUGAUGCUGUCCCCGGACGAUAUUGAACAAUGGUUCACUGAAGACCCAGGUCCAGAUGAAGCUCCCAGAAUGCCAGAGGCUGCUCCCCCCGUGGCCCCUGCACCAGCAGCUCCUACACCGGCGGCCCCUGCACCAGCCCCCUCCUGGCCCCUGUCAUCUUCUGUCCCUUCCCAGAAAACCUACCAGGGCAGCUACGGUUUCCGUCUGGGCUUCUUGCAUUCUGGGACAGCCAAGUCUGUGACUUGCACGUACUCCCCUGCCCUCAACAAGAUGUUUUGCCAACUGGCCAAGACCUGCCCUGUGCAGCUGUGGGUUGAUUCCACACCCCCGCCCGGCACCCGCGUCCGCGCCAUGGCCAUCUACAAGCAGUCACAGCACAUGACGGAGGUUGUGAGGCGCUGCCCCCACCAUGAGCGCUGCUCAGAUAGCGAUGGUCUGGCCCCUCCUCAGCAUCUUAUCCGAGUGGAAGGAAAUUUGCGUGUGGAGUAUUUGGAUGACAGAAACACUUUUCGACAUAGUGUGGUGGUGCCCUAUGAGCCGCCUGAGGUUGGCUCUGACUGUACCACCAUCCACUACAACUACAUGUGUAACAGUUCCUGCAUGGGCGGCAUGAACCGGAGGCCCAUCCUCACCAUCAUCACACUGGAAGACUCCAGUGGUAAUCUACUGGGACGGAACAGCUUUGAGGUGCGUGUUUGUGCCUGUCCUGGGAGAGACCGGCGCACAGAGGAAGAGAAUCUCCGCAAGAAAGGGGAGCCUCACCACGAGCUGCCCCCAGGGAGCACUAAGCGAGCACUGCCCAACAACACCAGCUCCUCUCCCCAGCCAAAGAAGAAACCACUGGAUGGAGAAUAUUUCACCCUUCAGAUCCGUGGGCGUGAGCGCUUCGAGAUGUUCCGAGAGCUGAAUGAGGCCUUGGAACUCAAGGAUGCCCAGGCUGGGAAGGAGCCAGGGGGGAGCAGGGCUCACUCCAGCCACCUGAAGUCCAAAAAGGGUCAGUCUACCUCCCGCCAUAAAAAACUCAUGUUCAAGACAGAAGGGCCUGACUCAGACUGA\n",
      "\n",
      "\n",
      "Translation (DNA --> RNA --> AA)\n",
      "MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD*\n",
      "\n",
      "\n",
      "Translation (implicit transcription, DNA --> AA)\n",
      "MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD*\n",
      "\n",
      "\n",
      "Complement translation\n",
      "YLLGVSLGSQLGGDSVLCKSLDTFDEGLLLQDRGNGRVRYLLNYDRGLL*LVTK*LLGPGLLRGSYGLRRGGHRGRGRRGCGRRGRGRGRTGDSRRQGRVFWMVPSMPKADPKNVRPCRFRH*TCMRGRELFYKTVDRFWTGHVDTQLRCGGGPWAQARYR*MFVSVVYCLQHSATGVVLATSLSLPDRGGVVE*AHLPLNAHLINLLSL*KAVSHHHGILGGLQPRLTWW*VMLMYTLSRTYPPYLASG*EW**CDLLRSPLDDPALSKLHAQTRTGPSLAACLLLLEAFFPLGVVLDGGPS*FARDGLLWSRRGVGFFFGDLPLIKWEV*APALAKLYKALDLLRNLEFLRVRPFLGPPSSRVRSVDFRFFPVRWRAVFFEYKFCLPGLSLT\n",
      "\n",
      "\n",
      "Reverse complement translation\n",
      "SV*VRPFCLEHEFFMAGGRLTLFGLQVAGVSPAPPWLLPSLGILEFQGLIQLSEHLEALTPTDLKGEIFSIQWFLLWLGRGAGVVGQCSLSAPWGQLVVRLPFLAEILFLCAPVSPRTGTNTHLKAVPSQ*ITTGVFQCDDGEDGPPVHAAHAGTVTHVVVVDGGTVRANLRRLIGHHHTMSKSVSVIQILHTQISFHSDKMLRRGQTIAI*AALMVGAAPHNLRHVL*LLVDGHGADAGAGRGCGINPQLHRAGLGQLAKHLVEGRGVRASHRLGCPRMQEAQTETVAALVGFLGRDRR*QGPGGGWCRGRRCRSCWCRGHGGSSLWHSGSFIWTWVFSEPLFNIVRGQHQIIHCLGRQGGQNVVFRK*FP*V*KCFLTQRGLDARI*LRLLH\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Parse FASTA\n",
    "seq_records = SeqIO.parse(\"data/ensembl_p53.fasta\", \"fasta\")\n",
    "\n",
    "data = {}  # {molecule_type: [seq_record, ...]}\n",
    "for rec in seq_records:\n",
    "#     print(rec)\n",
    "    # Parse the sequence type (peptide, cds, cdna, utr5, utr3, <x>_exon, intron_<x>)\n",
    "    seq_type = \"_\".join(rec.description.split(':')[0].split()[1:])\n",
    "    data[seq_type] = rec\n",
    "    \n",
    "\n",
    "print(data.keys())    \n",
    "    \n",
    "# for record_type in data:\n",
    "#     print(record_type, data[record_type], len(data[record_type].seq), sep=\"\\n\")\n",
    "#     print()\n",
    "\n",
    "\n",
    "# Turn a nucleotide sequence into a protein sequence\n",
    "print(\"CDS\\n{}\\n\\n\".format(data[\"cds\"].seq))\n",
    "\n",
    "print(\"Complement\\n{}\\n\\n\".format(data[\"cds\"].seq.complement()))\n",
    "\n",
    "print(\"Reverse complement\\n{}\\n\\n\".format(data[\"cds\"].seq.reverse_complement()))\n",
    "\n",
    "print(\"Transcription (DNA --> RNA)\\n{}\\n\\n\".format(data[\"cds\"].seq.transcribe()))\n",
    "\n",
    "print(\"Translation (DNA --> RNA --> AA)\\n{}\\n\\n\".format(data[\"cds\"].seq.transcribe().translate()))\n",
    "\n",
    "print(\"Translation (implicit transcription, DNA --> AA)\\n{}\\n\\n\".format(data[\"cds\"].seq.translate()))\n",
    "\n",
    "print(\"Complement translation\\n{}\\n\\n\".format(data[\"cds\"].seq.complement().transcribe().translate()))\n",
    "\n",
    "print(\"Reverse complement translation\\n{}\\n\\n\".format(data[\"cds\"].seq.reverse_complement().transcribe().translate()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "292a141f",
   "metadata": {},
   "source": [
    "# Map exons to the cDNA and CDS\n",
    "\n",
    "Find position of exons using two strategies:\n",
    "\n",
    "- Exact match (ex. regex)\n",
    "- Alignment (pairwise2 module)\n",
    "\n",
    "### Pairwise 2 module\n",
    "The module includes different methods. The names are **\\<algorithm_type\\>\\<match_parameter\\>\\<gap_parameter\\>**, ex. **localxx()**\n",
    "\n",
    "The match parameters are (CODE  DESCRIPTION):\n",
    "- **x**     No parameters. Identical characters have score of 1, otherwise 0\n",
    "- **m**     A match score is the score of identical chars, otherwise mismatch score\n",
    "- **d**     A dictionary returns the score of any pair of characters\n",
    "- **c**     A callback function returns scores\n",
    "\n",
    "The gap penalty parameters are (CODE  DESCRIPTION):\n",
    "- **x**     No gap penalties\n",
    "- **s**     Same open and extend gap penalties for both sequences\n",
    "- **d**     The sequences have different open and extend gap penalties\n",
    "- **c**     A callback function returns the gap penalties\n",
    "\n",
    "The align method return a list of tuple corresponing to alignments, each tuple has:\n",
    "    \n",
    "- seq1 aligned\n",
    "- seq2 aligned\n",
    "- score\n",
    "- alignment start position\n",
    "- alignment end position\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f3b82f1d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('ENST00000269305.9 ENSE00003753508 exon:protein_coding', (0, 114))\n",
      "('ENST00000269305.9 ENSE00002667911 exon:protein_coding', (114, 216))\n",
      "('ENST00000269305.9 ENSE00002419584 exon:protein_coding', (216, 238))\n",
      "('ENST00000269305.9 ENSE00003625790 exon:protein_coding', (238, 517))\n",
      "('ENST00000269305.9 ENSE00003518480 exon:protein_coding', (517, 701))\n",
      "('ENST00000269305.9 ENSE00003723991 exon:protein_coding', (701, 814))\n",
      "('ENST00000269305.9 ENSE00003712342 exon:protein_coding', (814, 924))\n",
      "('ENST00000269305.9 ENSE00003725258 exon:protein_coding', (924, 1061))\n",
      "('ENST00000269305.9 ENSE00003786593 exon:protein_coding', (1061, 1135))\n",
      "('ENST00000269305.9 ENSE00003545950 exon:protein_coding', (1135, 1242))\n",
      "('ENST00000269305.9 ENSE00002037735 exon:protein_coding', (1242, 2512))\n"
     ]
    }
   ],
   "source": [
    "# Use exact match against the cDNA (regex)\n",
    "matches = []\n",
    "for k in sorted(data.keys()):\n",
    "    if \"exon\" in k:\n",
    "        match = re.search(str(data[k].seq), str(data[\"cdna\"].seq))\n",
    "        matches.append((data[k].description, match.span() if match else (0, 0)))\n",
    "for match in sorted(matches, key=lambda x: x[1][0]):\n",
    "    print(match)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b1936dfa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('ENST00000269305.9 ENSE00002037735 exon:protein_coding', (0, 0))\n",
      "('ENST00000269305.9 ENSE00002667911 exon:protein_coding', (0, 0))\n",
      "('ENST00000269305.9 ENSE00003753508 exon:protein_coding', (0, 0))\n",
      "('ENST00000269305.9 ENSE00002419584 exon:protein_coding', (74, 96))\n",
      "('ENST00000269305.9 ENSE00003625790 exon:protein_coding', (96, 375))\n",
      "('ENST00000269305.9 ENSE00003518480 exon:protein_coding', (375, 559))\n",
      "('ENST00000269305.9 ENSE00003723991 exon:protein_coding', (559, 672))\n",
      "('ENST00000269305.9 ENSE00003712342 exon:protein_coding', (672, 782))\n",
      "('ENST00000269305.9 ENSE00003725258 exon:protein_coding', (782, 919))\n",
      "('ENST00000269305.9 ENSE00003786593 exon:protein_coding', (919, 993))\n",
      "('ENST00000269305.9 ENSE00003545950 exon:protein_coding', (993, 1100))\n"
     ]
    }
   ],
   "source": [
    "# Use exact match against the CDS (regex)\n",
    "matches = []\n",
    "for k in sorted(data.keys()):\n",
    "    if \"exon\" in k:\n",
    "        match = re.search(str(data[k].seq), str(data[\"cds\"].seq))\n",
    "        matches.append((data[k].description, match.span() if match else (0, 0)))\n",
    "for match in sorted(matches, key=lambda x: x[1][0]):\n",
    "    print(match)\n",
    "\n",
    "# Why some exons do not align with the cds?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f0b9e325",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{('W', 'F'): -1, ('L', 'R'): -1, ('I', 'I'): 6, ('Q', 'Q'): 6, ('W', 'N'): -1, ('V', 'I'): -1, ('H', 'T'): -1, ('H', 'P'): -1, ('W', 'V'): -1, ('Q', 'E'): -1, ('W', 'R'): -1, ('Q', 'A'): -1, ('H', 'H'): 6, ('H', 'D'): -1, ('L', 'N'): -1, ('Y', 'M'): -1, ('Y', 'I'): -1, ('Y', 'E'): -1, ('E', 'S'): -1, ('Y', 'A'): -1, ('Y', 'Y'): 6, ('T', 'C'): -1, ('E', 'C'): -1, ('Y', 'Q'): -1, ('E', 'G'): -1, ('V', 'A'): -1, ('C', 'C'): 6, ('M', 'R'): -1, ('P', 'T'): -1, ('V', 'E'): -1, ('P', 'P'): 6, ('I', 'T'): -1, ('K', 'S'): -1, ('R', 'G'): -1, ('I', 'P'): -1, ('R', 'C'): -1, ('A', 'T'): -1, ('K', 'K'): 6, ('A', 'P'): -1, ('V', 'M'): -1, ('I', 'D'): -1, ('K', 'C'): -1, ('K', 'G'): -1, ('R', 'S'): -1, ('F', 'Q'): -1, ('F', 'A'): -1, ('V', 'V'): 6, ('M', 'N'): -1, ('F', 'E'): -1, ('D', 'N'): -1, ('F', 'I'): -1, ('F', 'M'): -1, ('M', 'S'): -1, ('S', 'S'): 6, ('L', 'Q'): -1, ('W', 'E'): -1, ('W', 'A'): -1, ('W', 'M'): -1, ('H', 'S'): -1, ('W', 'I'): -1, ('S', 'C'): -1, ('L', 'A'): -1, ('L', 'E'): -1, ('W', 'Q'): -1, ('H', 'G'): -1, ('Q', 'N'): -1, ('H', 'C'): -1, ('L', 'M'): -1, ('W', 'Y'): -1, ('Y', 'N'): -1, ('E', 'P'): -1, ('Y', 'F'): -1, ('E', 'T'): -1, ('A', 'A'): 6, ('I', 'N'): -1, ('G', 'A'): -1, ('Y', 'V'): -1, ('E', 'D'): -1, ('W', 'H'): -1, ('Y', 'R'): -1, ('M', 'Q'): -1, ('P', 'S'): -1, ('R', 'H'): -1, ('A', 'C'): -1, ('R', 'D'): -1, ('K', 'P'): -1, ('L', 'D'): -1, ('K', 'T'): -1, ('V', 'N'): -1, ('M', 'A'): -1, ('K', 'H'): -1, ('V', 'R'): -1, ('P', 'C'): -1, ('M', 'E'): -1, ('A', 'S'): -1, ('T', 'T'): 6, ('R', 'T'): -1, ('I', 'G'): -1, ('R', 'P'): -1, ('K', 'D'): -1, ('I', 'C'): -1, ('F', 'R'): -1, ('F', 'V'): -1, ('L', 'C'): -1, ('F', 'F'): 6, ('D', 'A'): -1, ('F', 'N'): -1, ('W', 'D'): -1, ('L', 'P'): -1, ('Q', 'S'): -1, ('N', 'C'): -1, ('N', 'G'): -1, ('H', 'N'): -1, ('W', 'T'): -1, ('Q', 'G'): -1, ('W', 'P'): -1, ('Q', 'C'): -1, ('N', 'S'): -1, ('L', 'H'): -1, ('L', 'L'): 6, ('G', 'T'): -1, ('M', 'M'): 6, ('G', 'P'): -1, ('Y', 'K'): -1, ('Y', 'G'): -1, ('Y', 'C'): -1, ('E', 'A'): -1, ('E', 'E'): 6, ('Y', 'S'): -1, ('M', 'P'): -1, ('V', 'C'): -1, ('M', 'T'): -1, ('V', 'G'): -1, ('R', 'E'): -1, ('V', 'K'): -1, ('K', 'Q'): -1, ('R', 'A'): -1, ('I', 'R'): -1, ('N', 'A'): -1, ('V', 'S'): -1, ('M', 'D'): -1, ('M', 'H'): -1, ('K', 'A'): -1, ('R', 'Q'): -1, ('K', 'E'): -1, ('F', 'S'): -1, ('I', 'K'): -1, ('D', 'P'): -1, ('D', 'T'): -1, ('I', 'M'): -1, ('F', 'C'): -1, ('W', 'L'): -1, ('F', 'G'): -1, ('F', 'K'): -1, ('F', 'T'): -1, ('D', 'D'): 6, ('Q', 'T'): -1, ('W', 'G'): -1, ('Q', 'P'): -1, ('W', 'C'): -1, ('W', 'K'): -1, ('H', 'Q'): -1, ('Q', 'D'): -1, ('W', 'W'): 6, ('V', 'L'): -1, ('L', 'G'): -1, ('W', 'S'): -1, ('L', 'K'): -1, ('N', 'P'): -1, ('H', 'E'): -1, ('N', 'T'): -1, ('H', 'A'): -1, ('Y', 'L'): -1, ('Y', 'H'): -1, ('G', 'S'): -1, ('Y', 'D'): -1, ('V', 'Q'): -1, ('L', 'T'): -1, ('G', 'G'): 6, ('G', 'C'): -1, ('E', 'N'): -1, ('Y', 'T'): -1, ('Y', 'P'): -1, ('R', 'N'): -1, ('V', 'D'): -1, ('K', 'R'): -1, ('V', 'H'): -1, ('I', 'Q'): -1, ('V', 'P'): -1, ('M', 'C'): -1, ('K', 'N'): -1, ('V', 'T'): -1, ('M', 'G'): -1, ('T', 'S'): -1, ('I', 'E'): -1, ('M', 'K'): -1, ('I', 'A'): -1, ('N', 'N'): 6, ('R', 'R'): 6, ('F', 'P'): -1, ('L', 'I'): -1, ('I', 'S'): -1, ('D', 'S'): -1, ('L', 'S'): -1, ('I', 'H'): -1, ('F', 'D'): -1, ('D', 'C'): -1, ('F', 'H'): -1, ('D', 'G'): -1, ('F', 'L'): -1}\n",
      "ENSE00002667911_exon start:  28 end: 102 score: 444.0 score_average:6.00 \n",
      "ENSE00002037735_exon start:  60 end:1211 score:1124.0 score_average:0.98 \n",
      "ENSE00002419584_exon start:  74 end:  96 score: 132.0 score_average:6.00 \n",
      "ENSE00003625790_exon start:  96 end: 375 score:1674.0 score_average:6.00 \n",
      "ENSE00003753508_exon start: 170 end: 280 score: 198.0 score_average:1.80 \n",
      "ENSE00003518480_exon start: 375 end: 559 score:1104.0 score_average:6.00 \n",
      "ENSE00003723991_exon start: 559 end: 672 score: 678.0 score_average:6.00 \n",
      "ENSE00003712342_exon start: 672 end: 782 score: 660.0 score_average:6.00 \n",
      "ENSE00003725258_exon start: 782 end: 919 score: 822.0 score_average:6.00 \n",
      "ENSE00003786593_exon start: 919 end: 993 score: 444.0 score_average:6.00 \n",
      "ENSE00003545950_exon start: 993 end:1100 score: 642.0 score_average:6.00 \n"
     ]
    }
   ],
   "source": [
    "# Use alignments against the CDS \n",
    "# Provide some indexes to evaluate alignment quality\n",
    "# The score is given by combining match (substitution matrix) and penalty (gap) scores\n",
    "\n",
    "alignments = []\n",
    "for k in sorted(data.keys()):\n",
    "    if \"exon\" in k:\n",
    "        \n",
    "        # Global alignment, match 1, non-match 0, no gap penalties\n",
    "#         alignment = pairwise2.align.globalxx(data[\"cds\"].seq, data[k].seq, one_alignment_only=True)[0]\n",
    "        \n",
    "    \n",
    "        # Local alignment\n",
    "        \n",
    "        # match 1, non-match 0, no gap penalties \n",
    "#         alignment = pairwise2.align.localxx(data[\"cds\"].seq, data[k].seq, one_alignment_only=True)[0]\n",
    "        \n",
    "        # match 1, non-match 0, gap open -10, gap extention -0.5\n",
    "#         alignment = pairwise2.align.localxs(data[\"cds\"].seq, data[k].seq, -10, -0.5, one_alignment_only=True)[0]\n",
    "\n",
    "        # identity matrix (match 6, non-match -1), gap open -200, gap extention -0.5\n",
    "        alignment = pairwise2.align.localds(data[\"cds\"].seq, data[k].seq, MatrixInfo.ident, -200, -0.5, one_alignment_only=True)[0]\n",
    "\n",
    "        seq1_aligned, seq2_aligned, score, alignment_begin, alignment_end = alignment\n",
    "        \n",
    "        # Calculate average score per alignment position\n",
    "        score_average = score / float(alignment_end - alignment_begin)\n",
    "\n",
    "        alignments.append((k, alignment_begin, alignment_end, score, score_average, alignment))\n",
    "\n",
    "        \n",
    "print(MatrixInfo.ident)\n",
    "        \n",
    "# Sort by starting position of the alignment\n",
    "for k, alignment_begin, alignment_end, score, score_average, alignment in sorted(alignments, \n",
    "                                                                                 key=lambda ele: ele[1]):\n",
    "\n",
    "#     print(pairwise2.format_alignment(*alignment))\n",
    "    print(\"{} start:{:>4} end:{:>4} score:{:>6} score_average:{:>4.2f} \".format(k, alignment_begin, alignment_end, \n",
    "                                                                     score, score_average))\n",
    "   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d652d48c",
   "metadata": {},
   "source": [
    "## Questions\n",
    "\n",
    "- Why even when using high penalties I don't get a good alignment for the CDS/UTR overlapping exon?\n",
    "\n",
    "- Why the non overlapping exon has non-zero average score?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d76fc69d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ENSE00002037735_exon 0   129  197  47.5   366  37\n",
      "ENSE00002667911_exon 0   351  374  26.0    33   7\n",
      "ENSE00003753508_exon 0     8   22  27.0    37   7\n",
      "ENSE00002037735_exon 1   367  393 137.0    26  26\n",
      "ENSE00002667911_exon 1     9   33 130.0    24  24\n",
      "ENSE00003753508_exon 1   152  157  20.0    35   4\n",
      "ENSE00002037735_exon 2    74  221  45.0   345  45\n",
      "ENSE00002667911_exon 2    79  108  28.0    32  10\n",
      "ENSE00003753508_exon 2    98  133  27.5    34  12\n"
     ]
    }
   ],
   "source": [
    "# Translate the exons\n",
    "# Identify those which are coding for protein fragments\n",
    "# Align translated exons against the protein (AA) sequence\n",
    "# Note, the DNA of the exons can have a non-zero phase\n",
    "\n",
    "# ENSE00003753508 (1, 3'-UTR, non coding)\n",
    "# ENSE00002667911 (2, 3'-UTR + coding)\n",
    "# ...\n",
    "# ENSE00002037735 (11, 5'-UTR + coding)\n",
    "\n",
    "\n",
    "# Use alignments (different parameters)\n",
    "alignments = []\n",
    "for k in sorted(data.keys()):\n",
    "    if \"exon\" in k:\n",
    "        \n",
    "        # Try all possible reading frames\n",
    "        for i in range(0, 3):\n",
    "        \n",
    "            # Translate the exon (remove stop codon symbols)\n",
    "            exon_peptide = data[k].seq[i:].translate(stop_symbol=\"\")\n",
    "\n",
    "            alignment = pairwise2.align.localds(data[\"peptide\"].seq, exon_peptide, MatrixInfo.blosum62, -10, -0.5, one_alignment_only=True)[0]\n",
    "\n",
    "            seq1_aligned, seq2_aligned, score, alignment_begin, alignment_end = alignment\n",
    "            \n",
    "            \n",
    "            \n",
    "            identity = 0\n",
    "            matches = 0\n",
    "            for l1, l2 in zip(seq1_aligned, seq2_aligned):\n",
    "                if l1 != '-' and l2 != '-':\n",
    "                    matches += 1\n",
    "                    if l1 == l2:\n",
    "                        identity += 1\n",
    "\n",
    "            alignments.append((k, i, alignment_begin, alignment_end, score, identity, matches, alignment))\n",
    "\n",
    "#             print(seq1_aligned, seq2_aligned, identity, matches, sep=\"\\n\")\n",
    "#             break\n",
    "#         break\n",
    "\n",
    "# Print the Blosum62 matrix\n",
    "# print(sorted(MatrixInfo.blosum62.items(), key=lambda x: x[1], reverse=True))\n",
    "        \n",
    "# Sort by starting position of the alignment\n",
    "for k, i, alignment_begin, alignment_end, score, identity, matches, alignment in sorted(alignments, \n",
    "                                                                                 key=lambda ele: ele[1]):\n",
    "    if k in ['ENSE00002037735_exon', 'ENSE00003753508_exon', 'ENSE00002667911_exon']:\n",
    "        \n",
    "        \n",
    "        # Exon id, start, end, score, score_average, alignment_length, exon_aligment_coverage\n",
    "        print(\"{} {} {:>5}{:>5}{:>6.1f}{:>6}{:>4}\".format(k, i, alignment_begin, alignment_end, \n",
    "                                                score, \n",
    "                                                matches, \n",
    "                                                identity))\n",
    "\n",
    "# #         print(alignment)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}