notebooks/07_module_biopython.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4f931802",
   "metadata": {},
   "source": [
    "# Table of Content <a id='toc'></a>\n",
    "\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[Module 7 - Biopython](#0)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Broad examples of what we can do with Biopython](#1)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Main concept of Biopython](#2)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Objects](#3)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Help - important!](#4)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Accessing Online Databases](#17)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Bio.SeqIO Module](#11)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Reading sequence records from files](#12)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Writing to sequence files](#16)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[`Seq` objects and 'Biological' methods](#10)\n",
    "\n",
    "\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Now, you can try to solve the exercises.](#18)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[Additional Theory](#19)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[On `SeqRecord` iterators and processing large sequence files with `SeqIO`](#20)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[More on `SeqFeature` objects](#21)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b6d785e",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "# Module 7 - Biopython <a id='0'></a>\n",
    "----------------------------------\n",
    "\n",
    "<br>\n",
    "\n",
    "Let's consider a typical task in bioinformatics: we are interested in finding the GC% for some sequences. For this we would need i) open the file, ii) parse and extract the sequence information, and iii) calculate and report their GC%. In a very simple example, we could write something like following: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c86d519",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_seq = {}\n",
    "identifier = None\n",
    "with open(\"data/my_sequences.fa\") as input_file:\n",
    "    for line in input_file:\n",
    "        if line[0] == \">\":\n",
    "            identifier, *description = line.strip()[1:].split(\" \", 1)\n",
    "            my_seq[identifier] = \"\"\n",
    "        else:\n",
    "            my_seq[identifier] += line.strip().upper()\n",
    "for identifier in my_seq:\n",
    "    gc_count = my_seq[identifier].count(\"G\") + my_seq[identifier].count(\"C\")\n",
    "    gc_pc = gc_count / len(my_seq[identifier]) * 100\n",
    "    print(\"GC% for {} is {:.2f}\".format(identifier, gc_pc))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70ba9a69",
   "metadata": {},
   "source": [
    "Now let's see how it would look like if we used Biopython:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6719291b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "from Bio.SeqUtils import GC\n",
    "for rec in SeqIO.parse(\"data/my_sequences.fa\", \"fasta\"):\n",
    "    print(\"GC% for {} is {:.2f}\".format(rec.id, GC(rec.seq)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a6aa630",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "## Broad examples of what we can do with Biopython <a id='1'></a>\n",
    "\n",
    "- Sequence analysis\n",
    "  - Motif\n",
    "  - Search: HMMs, Alignments\n",
    "  - Restriction\n",
    "- Structures\n",
    "  - SCOP\n",
    "  - PDB\n",
    "- Database query\n",
    "- Phylogeny\n",
    "- Pathway\n",
    "- And more ...\n",
    "\n",
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "## Main concept of Biopython <a id='2'></a>\n",
    "\n",
    "- I/O interface and parsing abilities for bioinformatics files/DBs\n",
    "  - Blast\n",
    "  - FASTA\n",
    "  - PubMed\n",
    "  - SwissProt\n",
    "- Efficient and practical data-structures for bioinformatics data\n",
    " - Sequences\n",
    " - Alignments\n",
    " - Structures\n",
    "- Methods implementing bioinformatics analysis\n",
    " - Translation\n",
    " - Classification\n",
    " - Phylogeny trees\n",
    "- Interface to common bioinformatics programs\n",
    "  - Standalone Blast\n",
    "  - ClustalW\n",
    "  - EMBOSS\n",
    "\n",
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "## Objects <a id='3'></a>\n",
    "\n",
    "In Biopython there are many modules and each module contains several major new data types. Objects created with these types serve specific purposes in the above mentioned functionalities. We will focus on sequence analysis; some new objects we will discover are:\n",
    "\n",
    "- `Seq`\n",
    "- `SeqRecord`\n",
    "\n",
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "## Help - important! <a id='4'></a>\n",
    "\n",
    "- A relatively detailed tutorial: http://biopython.org/DIST/docs/tutorial/Tutorial.html\n",
    "- Help on certain concepts and modules: https://biopython.org/wiki/Category%3AWiki_Documentation\n",
    "\n",
    "In this notebook we will cover basics of the Biopython package:\n",
    "    - Seq object and alphabets\n",
    "    - Bio.SeqIO module and SeqRecord object\n",
    "    - Interacting with external DBs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9824219a",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "## Accessing Online Databases <a id='17'></a>\n",
    "\n",
    "Biophyton introduces multiple means to interact with commonly used bioinformatics databases over the internet. We will cover the basics of accessing online databases with two databases:\n",
    "\n",
    "1. NCBI's Entrez\n",
    "2. UniProt/SwissProt's ExPASy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da3e3e5b",
   "metadata": {},
   "source": [
    "The Entrez can be queried using the `Bio.Entrez` module. There many different functionalities, but today we will only cover the `.efetch()` method, which we can use to fetch records from Entrez over the internet.\n",
    "\n",
    "Let's discover this functionality with an example.\n",
    "\n",
    "<br>\n",
    "\n",
    "Imagine that you used blast using some sequence as query and got significant hits with the following Entrez sequence ids :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2880f0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "hits_ids = ['KAJ0778808.1',\n",
    "             'PWA60940.1',\n",
    "             'KAI7755456.1',\n",
    "             'GEU62861.1',\n",
    "             'AIG92846.1',\n",
    "             'AIW00681.1',\n",
    "             'XP_024994640.1',\n",
    "             'KAI3824653.1',\n",
    "             'KAD3641182.1',\n",
    "             'KAI7751279.1']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c384db52",
   "metadata": {},
   "source": [
    "We would like to gather some information regarding these sequences.\n",
    "\n",
    "First we will just check the connection and look-up the accessible databases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65fcf4bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import Entrez\n",
    "from Bio import SeqIO\n",
    "import time \n",
    "\n",
    "Entrez.email = \"...\"  # write you e-mail here. Always tell NCBI who you are\n",
    "stream = Entrez.einfo()\n",
    "result = stream.read()\n",
    "stream.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcc419ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2825848",
   "metadata": {},
   "source": [
    "As a result we received some XML  which contains a `<DbList>` which contains a number of `<DbName>`.\n",
    "\n",
    "The batabase we are interested in is `protein`. \n",
    "\n",
    "Let's fetch the first ID:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a07dfb74",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "pid = hits_ids[0] # we will fetch the 1st ID in our list\n",
    "stream = Entrez.efetch(db=\"protein\", id=pid, rettype=\"gb\", retmode=\"text\")\n",
    "\n",
    "## we receive the data as a text stream in the genbank forma, which we need to interpret:\n",
    "record = SeqIO.read(stream, \"genbank\")\n",
    "\n",
    "stream.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57f07da9",
   "metadata": {},
   "source": [
    "What we received is a `SeqRecord` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8ee2b3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "record"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00d46fa0",
   "metadata": {},
   "source": [
    "Biopython often stores complex objects using custom types which store their elements in different attributes.\n",
    "\n",
    "Let's explore this SeqRecord object (you can also browse the [online documentation](https://biopython.org/docs/latest/api/Bio.SeqRecord.html?highlight=seqrecord#module-Bio.SeqRecord))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed46bc9b",
   "metadata": {},
   "source": [
    "`SeqRecord` contains many attributes, some of which can be rather complex. We will further examine most commonly used ones:\n",
    "\n",
    "- seq\n",
    "- id\n",
    "- name\n",
    "- description\n",
    "- annotations\n",
    "- features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7bdbced",
   "metadata": {},
   "outputs": [],
   "source": [
    "print( record.seq )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b15d4e8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "record.id , record.description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "924caad7",
   "metadata": {},
   "outputs": [],
   "source": [
    "record.annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82211846",
   "metadata": {},
   "outputs": [],
   "source": [
    "record.features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eda453a",
   "metadata": {},
   "source": [
    "Now we can download the records for all of our hits, and save the results in genbank files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dce2dee",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import os\n",
    "\n",
    "i = 0\n",
    "for pid in hits_ids:\n",
    "    i+=1\n",
    "    \n",
    "    print(f\"{i}/{len(hits_ids)}\")\n",
    "        \n",
    "    stream = Entrez.efetch(db=\"protein\", id=pid, rettype=\"gb\", retmode=\"text\")\n",
    "    record = SeqIO.read(stream, \"genbank\")\n",
    "    stream.close()\n",
    "        \n",
    "    ## writing the SeqRecord to a genbank file\n",
    "    SeqIO.write(record , pid + '.gb' , \"genbank\")\n",
    "    \n",
    "    ## wait to be sure to respect ncbi API guidelines\n",
    "    time.sleep(0.3)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f3dd02c",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "## Bio.SeqIO Module  <a id='11'></a>\n",
    "\n",
    "The `Bio.SeqIO` module provides a simple and uniform interface to read and parse from and write to various sequence file formats (including multiple sequence alignments). The `Bio.SeqIO` module supports a large number of sequence file formats, including *fasta*, *fastq*, and *genbank* and many more, which can be found [here](https://biopython.org/wiki/SeqIO). All sequences in this method will be accessed as `SeqRecord` objects.\n",
    "\n",
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "### Reading sequence records from files <a id='12'></a>\n",
    "\n",
    "The main function of the module is `.parse()`, which takes a file handle (or filename) and format name, and returns a `SeqRecord` iterator. An iterator is an object that can be iterated upon, meaning that you can traverse through all the values, one by one. In this case `Bio.SeqIO.parse()` method returns an iterator of `SeqRecord` objects extracted and parsed from the input file. We can then iterate over the records and process them in a very efficient manner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45d52734",
   "metadata": {},
   "outputs": [],
   "source": [
    "## reading from a single file\n",
    "\n",
    "records= []\n",
    "\n",
    "## SeqIO.parse tries to look for several records in a single file, so we need to put it in a for loop\n",
    "for rec  in SeqIO.parse( hits_ids[0] + '.gb' , format = 'genbank' ) :\n",
    "    records.append( rec )\n",
    "\n",
    "rec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9ace9f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "## reading several files\n",
    "\n",
    "records = []\n",
    "\n",
    "for h in hits_ids:\n",
    "    records += list( SeqIO.parse( h + \".gb\", \"genbank\") ) # using the list() function to go through all records\n",
    "records"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3426abd2",
   "metadata": {},
   "source": [
    "#### micro-exercise\n",
    "\n",
    "Similarly parse the **example.fa** file and compare the `SeqRecord` objects to the above cell. Which attributes are identical? What are the differences? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b12ff34a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "188388c9",
   "metadata": {},
   "source": [
    "We can also create a `SeqRecord` directly. We need to import the class from its module, `Bio.SeqRecord`, first. In order to create its sequence we will also need the `Bio.Seq.Seq` or a similar `Seq` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edaf7a6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "\n",
    "record = SeqRecord(Seq(\"ACGGCTATCTGAGGACTACGAGCATCATCGAG\"),\n",
    "                   id=\"my_seq_ID\", name=\"made-up sequence\",\n",
    "                   description=\"Just some randomly typed ATGCs\")\n",
    "print(record)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c487a27",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "afbe4dd2",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "### Writing to sequence files <a id='16'></a>\n",
    "So far we have seen how to parse sequence files and work with sequence records. In many applications, we will also need to write our processed sequence records back into a standard sequence file. For this purpose, the `Bio.SeqIO` module has a `.write()` function. It can be thought as the reverse of the `.parse()` method we have learnt above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff514f28",
   "metadata": {},
   "outputs": [],
   "source": [
    "SeqIO.write(records, 'records.fasta', \"fasta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6940f4c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "! head -n 20 records.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a1095c5",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "## `Seq` objects and 'Biological' methods <a id='10'></a>\n",
    "\n",
    "Seq object are a flexible encapsulator for _sequence-like_ strings. \n",
    "\n",
    "In its essence, it behaves very much like a classical `str`, but it also has some\n",
    "specific 'biological' methods, such has:\n",
    "\n",
    "- complement()\n",
    "- reverse_complement()\n",
    "- transcribe()\n",
    "- back_transcribe()\n",
    "- translate()\n",
    "\n",
    "These methods do not *require* any other extra arguments and they return a new `Seq` object."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "835575bc",
   "metadata": {},
   "source": [
    "#### Complementing and reverse complementing sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb80d902",
   "metadata": {},
   "outputs": [],
   "source": [
    "intronless_dna = Seq(\n",
    "    'CACCTCTGGAGCGGACTTATTTACCAAGCATTGGAGGAATATCGTAGGTAAAAATGCCTA'\n",
    "    'TAGGATCCAAAGAGAGGCCAACATTTTTTGAAATTTTTAAGACACGCTGCAACAAAGCA')\n",
    "print(\"intronless_dna:\", intronless_dna)\n",
    "\n",
    "\n",
    "# If we need to find its complement sequence:\n",
    "\n",
    "compl_seq = intronless_dna.complement()\n",
    "print()\n",
    "print(\"compl_seq:\", compl_seq)\n",
    "\n",
    "# If we would need to find the reverse complement of our sequence, we can simply\n",
    "\n",
    "rev_strand = intronless_dna.reverse_complement()\n",
    "print()\n",
    "print(\"rev_strand:\", rev_strand)\n",
    "print(\"alternatively:\", compl_seq[::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cf7f720",
   "metadata": {},
   "source": [
    "#### Transcription and reverse transcription of sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0d56944",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's now transcribe this piece of DNA into RNA. This is an intronless DNA sequence, so\n",
    "\n",
    "rna_seq = intronless_dna.transcribe()\n",
    "print(\"RNA:\", rna_seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1829dea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also reverse transcribe an RNA sequence to cDNA...\n",
    "cdna_seq = rna_seq.back_transcribe()\n",
    "cdna_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25b9ed49",
   "metadata": {},
   "source": [
    "#### Translation of sequences"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "882657e0",
   "metadata": {},
   "source": [
    "Up to here, we have described an intronless DNA sequence, and 'transcribed' it into an RNA sequence, `rna_seq`. If we know where the CDS starts, we can also translate the CDS into a protein sequence. In this example, the start codon can be found at 53. position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78abb8de",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cds_seq = rna_seq[53:]\n",
    "cds_seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fbbadee",
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_seq = cds_seq.translate()\n",
    "protein_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd4bcf43",
   "metadata": {},
   "source": [
    "The `.translate()` method can take extra arguments to finetune the translation specifics:\n",
    "\n",
    "- **table** - Which codon table to use? Tables are based on [NCBI's Genetic Codes](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi). Default is the \"Standard\" table.\n",
    "\n",
    "\n",
    "- **stop_symbol** - Single character string, what to use for terminators.  Default is the asterisk, *\n",
    "\n",
    "\n",
    "- **to_stop** - Boolean, defaults to False meaning do a full                                                        translation continuing on past any stop codons (translated as the                                                  specified stop_symbol).  If True, translation is terminated at                                                      the first in frame stop codon (and the stop_symbol is not                                                          appended to the returned protein sequence).\n",
    "\n",
    "\n",
    "- **cds** - Boolean, indicates this is a complete CDS.  If True,                                                    this checks the sequence starts with a valid alternative start                                                      codon (which will be translated as methionine, M), that the                                                        sequence length is a multiple of three, and that there is a                                                        single in frame stop codon at the end (this will be excluded                                                        from the protein sequence, regardless of the to_stop option).                                                      If these tests fail, an exception is raised.\n",
    "\n",
    "\n",
    "- **gap** - Single character string to denote symbol used for gaps. Defaults to the minus sign."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41751698",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Codon tables can be defined by string names (NCBI genetic code table name)\n",
    "# or integers (NCBI table ids)\n",
    "codon_tables = [\n",
    "    \"Vertebrate Mitochondrial\", # 2. The Vertebrate Mitochondrial Code\n",
    "    6                           # 6. The Ciliate, Dasycladacean and Hexamita Nuclear Code\n",
    "]\n",
    "\n",
    "# to_stop is boolean, either True or False\n",
    "to_stops = [True, False]\n",
    "\n",
    "protein = cds_seq.translate(table=\"Vertebrate Mitochondrial\",to_stop=True)\n",
    "print(protein)\n",
    "protein = cds_seq.translate(table=\"Vertebrate Mitochondrial\",to_stop=False)\n",
    "print(protein)\n",
    "protein = cds_seq.translate(table=6,to_stop=True)\n",
    "print(protein)\n",
    "protein = cds_seq.translate(table=6,to_stop=False)\n",
    "print(protein)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e81a9fed",
   "metadata": {},
   "source": [
    "A potential source of error may come from the fact that two `Seq` objects may represent different types of biological sequences (protein and cDNA for instance).\n",
    "\n",
    "Indeed, nothing prevents us from erroneously applying `Seq` methods to types of sequences that would 'biologically' not support these actions. For example, if we try to reverse-complement a protein sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f67220e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_seq.reverse_complement() # returns a meaningless sequence by trying to interpret AA as IUPAC nts."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13585854",
   "metadata": {},
   "source": [
    "Or, we can concatenate two incompatible `Seq` objects together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09e909ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "print( cds_seq + protein_seq )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cc5a893",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "## Aligning sequences <a id='16'></a>\n",
    "    \n",
    "One may use any external software to perform multiple sequence alignment, \n",
    "but sometimes it may be more convenient to do it directly within tour code.\n",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4594e97",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import Align\n",
    "\n",
    "aligner = Align.PairwiseAligner()\n",
    "\n",
    "aln = aligner.align( \"ATCGGCCGAT\" , \"TCGGACTAAGGCAGAT\" )[0]\n",
    "print(aln.score)\n",
    "print(aln)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b8965fc",
   "metadata": {},
   "source": [
    "The default arguments of the aligner can be consulted and changed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b41b2f77",
   "metadata": {},
   "outputs": [],
   "source": [
    "print( \"mode:\" , aligner.mode )\n",
    "print( \"match_score:\" , aligner.match_score )\n",
    "print( \"mismatch_score:\" , aligner.mismatch_score )\n",
    "print( \"open_gap_score:\" , aligner.open_gap_score )\n",
    "print( \"extend_gap_score:\" , aligner.extend_gap_score )\n",
    "print( \"end_gap_score:\" , aligner.end_gap_score )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc388c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "aligner.mode = 'global'\n",
    "aligner.match_score = 2\n",
    "aligner.open_gap_score = -1\n",
    "\n",
    "aln = aligner.align( \"ATCGGCCGAT\" , \"TCGGACTAAGGCAGAT\" )[0]\n",
    "print(aln.score)\n",
    "print(aln)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dcea81c",
   "metadata": {},
   "source": [
    "Now for example, we could compute the alignment score of each of our records with another reference sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e57d15a",
   "metadata": {},
   "outputs": [],
   "source": [
    "referenceSeq = 'MSNFQLSSVSFSSPILPSAVDDNSDTKLDVIRNTVNFSASIWGDQFLTFEEPDDLEMEKKVVEVLKEEVRKELVIKGSSNESLQHMKLIELIDVVQRLGVAYHFEEEIEEALKHIYVTYGEKWVDLNNLHNLSLWFRLLRQQGFNVSSGIFKYHMYEKGNFKESLCEDAQGLLALYEASYMRVEGEKVLDDALEFTKTHLAIIAKDPSCDSLLRAQIQEALKQPLRKRLPRLEAVRYIPIYQQDVSHSEVLLKLAKSDFNVLQSMHKEELSQICKWWKDLDMQKKLPYVRDRLIEGYFWILGIYFEPQHSHTRIFLMKTSMWLIVLDDTYDNYGTYEELKIFTDAVQRWSMSCLDLLPEYMKLIYQELLNHHQEMEESLEKEGKTYQIHYVIEMAKEVLENNLVEAKWLKEGYMPTLEEYMSVSMKTCTYGLMIARSFVGRVDNMVTEDTFKWVATYPPIVKAACLVLRLMDDITTHKEEQERGHVASSIECYQKETGASEKEACEFISNMVEDAWKVINRESLRPTDIPFSLLPSTINFARACDVIYKVNDSYTHARKEMINHIKSLLVHPLAI'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7e261ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "alignments = {}\n",
    "\n",
    "aligner = Align.PairwiseAligner()\n",
    "\n",
    "for rec in  records :\n",
    "    alignments[ rec.id ] = aligner.align( referenceSeq , rec.seq ).score\n",
    "    \n",
    "alignments"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bb5ba02",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "### Now, you can try to solve the exercises. <a id='18'></a>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fc58552",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "# Additional Theory <a id='19'></a>\n",
    "-----------------------------\n",
    "\n",
    "The Biopython module is 'huge'! Under this section, we have included some more information on basic functionality which is not needed for the exercises."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a16c7c4e",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "### On `SeqRecord` iterators and processing large sequence files with `SeqIO` <a id='20'></a>\n",
    "\n",
    "We have learnt how we can iterate over the sequence records contained in an input file using the `.parse()` method of `SeqIO`. There are some important limitations that come with iterators. Below, we will explore some of there limitations, how we can overcome them but at a cost in memory usage. Finally, some alternatives if we need to process large files without iterators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4812a313",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "990270b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A word about iterators. Once we iterate over all values an iterator gets empty.\n",
    "# That means, we can not access the contents any longer.\n",
    "\n",
    "# Let's create an iterator of sequence records:\n",
    "records = SeqIO.parse(\"data/example.fa\", \"fasta\")\n",
    "\n",
    "# We can access the next item (in this case, the first one) via core next() method:\n",
    "print(\"First record is\", next(records).id)\n",
    "\n",
    "# Now, let's iterate over all values with a for loop:\n",
    "for record in records:\n",
    "    print(\"first use\", record.seq[:5])\n",
    "print(\"We have reached the end of the iterator\")\n",
    "\n",
    "# At this point, we have reached the end of the iterator. We can not use it again:\n",
    "for record in records:\n",
    "    print(\"second use\", record.seq[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a87c4f6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# If we force the iterator with next(), it will not return a None, but a StopIteration exception\n",
    "next(records)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edac1a74",
   "metadata": {},
   "source": [
    "Therefore, sometimes it is useful to convert the iterator into a list or a dictionary, given that it is not too big to fit into the memory. Iterators can be easily converted into lists, simply casting a `list()` function. For conversion into a dictionary, the `Bio.SeqIO` provides a specialized method called `.to_dict()`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2b904bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Converting the sequence record iterator into a list\n",
    "\n",
    "records = list(SeqIO.parse(\"data/example.fa\", \"fasta\"))\n",
    "# Now, let's loop over all values with for:\n",
    "for record in records:\n",
    "    print(\"first use\", record.seq[:5])\n",
    "# And again:\n",
    "for record in records:\n",
    "    print(\"second use\", record.seq[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a48aafa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Converting the sequence record iterator into a dictionary\n",
    "\n",
    "records = SeqIO.parse(\"data/example.fa\", \"fasta\")\n",
    "rec_dict = SeqIO.to_dict(records)\n",
    "print(\"Listing all items:\")\n",
    "for rec_id, rec in rec_dict.items():\n",
    "    print(rec_id, rec.seq[:5])\n",
    "print()\n",
    "\n",
    "specific_rec_id = \"ADY55933.1\"\n",
    "print(\"Accessing to a specific record:\", specific_rec_id)\n",
    "print(specific_rec_id, rec_dict[specific_rec_id].seq[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc649a59",
   "metadata": {},
   "source": [
    "For larger files, it isn’t possible to hold everything in memory, so `Bio.SeqIO.to_dict()` is not suitable. In these situations we can use the `Bio.SeqIO.index()` function. This will not populate a dictionary, rather index the input file such that we can access records in a an arbitrary order. This should be used with very large files; for small files it will be slower than the `.to_dict()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ee64ad3",
   "metadata": {},
   "outputs": [],
   "source": [
    "rec_dict = SeqIO.index(\"data/example.fa\", \"fasta\")\n",
    "print(specific_rec_id, rec_dict[specific_rec_id].seq[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "460eef48",
   "metadata": {},
   "source": [
    "\n",
    "[back to the toc](#toc)\n",
    "\n",
    "<br>\n",
    "\n",
    "### More on `SeqFeature` objects <a id='21'></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8f542de",
   "metadata": {},
   "source": [
    "`SeqFeature` object describes a feature along a sequence. Some of its useful attributes are:\n",
    "\n",
    "- `.location` which tells us the 0-based coordinates of the feature\n",
    "- `.type` which holds the type or the annotation of the feature\n",
    "- `.qualifiers` is an ordered dictionary (the keys are always returned in same order) that holds qualifying information about the feature. These are analogous to the qualifiers from a GenBank feature table.\n",
    "\n",
    "`SeqFeature` is defined within its own module: `Bio.SeqFeature`. For more detail on this object and the module, you can start with its [tutorial](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc37) and [documentation](https://biopython.org/docs/1.74/api/Bio.SeqFeature.html).\n",
    "\n",
    "Let's extract some useful information for our sequence record from its features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "608572f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "for feature in rec.features:\n",
    "    print(\"From {} to {}, there is a '{}' feature\".format(\n",
    "        feature.location.start, feature.location.end, feature.type))\n",
    "    for key, val in feature.qualifiers.items():\n",
    "        print(\"    {} is {}\".format(key, val))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4289e184",
   "metadata": {},
   "source": [
    "When we slice a `SeqRecord` that contains `SeqFeature`s, we have learnt that only those features that completely fall within the slice will be kept. Then the positions will be offsetted relative to the new sliced sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9083744b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's grab the last record in our example file\n",
    "records = list(SeqIO.parse(\"data/example.gp\", \"genbank\"))\n",
    "last_record = records[-1]\n",
    "\n",
    "# We have seen above that the zinc binding site is from 47 to 178 (0-based)\n",
    "zinc_binding = last_record[47:178]\n",
    "incomplete_overlap = last_record[30:100]\n",
    "print(zinc_binding)\n",
    "print()\n",
    "\n",
    "# Let's have a closer look at the single feature included within the sliced portion\n",
    "feat = zinc_binding.features[0]\n",
    "print(\"The feature is of '{}' type\".format(feat.type))\n",
    "print(\"It is from {} to {}\".format(feat.location.start, feat.location.end))\n",
    "print(\"Its location is {}\".format(feat.location)) # We have the positions of 4 residues that make contact with Zn\n",
    "print(\"Its qualifiers are: {}\".format(feat.qualifiers))\n",
    "print()\n",
    "\n",
    "# On the other hand, the incompletely overlapping part (look at the features)\n",
    "print(incomplete_overlap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e490eda8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Environment (conda_py311)",
   "language": "python",
   "name": "conda_py311"
  },
  "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}