chp4-5-seq-feature.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Seq annotation objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "my_snp = 4350\n",
    "record = SeqIO.read(\"files/NC_005816.gb\", \"genbank\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "source ['taxon:229193']\n",
      "gene ['GeneID:2767712']\n",
      "CDS ['GI:45478716', 'GeneID:2767712']\n"
     ]
    }
   ],
   "source": [
    "for feature in record.features:\n",
    "    if my_snp in feature:\n",
    "        print(\"%s %s\" % (feature.type, feature.qualifiers.get(\"db_xref\")))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Seq import Seq\n",
    "from Bio.SeqFeature import SeqFeature, FeatureLocation\n",
    "example_parent = Seq(\"ACCGAGACGGCAAAGGCTAGCATAGGTATGAGACTTCCTTCCTGCCAGTGCTGAGGAACTGGGAGCCTAC\")\n",
    "example_feature = SeqFeature(FeatureLocation(5, 18), type=\"gene\", strand=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AGCCTTTGCCGTC\n"
     ]
    }
   ],
   "source": [
    "feature_seq = example_parent[example_feature.location.start:example_feature.location.end].reverse_complement()\n",
    "print(feature_seq)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'GACGGCAAAGGCT'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simple_seq = \"ACCGAGACGGCAAAGGCTAGCATAGGTATGAGACTTCCTTCCTGCCAGTGCTGAGGAACTGGGAGCCTAC\"\n",
    "simple_seq[5:18]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('AGCCTTTGCCGTC')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Seq(simple_seq)[5:18].reverse_complement()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AGCCTTTGCCGTC\n"
     ]
    }
   ],
   "source": [
    "feature_seq = example_feature.extract(example_parent)\n",
    "print(feature_seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slicing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG', IUPACAmbiguousDNA()), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=['Project:58037'])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "record = SeqIO.read(\"files/NC_005816.gb\", \"genbank\")\n",
    "record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9609"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(record)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(9609), strand=1), type='source'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(1954), strand=1), type='repeat_region'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(86), ExactPosition(1109), strand=1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(86), ExactPosition(1109), strand=1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(86), ExactPosition(959), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(BeforePosition(110), ExactPosition(209), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(437), ExactPosition(812), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(1105), ExactPosition(1888), strand=1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(1105), ExactPosition(1888), strand=1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(1108), ExactPosition(1885), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(1366), AfterPosition(1669), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(1432), ExactPosition(1456), strand=1), type='misc_feature'),\n",
       " SeqFeature(CompoundLocation([FeatureLocation(ExactPosition(1435), ExactPosition(1459), strand=1), FeatureLocation(ExactPosition(1618), ExactPosition(1621), strand=1)], 'order'), type='misc_feature', location_operator='order'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(1606), ExactPosition(1624), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(2924), ExactPosition(3119), strand=1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(2924), ExactPosition(3119), strand=1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(2924), ExactPosition(3107), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(3485), ExactPosition(3857), strand=1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(3485), ExactPosition(3857), strand=1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(3497), ExactPosition(3626), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(4342), ExactPosition(4780), strand=1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(4342), ExactPosition(4780), strand=1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(4814), ExactPosition(5888), strand=-1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(4814), ExactPosition(5888), strand=-1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(5909), ExactPosition(5911), strand=1), type='variation'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(5933), ExactPosition(5933), strand=1), type='variation'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(5933), ExactPosition(5933), strand=1), type='variation'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(5947), ExactPosition(5948), strand=1), type='variation'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(6004), ExactPosition(6421), strand=1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(6004), ExactPosition(6421), strand=1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(6524), ExactPosition(6525), strand=1), type='variation'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(6663), ExactPosition(7602), strand=1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(6663), ExactPosition(7602), strand=1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(6663), ExactPosition(7599), strand=1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(7788), ExactPosition(8088), strand=-1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(7788), ExactPosition(8088), strand=-1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(7836), ExactPosition(7995), strand=-1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(8087), ExactPosition(8360), strand=-1), type='gene'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(8087), ExactPosition(8360), strand=-1), type='CDS'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(8090), AfterPosition(8357), strand=-1), type='misc_feature'),\n",
       " SeqFeature(FeatureLocation(ExactPosition(8529), ExactPosition(8529), strand=1), type='variation')]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "record.features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: gene\n",
      "location: [4342:4780](+)\n",
      "qualifiers:\n",
      "    Key: db_xref, Value: ['GeneID:2767712']\n",
      "    Key: gene, Value: ['pim']\n",
      "    Key: locus_tag, Value: ['YP_pPCP05']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(record.features[20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: CDS\n",
      "location: [4342:4780](+)\n",
      "qualifiers:\n",
      "    Key: codon_start, Value: ['1']\n",
      "    Key: db_xref, Value: ['GI:45478716', 'GeneID:2767712']\n",
      "    Key: gene, Value: ['pim']\n",
      "    Key: locus_tag, Value: ['YP_pPCP05']\n",
      "    Key: note, Value: ['similar to many previously sequenced pesticin immunity protein entries of Yersinia pestis plasmid pPCP, e.g. gi| 16082683|,ref|NP_395230.1| (NC_003132) , gi|1200166|emb|CAA90861.1| (Z54145 ) , gi|1488655| emb|CAA63439.1| (X92856) , gi|2996219|gb|AAC62543.1| (AF053945) , and gi|5763814|emb|CAB531 67.1| (AL109969)']\n",
      "    Key: product, Value: ['pesticin immunity protein']\n",
      "    Key: protein_id, Value: ['NP_995571.1']\n",
      "    Key: transl_table, Value: ['11']\n",
      "    Key: translation, Value: ['MGGGMISKLFCLALIFLSSSGLAEKNTYTAKDILQNLELNTFGNSLSHGIYGKQTTFKQTEFTNIKSNTKKHIALINKDNSWMISLKILGIKRDEYTVCFEDFSLIRPPTYVAIHPLLIKKVKSGNFIVVKEIKKSIPGCTVYYH']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(record.features[21])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[SeqFeature(FeatureLocation(ExactPosition(42), ExactPosition(480), strand=1), type='gene'), SeqFeature(FeatureLocation(ExactPosition(42), ExactPosition(480), strand=1), type='CDS')]\n"
     ]
    }
   ],
   "source": [
    "sub_record = record[4300:4800]\n",
    "print(sub_record.features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('ATAAATAGATTATTCCAAATAATTTATTTATGTAAGAACAGGATGGGAGGGGGA...TTA', IUPACAmbiguousDNA()), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sub_record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: gene\n",
      "location: [42:480](+)\n",
      "qualifiers:\n",
      "    Key: db_xref, Value: ['GeneID:2767712']\n",
      "    Key: gene, Value: ['pim']\n",
      "    Key: locus_tag, Value: ['YP_pPCP05']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sub_record.features[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: CDS\n",
      "location: [42:480](+)\n",
      "qualifiers:\n",
      "    Key: codon_start, Value: ['1']\n",
      "    Key: db_xref, Value: ['GI:45478716', 'GeneID:2767712']\n",
      "    Key: gene, Value: ['pim']\n",
      "    Key: locus_tag, Value: ['YP_pPCP05']\n",
      "    Key: note, Value: ['similar to many previously sequenced pesticin immunity protein entries of Yersinia pestis plasmid pPCP, e.g. gi| 16082683|,ref|NP_395230.1| (NC_003132) , gi|1200166|emb|CAA90861.1| (Z54145 ) , gi|1488655| emb|CAA63439.1| (X92856) , gi|2996219|gb|AAC62543.1| (AF053945) , and gi|5763814|emb|CAB531 67.1| (AL109969)']\n",
      "    Key: product, Value: ['pesticin immunity protein']\n",
      "    Key: protein_id, Value: ['NP_995571.1']\n",
      "    Key: transl_table, Value: ['11']\n",
      "    Key: translation, Value: ['MGGGMISKLFCLALIFLSSSGLAEKNTYTAKDILQNLELNTFGNSLSHGIYGKQTTFKQTEFTNIKSNTKKHIALINKDNSWMISLKILGIKRDEYTVCFEDFSLIRPPTYVAIHPLLIKKVKSGNFIVVKEIKKSIPGCTVYYH']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sub_record.features[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Record slicing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "record = SeqIO.read(\"files/NC_005816.gb\", \"genbank\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we're shifting the origin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('GATACGCAGTCATATTTTTTACACAATTCTCTAATCCCGACAAGGTCGTAGGTC...GGA', IUPACAmbiguousDNA()), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shifted = record[2000:] + record[:2000]\n",
    "shifted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "41"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(record.features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "40"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(shifted.features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which feature disappeared"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "flag = 1000\n",
    "if flag in record.features:\n",
    "    print(\"%s %s\" % (feature.type, feature.qualifiers.get(\"db_xref\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That didn't help.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# shifted.features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"repeat_region\" is carried to the end of the sequence and \"source\" is lost!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reverse complementing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NC_005816.1 9609 41 1 13\n"
     ]
    }
   ],
   "source": [
    "record = SeqIO.read(\"files/NC_005816.gb\", \"genbank\")\n",
    "print(\"%s %i %i %i %i\" % (record.id, len(record), len(record.features), \n",
    "                          len(record.dbxrefs), len(record.annotations)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'molecule_type': 'DNA',\n",
       " 'topology': 'circular',\n",
       " 'data_file_division': 'BCT',\n",
       " 'date': '21-JUL-2008',\n",
       " 'accessions': ['NC_005816'],\n",
       " 'sequence_version': 1,\n",
       " 'gi': '45478711',\n",
       " 'keywords': [''],\n",
       " 'source': 'Yersinia pestis biovar Microtus str. 91001',\n",
       " 'organism': 'Yersinia pestis biovar Microtus str. 91001',\n",
       " 'taxonomy': ['Bacteria',\n",
       "  'Proteobacteria',\n",
       "  'Gammaproteobacteria',\n",
       "  'Enterobacteriales',\n",
       "  'Enterobacteriaceae',\n",
       "  'Yersinia'],\n",
       " 'references': [Reference(title='Genetics of metabolic variations between Yersinia pestis biovars and the proposal of a new biovar, microtus', ...),\n",
       "  Reference(title='Complete genome sequence of Yersinia pestis strain 91001, an isolate avirulent to humans', ...),\n",
       "  Reference(title='Direct Submission', ...),\n",
       "  Reference(title='Direct Submission', ...)],\n",
       " 'comment': 'PROVISIONAL REFSEQ: This record has not yet been subject to final\\nNCBI review. The reference sequence was derived from AE017046.\\nCOMPLETENESS: full length.'}"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "record.annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TESTING 9609 41 0 0\n"
     ]
    }
   ],
   "source": [
    "rc = record.reverse_complement(id=\"TESTING\")\n",
    "print(\"%s %i %i %i %i\" % (rc.id, len(rc), len(rc.features), len(rc.dbxrefs), len(rc.annotations)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the sequence of both records, first 10 and last 10 nt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('CAGGGGTCGGGTTCGTTACA', IUPACAmbiguousDNA())"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rc.seq[:10] + rc.seq[-10:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('TGTAACGAACCCGACCCCTG', IUPACAmbiguousDNA())"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "record.seq[:10] + record.seq[-10:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SeqIO - multiple records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gi|2765658|emb|Z78533.1|CIZ78533\n",
      "Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', SingleLetterAlphabet())\n",
      "740\n",
      "gi|2765657|emb|Z78532.1|CCZ78532\n",
      "Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAACAG...GGC', SingleLetterAlphabet())\n",
      "753\n",
      "gi|2765656|emb|Z78531.1|CFZ78531\n",
      "Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAGCAG...TAA', SingleLetterAlphabet())\n",
      "748\n",
      "gi|2765655|emb|Z78530.1|CMZ78530\n",
      "Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAAACAACAT...CAT', SingleLetterAlphabet())\n",
      "744\n",
      "gi|2765654|emb|Z78529.1|CLZ78529\n",
      "Seq('ACGGCGAGCTGCCGAAGGACATTGTTGAGACAGCAGAATATACGATTGAGTGAA...AAA', SingleLetterAlphabet())\n",
      "733\n"
     ]
    }
   ],
   "source": [
    "count=0\n",
    "for seq_record in SeqIO.parse(\"files/ls_orchid.fasta\", \"fasta\"):\n",
    "    if count==5: break\n",
    "    print(seq_record.id)\n",
    "    print(repr(seq_record.seq))\n",
    "    print(len(seq_record))\n",
    "    count += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterator/generator is different than lists!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'generator' object is not subscriptable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-29-2e0fd57b377c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mseq_record2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSeqIO\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"files/ls_orchid.fasta\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"fasta\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mseq_record2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m9\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: 'generator' object is not subscriptable"
     ]
    }
   ],
   "source": [
    "seq_record2 = SeqIO.parse(\"files/ls_orchid.fasta\", \"fasta\")\n",
    "seq_record2[9]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAGCAG...GAG', SingleLetterAlphabet()), id='gi|2765647|emb|Z78522.1|CMZ78522', name='gi|2765647|emb|Z78522.1|CMZ78522', description='gi|2765647|emb|Z78522.1|CMZ78522 C.macranthum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(seq_record2)[10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterator will provide each record in a loop. Or you can use `next()` to get the next record."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Cypripedium irapeanum', 'Cypripedium californicum', 'Cypripedium fasciculatum', 'Cypripedium margaritaceum', 'Cypripedium lichiangense', 'Cypripedium yatabeanum', 'Cypripedium guttatum', 'Cypripedium acaule', 'Cypripedium formosanum', 'Cypripedium himalaicum', 'Cypripedium macranthon', 'Cypripedium calceolus', 'Cypripedium segawai', 'Cypripedium parviflorum var. pubescens', 'Cypripedium reginae', 'Cypripedium flavum', 'Cypripedium passerinum', 'Mexipedium xerophyticum', 'Phragmipedium schlimii', 'Phragmipedium besseae', 'Phragmipedium wallisii', 'Phragmipedium exstaminodium', 'Phragmipedium caricinum', 'Phragmipedium pearcei', 'Phragmipedium longifolium', 'Phragmipedium lindenii', 'Phragmipedium lindleyanum', 'Phragmipedium sargentianum', 'Phragmipedium kaiteurum', 'Phragmipedium czerwiakowianum', 'Phragmipedium boissierianum', 'Phragmipedium caudatum', 'Phragmipedium warszewiczianum', 'Paphiopedilum micranthum', 'Paphiopedilum malipoense', 'Paphiopedilum delenatii', 'Paphiopedilum armeniacum', 'Paphiopedilum emersonii', 'Paphiopedilum niveum', 'Paphiopedilum godefroyae', 'Paphiopedilum bellatulum', 'Paphiopedilum concolor', 'Paphiopedilum fairrieanum', 'Paphiopedilum druryi', 'Paphiopedilum tigrinum', 'Paphiopedilum hirsutissimum', 'Paphiopedilum barbigerum', 'Paphiopedilum henryanum', 'Paphiopedilum charlesworthii', 'Paphiopedilum villosum', 'Paphiopedilum exul', 'Paphiopedilum insigne', 'Paphiopedilum gratrixianum', 'Paphiopedilum primulinum', 'Paphiopedilum victoria', 'Paphiopedilum victoria', 'Paphiopedilum glaucophyllum', 'Paphiopedilum supardii', 'Paphiopedilum kolopakingii', 'Paphiopedilum sanderianum', 'Paphiopedilum lowii', 'Paphiopedilum dianthum', 'Paphiopedilum parishii', 'Paphiopedilum haynaldianum', 'Paphiopedilum adductum', 'Paphiopedilum stonei', 'Paphiopedilum philippinense', 'Paphiopedilum rothschildianum', 'Paphiopedilum glanduliferum', 'Paphiopedilum glanduliferum', 'Paphiopedilum sukhakulii', 'Paphiopedilum wardii', 'Paphiopedilum ciliolare', 'Paphiopedilum dayanum', 'Paphiopedilum hennisianum', 'Paphiopedilum callosum', 'Paphiopedilum tonsum', 'Paphiopedilum javanicum', 'Paphiopedilum fowliei', 'Paphiopedilum schoseri', 'Paphiopedilum bougainvilleanum', 'Paphiopedilum hookerae', 'Paphiopedilum papuanum', 'Paphiopedilum mastersianum', 'Paphiopedilum argus', 'Paphiopedilum venustum', 'Paphiopedilum acmodontum', 'Paphiopedilum urbanianum', 'Paphiopedilum appletonianum', 'Paphiopedilum lawrenceanum', 'Paphiopedilum bullenianum', 'Paphiopedilum superbiens', 'Paphiopedilum purpuratum', 'Paphiopedilum barbatum']\n"
     ]
    }
   ],
   "source": [
    "all_species = []\n",
    "for seq_record in SeqIO.parse(\"files/ls_orchid.gbk\", \"genbank\"):\n",
    "    all_species.append(seq_record.annotations[\"organism\"])\n",
    "print(all_species)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now, list comprehension way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Cypripedium irapeanum', 'Cypripedium californicum', 'Cypripedium fasciculatum', 'Cypripedium margaritaceum', 'Cypripedium lichiangense', 'Cypripedium yatabeanum', 'Cypripedium guttatum', 'Cypripedium acaule', 'Cypripedium formosanum', 'Cypripedium himalaicum', 'Cypripedium macranthon', 'Cypripedium calceolus', 'Cypripedium segawai', 'Cypripedium parviflorum var. pubescens', 'Cypripedium reginae', 'Cypripedium flavum', 'Cypripedium passerinum', 'Mexipedium xerophyticum', 'Phragmipedium schlimii', 'Phragmipedium besseae', 'Phragmipedium wallisii', 'Phragmipedium exstaminodium', 'Phragmipedium caricinum', 'Phragmipedium pearcei', 'Phragmipedium longifolium', 'Phragmipedium lindenii', 'Phragmipedium lindleyanum', 'Phragmipedium sargentianum', 'Phragmipedium kaiteurum', 'Phragmipedium czerwiakowianum', 'Phragmipedium boissierianum', 'Phragmipedium caudatum', 'Phragmipedium warszewiczianum', 'Paphiopedilum micranthum', 'Paphiopedilum malipoense', 'Paphiopedilum delenatii', 'Paphiopedilum armeniacum', 'Paphiopedilum emersonii', 'Paphiopedilum niveum', 'Paphiopedilum godefroyae', 'Paphiopedilum bellatulum', 'Paphiopedilum concolor', 'Paphiopedilum fairrieanum', 'Paphiopedilum druryi', 'Paphiopedilum tigrinum', 'Paphiopedilum hirsutissimum', 'Paphiopedilum barbigerum', 'Paphiopedilum henryanum', 'Paphiopedilum charlesworthii', 'Paphiopedilum villosum', 'Paphiopedilum exul', 'Paphiopedilum insigne', 'Paphiopedilum gratrixianum', 'Paphiopedilum primulinum', 'Paphiopedilum victoria', 'Paphiopedilum victoria', 'Paphiopedilum glaucophyllum', 'Paphiopedilum supardii', 'Paphiopedilum kolopakingii', 'Paphiopedilum sanderianum', 'Paphiopedilum lowii', 'Paphiopedilum dianthum', 'Paphiopedilum parishii', 'Paphiopedilum haynaldianum', 'Paphiopedilum adductum', 'Paphiopedilum stonei', 'Paphiopedilum philippinense', 'Paphiopedilum rothschildianum', 'Paphiopedilum glanduliferum', 'Paphiopedilum glanduliferum', 'Paphiopedilum sukhakulii', 'Paphiopedilum wardii', 'Paphiopedilum ciliolare', 'Paphiopedilum dayanum', 'Paphiopedilum hennisianum', 'Paphiopedilum callosum', 'Paphiopedilum tonsum', 'Paphiopedilum javanicum', 'Paphiopedilum fowliei', 'Paphiopedilum schoseri', 'Paphiopedilum bougainvilleanum', 'Paphiopedilum hookerae', 'Paphiopedilum papuanum', 'Paphiopedilum mastersianum', 'Paphiopedilum argus', 'Paphiopedilum venustum', 'Paphiopedilum acmodontum', 'Paphiopedilum urbanianum', 'Paphiopedilum appletonianum', 'Paphiopedilum lawrenceanum', 'Paphiopedilum bullenianum', 'Paphiopedilum superbiens', 'Paphiopedilum purpuratum', 'Paphiopedilum barbatum']\n"
     ]
    }
   ],
   "source": [
    "all_species = [seq_record.annotations[\"organism\"] for seq_record in \n",
    "SeqIO.parse(\"files/ls_orchid.gbk\", \"genbank\")]\n",
    "print(all_species)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## List comprehension\n",
    "\n",
    "Let's see some examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.39272572286499"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([1/x for x in range(1,1000000)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['s', 's', 'a', '.', 'D', 's', 'S', 'r', '.', 'l', '.']"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "par=\"This is a sentencE. AnD this IS another one. Final sentence.\"\n",
    "\n",
    "sent = par.split() #par.split(\". \")\n",
    "\n",
    "[x[-1] for x in sent]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['t', 's', 's', ' ', 't', ' ', 's', 's', 't', 'c']"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sent2 = \"This is a Title Case Sentence.\"\n",
    "\n",
    "[sent2[x].lower() for x in range(0,len(sent2),3)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4, 2, 1, 5, 4, 8]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[len(x.replace(\".\",\"\")) for x in sent2.split()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[28, 27, 28, 28, 26, 27]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "years_of_birth = [1990, 1991, 1990, 1990, 1992, 1991]\n",
    "ages = [2018 - year for year in years_of_birth]\n",
    "ages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading from compressed file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "67518\n"
     ]
    }
   ],
   "source": [
    "import gzip\n",
    "from Bio import SeqIO\n",
    "# with gzip.open(\"files/ls_orchid.gbk.gz\", \"rt\") as handle:\n",
    "handle = gzip.open(\"files/ls_orchid.gbk.gz\", \"rt\")\n",
    "print(sum([len(r) for r in SeqIO.parse(handle, \"gb\")]))  # used []\n",
    "handle.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "## (len(r) for r in SeqIO.parse(handle, \"gb\")) why is this a generator?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading from FastQ file\n",
    "\n",
    "Let's use iterator approach to count number of records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000\n"
     ]
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "handle = gzip.open(\"files/SRR3579118_tiny_1.fastq.gz\", \"rt\")\n",
    "count=0\n",
    "for i in SeqIO.parse(handle, \"fastq\"):\n",
    "    count += 1\n",
    "print(count)\n",
    "handle.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can be more eloborate, let's use list comprehension or generator approach. Generator approach will take up minimal memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000\n"
     ]
    }
   ],
   "source": [
    "import gzip\n",
    "from Bio import SeqIO\n",
    "handle = gzip.open(\"files/SRR3579118_tiny_1.fastq.gz\", \"rt\")\n",
    "print ( sum([1 for i in SeqIO.parse(handle, \"fastq\")]) )\n",
    "handle.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# alternatives\n",
    "# OR sum(1 for i in SeqIO.parse(handle, \"fastq\"))   : generator, notice no []\n",
    "# OR len([1 for i in SeqIO.parse(handle, \"fastq\")])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's retrieve first 10 records via `list()` so that we can slice a range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "handle = gzip.open(\"files/SRR3579118_tiny_1.fastq.gz\", \"rt\")\n",
    "fastq_rec = list(SeqIO.parse(handle, \"fastq\"))[:10]\n",
    "handle.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NGAGAAAATAACTTTATTTCATTGTGGGGAGCGGGCCGATGTCCAGCCTCA\n",
      "NCCGCTCGCAATCACCCAGATTTCAAGAGCGTGGGTGGCGCCCCGAGAGCC\n",
      "NCCTGAACCACACTTCAACCTTAAGACCACTGGTGGTGTTATTTCAAAGCC\n",
      "NGGCAGGCATACTCATCTTTTTCAGTGGGGGTGAATTCAGTGTAGTACAAG\n",
      "TTTGGATGTGTCTGGAGTCTTGGAAGCTTGACTACACTACGTTCTCCTACA\n",
      "CCCCTCTCCCCCTGCCTCTGTTCGCAGGAGAGCGCGCGCGCGCGCGCGCGC\n",
      "CCACTCCCAATAAATCACAGTCAAAATAAATGAAGAGCTCAAGATGACATC\n",
      "TGACNGATGACGCCCCCGCGCCACGCCGCTCAGCGGGATACGCTTCTTGGN\n",
      "CGCTGGAGTGCAGTGGCTATTCACAGGCGCGATCCCACTACTGATCAGCAC\n",
      "CCCTCGTGGACTTTAGGATGCCCTCAGCTATCACTGCCTCGCTGGGAAAGA\n"
     ]
    }
   ],
   "source": [
    "for n in fastq_rec:\n",
    "    print(n.seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2, 26, 27, 32, 33, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 36, 38, 38, 38, 38, 37, 37, 38, 38, 38, 38, 38, 38, 38]\n",
      "[2, 26, 27, 32, 33, 36, 34, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 35, 37, 38, 38, 38, 38, 38, 38, 38, 35, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 36, 38]\n",
      "[2, 25, 27, 32, 32, 31, 29, 29, 29, 16, 26, 28, 14, 16, 26, 36, 37, 36, 38, 38, 38, 38, 31, 37, 35, 16, 16, 34, 37, 29, 31, 37, 33, 33, 25, 37, 24, 35, 37, 38, 38, 38, 38, 37, 38, 36, 38, 38, 38, 38, 38]\n",
      "[2, 26, 27, 32, 33, 37, 38, 38, 38, 38, 38, 38, 38, 38, 37, 38, 38, 38, 38, 38, 38, 38, 37, 31, 37, 37, 38, 38, 38, 38, 38, 38, 38, 15, 36, 37, 38, 34, 38, 38, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 37]\n",
      "[25, 31, 33, 33, 31, 38, 38, 36, 38, 35, 35, 26, 37, 34, 38, 36, 31, 38, 38, 38, 38, 38, 38, 16, 37, 38, 15, 36, 38, 38, 38, 35, 38, 35, 38, 16, 36, 37, 38, 38, 38, 38, 38, 38, 37, 37, 38, 38, 38, 38, 38]\n",
      "[32, 28, 33, 33, 32, 35, 16, 16, 31, 37, 38, 38, 38, 38, 38, 31, 38, 38, 38, 38, 37, 38, 37, 38, 38, 37, 38, 38, 35, 38, 36, 37, 37, 38, 38, 38, 35, 38, 38, 38, 24, 33, 38, 38, 38, 38, 38, 38, 34, 35, 2]\n",
      "[32, 31, 33, 32, 33, 38, 37, 15, 31, 36, 16, 16, 37, 37, 38, 36, 37, 38, 36, 29, 28, 36, 31, 37, 38, 38, 38, 35, 31, 37, 38, 25, 16, 25, 34, 16, 27, 31, 34, 35, 16, 31, 25, 16, 33, 25, 35, 16, 37, 38, 38]\n",
      "[25, 29, 27, 32, 2, 31, 26, 36, 35, 38, 38, 38, 38, 38, 38, 38, 38, 38, 14, 36, 34, 36, 27, 37, 34, 36, 33, 38, 27, 36, 35, 36, 38, 16, 34, 34, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 37, 38, 2]\n",
      "[18, 25, 30, 32, 33, 34, 38, 33, 38, 38, 38, 38, 38, 36, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 37, 37]\n",
      "[33, 33, 33, 33, 33, 36, 37, 36, 38, 38, 38, 38, 38, 38, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 33, 37, 38, 38, 38, 38, 38, 38, 38, 37, 38, 38, 33, 29, 31, 38, 38, 38]\n"
     ]
    }
   ],
   "source": [
    "for n in fastq_rec:\n",
    "    print(n.letter_annotations['phred_quality'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Average phred quality for first 10 sequences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[36, 36, 31, 35, 35, 34, 30, 33, 36, 36]"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[sum(n.letter_annotations['phred_quality'])//51 for n in fastq_rec]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is average quality per coordinate? (see bottom of the notebook for solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading from net"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "O23729\n",
      "CHS3_BROFI\n",
      "RecName: Full=Chalcone synthase 3; EC=2.3.1.74; AltName: Full=Naringenin-chalcone synthase 3;\n",
      "Seq('MAPAMEEIRQAQRAEGPAAVLAIGTSTPPNALYQADYPDYYFRITKSEHLTELK...GAE', ProteinAlphabet())\n",
      "Length 394\n",
      "['Acyltransferase', 'Flavonoid biosynthesis', 'Transferase']\n"
     ]
    }
   ],
   "source": [
    "from Bio import ExPASy\n",
    "from Bio import SeqIO\n",
    "with ExPASy.get_sprot_raw(\"O23729\") as handle:\n",
    "    seq_record = SeqIO.read(handle, \"swiss\")\n",
    "print(seq_record.id)\n",
    "print(seq_record.name)\n",
    "print(seq_record.description)\n",
    "print(repr(seq_record.seq))\n",
    "print(\"Length %i\" % len(seq_record))\n",
    "print(seq_record.annotations[\"keywords\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequence as dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "orchid_dict = SeqIO.to_dict(SeqIO.parse(\"files/ls_orchid.gbk\", \"genbank\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "94"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(orchid_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P.supardii 5.8S rRNA gene and ITS1 and ITS2 DNA\n"
     ]
    }
   ],
   "source": [
    "seq_record = orchid_dict[\"Z78475.1\"]\n",
    "print(seq_record.description)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACATAATAATTGATCGAGTTAATCTGGAGGATCAGTTTACTTTGGTCACCCATGGGCATCTGCTCTTGCAGTGACCTGGATTTGCCATCGAGCCTCCTTGGGAGCTTTCTTGCTGGCGATCTAAACCCGTCCCGGCGCAGTTTTGCGCCAAGTCATATGACACATAATTGGAAGGGGGTGGCATGCTGCCTTGACCCTCCCCAAATTATTTTTTTGACAACTCTCAGCAACGGATATCTCGGCTCTTNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNATCAGGCCAAGGGCACGCCTGCCTGGGCATTGCGAGTCATATCTCTCCCTTAATGAGGCTGTCCATACATACTGTTCAGCCAATGCGGATGTGAGTTTGGCCCCTTGTTCTTTGGTACGGGGGGTCTAAGAGCTGCATGGGCTTTTGATGGTCCAAAATACGGCAAGAGGTGGACGAACTATGCTACAACAAAATTGTTGTGCGAATGCCCCGGGTTGTCGTATTAGATGGGCCAGCATAATCTAAAGACCCTTTTGAACCCCATTGGAGGCCCATCAACCCATGATCAGTTGACGGCCATTTGGTTGCGACCCAGGTCAGGT\n"
     ]
    }
   ],
   "source": [
    "print(seq_record.seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "first 20 nt of first 10 records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['CGTAACAAGGTTTCCGTAGG',\n",
       " 'CGTAACAAGGTTTCCGTAGG',\n",
       " 'CGTAACAAGGTTTCCGTAGG',\n",
       " 'CGTAACAAGGTTTCCGTAGG',\n",
       " 'ACGGCGAGCTGCCGAAGGAC',\n",
       " 'CGTAACAAGGTTTCCGTAGG',\n",
       " 'CGTAACAAGGTTTCCGTAGG',\n",
       " 'TGTTGAGATAGCAGAATATA',\n",
       " 'CGTAACAAGGTTTCCGTAGG',\n",
       " 'CGTAACCAGGTTTCCGTAGG']"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[str(orchid_dict[x].seq[:20]) for x in orchid_dict.keys()][:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> dictionaries does not allow duplicate keys. `to_dict` will warn if there's a case of duplicate sequence record keys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dictionary keys might be ugly since they get whatever ID is available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['gi|2765658|emb|Z78533.1|CIZ78533',\n",
       " 'gi|2765657|emb|Z78532.1|CCZ78532',\n",
       " 'gi|2765656|emb|Z78531.1|CFZ78531',\n",
       " 'gi|2765655|emb|Z78530.1|CMZ78530',\n",
       " 'gi|2765654|emb|Z78529.1|CLZ78529',\n",
       " 'gi|2765652|emb|Z78527.1|CYZ78527',\n",
       " 'gi|2765651|emb|Z78526.1|CGZ78526',\n",
       " 'gi|2765650|emb|Z78525.1|CAZ78525',\n",
       " 'gi|2765649|emb|Z78524.1|CFZ78524',\n",
       " 'gi|2765648|emb|Z78523.1|CHZ78523']"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orchid_dict = SeqIO.to_dict(SeqIO.parse(\"files/ls_orchid.fasta\", \"fasta\"))\n",
    "list(orchid_dict.keys())[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`to_dict` can parse necessary bits by custom function supplied with `key_function` argument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_accession(record):\n",
    "    \"\"\"\"Given a SeqRecord, return the accession number as a string.\n",
    "    e.g. \"gi|2765613|emb|Z78488.1|PTZ78488\" -> \"Z78488.1\"\n",
    "    \"\"\"\n",
    "    parts = record.id.split(\"|\")\n",
    "    assert len(parts) == 5 and parts[0] == \"gi\" and parts[2] == \"emb\"\n",
    "    return parts[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Z78533.1',\n",
       " 'Z78532.1',\n",
       " 'Z78531.1',\n",
       " 'Z78530.1',\n",
       " 'Z78529.1',\n",
       " 'Z78527.1',\n",
       " 'Z78526.1',\n",
       " 'Z78525.1',\n",
       " 'Z78524.1',\n",
       " 'Z78523.1']"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orchid_dict = SeqIO.to_dict(SeqIO.parse(\"files/ls_orchid.fasta\", \"fasta\"), key_function=get_accession)\n",
    "list(orchid_dict.keys())[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SEGUID checksum as key\n",
    "\n",
    "It's possible to use checksum of the sequence to generate unique keys. SEGUID checksum is one of the checksum functions to be used. So, sequences will be indexed by their nt sequence, not by ID parsed from the record."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('Z78533.1', 'JUEoWn6DPhgZ9nAyowsgtoD9TTo'),\n",
       " ('Z78532.1', 'MN/s0q9zDoCVEEc+k/IFwCNF2pY'),\n",
       " ('Z78531.1', 'xN45pACrTnmBH8a8Y9cWSgoLrwE'),\n",
       " ('Z78530.1', 'yMhI5UUQfFOPcoJXb9B19XUyYlY'),\n",
       " ('Z78529.1', 's1Pnjq9zoSHoI/CG9jQr4GyeMZY'),\n",
       " ('Z78527.1', 'MRf6S1OYhtbdPVS845oCmLTqMgo'),\n",
       " ('Z78526.1', 'QCDzCtL6AKuc+h4UQDD6wFjz3Vs'),\n",
       " ('Z78525.1', 'Zh/FImuuDRmmM/5fXaCYAHS7wo0'),\n",
       " ('Z78524.1', 'hKw0C1fPNpi2KUM6iV0/8IadTX0'),\n",
       " ('Z78523.1', '71sZ82r6eAOBpwZs2solBr/biUs')]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "from Bio.SeqUtils.CheckSum import seguid\n",
    "\n",
    "[(record.id,seguid(record.seq)) for record in SeqIO.parse(\"files/ls_orchid.gbk\", \"genbank\")][:10]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> `seguid` is a generic function, not just for sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'sdV4ERHYT3s/5FoIUuWXWM16h+U'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seguid(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quote from book:\n",
    "\n",
    "> `Bio.SeqIO.to_dict()` function’s `key_function` argument expects a function which turns a SeqRecord into a string. We can’t use the `seguid()` function directly because it expects to be given a Seq object (or a string). However, we can use Python’s `lambda` feature to create a “one off” function to give to `Bio.SeqIO.to_dict()` instead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "from Bio.SeqUtils.CheckSum import seguid\n",
    "seguid_dict = SeqIO.to_dict(SeqIO.parse(\"files/ls_orchid.gbk\", \"genbank\"),\n",
    "                            lambda rec : seguid(rec.seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "JUEoWn6DPhgZ9nAyowsgtoD9TTo CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTGAATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGTGACCCTGATTTGTTGTTGGGCCGCCTCGGGAGCGTCCATGGCGGGTTTGAACCTCTAGCCCGGCGCAGTTTGGGCGCCAAGCCATATGAAAGCATCACCGGCGAATGGCATTGTCTTCCCCAAAACCCGGAGCGGCGGCGTGCTGTCGCGTGCCCAATGAATTTTGATGACTCTCGCAAACGGGAATCTTGGCTCTTTGCATCGGATGGAAGGACGCAGCGAAATGCGATAAGTGGTGTGAATTGCAAGATCCCGTGAACCATCGAGTCTTTTGAACGCAAGTTGCGCCCGAGGCCATCAGGCTAAGGGCACGCCTGCTTGGGCGTCGCGCTTCGTCTCTCTCCTGCCAATGCTTGCCCGGCATACAGCCAGGCCGGCGTGGTGCGGATGTGAAAGATTGGCCCCTTGTGCCTAGGTGCGGCGGGTCCAAGAGCTGGTGTTTTGATGGCCCGGAACCCGGCAAGAGGTGGACGGATGCTGGCAGCAGCTGCCGTGCGAATCCCCCATGTTGTCGTGCTTGTCGGACAGGCAGGAGAACCCTTCCGAACCCCAATGGAGGGCGGTTGACCGCCATTCGGATGTGACCCCAGGTCAGGCGGGGGCACCCGCTGAGTTTACGC\n",
      "MN/s0q9zDoCVEEc+k/IFwCNF2pY CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAACAGAATATATGATCGAGTGAATCTGGAGGACCTGTGGTAACTCAGCTCGTCGTGGCACTGCTTTTGTCGTGACCCTGCTTTGTTGTTGGGCCTCCTCAAGAGCTTTCATGGCAGGTTTGAACTTTAGTACGGTGCAGTTTGCGCCAAGTCATATAAAGCATCACTGATGAATGACATTATTGTCAGAAAAAATCAGAGGGGCAGTATGCTACTGAGCATGCCAGTGAATTTTTATGACTCTCGCAACGGATATCTTGGCTCTAACATCGATGAAGAACGCAGCTAAATGCGATAAGTGGTGTGAATTGCAGAATCCCGTGAACCATCGAGTCTTTGAACGCAAGTTGCGCTCGAGGCCATCAGGCTAAGGGCACGCCTGCCTGGGCGTCGTGTGTTGCGTCTCTCCTACCAATGCTTGCTTGGCATATCGCTAAGCTGGCATTATACGGATGTGAATGATTGGCCCCTTGTGCCTAGGTGCGGTGGGTCTAAGGATTGTTGCTTTGATGGGTAGGAATGTGGCACGAGGTGGAGAATGCTAACAGTCATAAGGCTGCTATTTGAATCCCCCATGTTGTTGTATTTTTTCGAACCTACACAAGAACCTAATTGAACCCCAATGGAGCTAAAATAACCATTGGGCAGTTGATTTCCATTCAGATGCGACCCCAGGTCAGGCGGGGCCACCCGCTGAGTTGAGGC\n"
     ]
    }
   ],
   "source": [
    "# keys are seguid and value is Seq object\n",
    "for record in list(seguid_dict.keys())[:2]:\n",
    "    print (record,seguid_dict[record].seq) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z78532.1\n",
      "C.californicum 5.8S rRNA gene and ITS1 and ITS2 DNA\n"
     ]
    }
   ],
   "source": [
    "record = seguid_dict[\"MN/s0q9zDoCVEEc+k/IFwCNF2pY\"]\n",
    "print(record.id)\n",
    "print(record.description)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converted 94 records\n"
     ]
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "records = SeqIO.parse(\"files/ls_orchid.gbk\", \"genbank\")\n",
    "count = SeqIO.write(records, \"files/my_example.fasta\", \"fasta\")\n",
    "print(\"Converted %i records\" % count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's even shorter way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converted 94 records\n"
     ]
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "count = SeqIO.convert(\"files/ls_orchid.gbk\", \"genbank\", \n",
    "                      \"files/my_example_convert.fasta\", \"fasta\")\n",
    "print(\"Converted %i records\" % count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "74"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "records = (rec for rec in SeqIO.parse(\"files/ls_orchid.gbk\", \"genbank\") \n",
    "           if len(rec) > 700)\n",
    "SeqIO.write(records, \"files/ls_orchid.tab\", \"tab\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two virus genbank files have been downloaded under `files` folder. Let's unzip them and then index them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "! gunzip files/gbvrl*.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 files to index\n"
     ]
    }
   ],
   "source": [
    "import glob\n",
    "from Bio import SeqIO\n",
    "files = glob.glob(\"files/gbvrl*.seq\")\n",
    "print(\"%i files to index\" % len(files))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['files/gbvrl2.seq', 'files/gbvrl1.seq']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indexing will generate an index file from which it can quickly lookup the key and then access to correct place in correct file. It's not required but usually index files are named `.idx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "133667 sequences indexed\n"
     ]
    }
   ],
   "source": [
    "gb_vrl = SeqIO.index_db(\"files/gbvrl.idx\", files, \"genbank\")\n",
    "print(\"%i sequences indexed\" % len(gb_vrl))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indexing might take time when performed for the first time. After it's complete, the lookup is quite fast.\n",
    "\n",
    "> You can use the index as a read only Python dictionary - without having to worry about which file the\n",
    "sequence comes from"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Equine encephalosis virus NS3 gene, complete cds, isolate: Kimron1\n"
     ]
    }
   ],
   "source": [
    "print(gb_vrl[\"AB811634.1\"].description)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indexing compressed files\n",
    "\n",
    "Normally, zip files are not random accessible. `bgzip`, installable from `tabix` package, can prepare random accessible compressed file. There's compressed file under `files` folder, *you don't need to run the following code*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/bin/sh: 1: bgzip: not found\n"
     ]
    }
   ],
   "source": [
    "! bgzip -c files/ls_orchid.gbk > files/ls_orchid.gbk.bgz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's prepare index file for the compressed file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "94"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "orchid_dict = SeqIO.index(\"files/ls_orchid.gbk.bgz\", \"genbank\")\n",
    "len(orchid_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "orchid_dict.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the index is ready to be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "94"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "orchid_dict = SeqIO.index_db(\"files/ls_orchid.gbk.bgz.idx\", \"files/ls_orchid.gbk.bgz\", \"genbank\")\n",
    "len(orchid_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'P.barbatum 5.8S rRNA gene and ITS1 and ITS2 DNA'"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orchid_dict['Z78439.1'].description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "orchid_dict.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "Question: what is average phred quality per coordinate in fastq?\n",
    "\n",
    "Answer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17.3\n",
      "28.0\n",
      "29.7\n",
      "32.3\n",
      "29.5\n",
      "35.4\n",
      "33.1\n",
      "31.5\n",
      "35.4\n",
      "35.2\n",
      "34.3\n",
      "33.6\n",
      "35.4\n",
      "35.1\n",
      "36.6\n",
      "36.7\n",
      "37.1\n",
      "37.8\n",
      "35.4\n",
      "36.9\n",
      "36.5\n",
      "37.6\n",
      "35.3\n",
      "34.5\n",
      "36.9\n",
      "35.4\n",
      "33.0\n",
      "37.1\n",
      "35.8\n",
      "36.8\n",
      "36.8\n",
      "36.0\n",
      "34.9\n",
      "30.9\n",
      "35.6\n",
      "33.0\n",
      "35.0\n",
      "36.3\n",
      "37.5\n",
      "37.7\n",
      "34.3\n",
      "36.8\n",
      "36.5\n",
      "35.6\n",
      "37.4\n",
      "35.9\n",
      "36.8\n",
      "35.1\n",
      "37.4\n",
      "37.4\n",
      "30.6\n"
     ]
    }
   ],
   "source": [
    "from Bio import SeqIO\n",
    "import gzip\n",
    "handle = gzip.open(\"files/SRR3579118_tiny_1.fastq.gz\", \"rt\")\n",
    "fastq_rec = list(SeqIO.parse(handle, \"fastq\"))[:10]\n",
    "handle.close()\n",
    "\n",
    "for i in range(0,51):\n",
    "    print(sum(n.letter_annotations['phred_quality'][i] for n in fastq_rec)/10)\n",
    "    \n",
    "#for coord in range(0,51):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "phred_qual_perbase=[0 for i in range(0,51)]\n",
    "\n",
    "for n in fastq_rec:\n",
    "    phred_temp=n.letter_annotations['phred_quality']\n",
    "    for i in range(0,51):\n",
    "        phred_qual_perbase[i] += phred_temp[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[17.3,\n",
       " 28.0,\n",
       " 29.7,\n",
       " 32.3,\n",
       " 29.5,\n",
       " 35.4,\n",
       " 33.1,\n",
       " 31.5,\n",
       " 35.4,\n",
       " 35.2,\n",
       " 34.3,\n",
       " 33.6,\n",
       " 35.4,\n",
       " 35.1,\n",
       " 36.6,\n",
       " 36.7,\n",
       " 37.1,\n",
       " 37.8,\n",
       " 35.4,\n",
       " 36.9,\n",
       " 36.5,\n",
       " 37.6,\n",
       " 35.3,\n",
       " 34.5,\n",
       " 36.9,\n",
       " 35.4,\n",
       " 33.0,\n",
       " 37.1,\n",
       " 35.8,\n",
       " 36.8,\n",
       " 36.8,\n",
       " 36.0,\n",
       " 34.9,\n",
       " 30.9,\n",
       " 35.6,\n",
       " 33.0,\n",
       " 35.0,\n",
       " 36.3,\n",
       " 37.5,\n",
       " 37.7,\n",
       " 34.3,\n",
       " 36.8,\n",
       " 36.5,\n",
       " 35.6,\n",
       " 37.4,\n",
       " 35.9,\n",
       " 36.8,\n",
       " 35.1,\n",
       " 37.4,\n",
       " 37.4,\n",
       " 0.0]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[phred_qual_perbase[i]/10 for i in range(0,51)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we should do this for all sequences in file, per line"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}