answers_day_5.ipynb

{
 "metadata": {
  "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.8.5"
  },
  "orig_nbformat": 4,
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.8.5 64-bit ('myenv': conda)"
  },
  "interpreter": {
   "hash": "f5c9dfabb21bd2a91b63810df81acbdc6b5e617e45414f0ef050ca96090c868f"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Day 5 Answers"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Google colab users only:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "pip install Biopython"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/single_seq.fasta\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/short_reads.fastq\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/blast_output.xml\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/short_read_blast_results.xml\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/short_fasta_reads.fasta\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/day5_1b_output.xml\n",
    "!wget https://raw.githubusercontent.com/agmcfarland/biopython_workshop/master/day5_2a_output.xml"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Everyone:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 1a"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "source": [
    "from Bio import SeqIO # from the lecture\n",
    "from Bio.Seq import Seq # from the lecture\n",
    "from Bio.Blast.NCBIWWW import qblast # from the lecture\n",
    "from Bio.Blast import NCBIXML # from the lecture\n",
    "import pandas as pd # from the lecture\n",
    "\n",
    "record = SeqIO.read('single_seq.fasta','fasta')\n",
    "\n",
    "record.seq = record.seq.translate()\n",
    "record"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "SeqRecord(seq=Seq('MNQKNPKDTQNFITSKKHVKEILNHTNISKQDNVIEIGSGKGHFTKELVKMSRS...LFH'), id='ERMA_STAAR', name='ERMA_STAAR', description='ERMA_STAAR Q6GKQ0 rRNA adenine N-6-methyltransferase OS=Staphylococcus aureus (strain MRSA252) OX=282458 GN=ermA1 PE=3 SV=1', dbxrefs=[])"
      ]
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 1b"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "source": [
    "blast_results = qblast(program='blastp', database='nr', sequence=record.seq, hitlist_size=4)\n",
    "\n",
    "output_file = open('day5_1b_output.xml', 'w')\n",
    "output_file.write(blast_results.read())\n",
    "output_file.close()\n",
    "blast_results.close()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 1c"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "source": [
    "# blast_record has to be read with NCBIXML each time you call the for loop\n",
    "blast_result_handle = open('day5_1b_output.xml')\n",
    "blast_record = NCBIXML.read(blast_result_handle)\n",
    "\n",
    "for alignment in blast_record.alignments:\n",
    "    print(alignment.title)\n",
    "    print('\\n')"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "gb|ETJ08129.1| rRNA adenine N-6-methyltransferase, partial [Streptococcus parasanguinis DORA_23_24]\n",
      "\n",
      "\n",
      "gb|EXM57106.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus DAR133]\n",
      "\n",
      "\n",
      "ref|WP_001072201.1| MULTISPECIES: 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Bacilli] >sp|P0A0H1.1| RecName: Full=rRNA adenine N-6-methyltransferase; AltName: Full=Erythromycin resistance protein; AltName: Full=Macrolide-lincosamide-streptogramin B resistance protein [Staphylococcus aureus subsp. aureus Mu50] >sp|P0A0H2.1| RecName: Full=rRNA adenine N-6-methyltransferase; AltName: Full=Erythromycin resistance protein; AltName: Full=Macrolide-lincosamide-streptogramin B resistance protein [Staphylococcus aureus subsp. aureus N315] >sp|P0A0H3.1| RecName: Full=rRNA adenine N-6-methyltransferase; AltName: Full=Erythromycin resistance protein; AltName: Full=Macrolide-lincosamide-streptogramin B resistance protein [Staphylococcus aureus] >sp|Q6GKQ0.1| RecName: Full=rRNA adenine N-6-methyltransferase; AltName: Full=Erythromycin resistance protein; AltName: Full=Macrolide-lincosamide-streptogramin B resistance protein [Staphylococcus aureus subsp. aureus MRSA252] >gb|EGO8568860.1| 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Enterococcus faecalis] >gb|EGP5183003.1| 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Enterococcus hirae] >gb|EHS14301.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus subsp. aureus IS-99] >gb|EHS72760.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus subsp. aureus IS-125] >gb|EHS75181.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus subsp. aureus IS-189] >gb|ENK62286.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus M0562] >gb|EUY42749.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus M0406] >gb|MBD4026367.1| 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Staphylococcus aureus] >emb|CBA13539.1| adenine methylase [Staphylococcus rostri] >emb|CBY88979.1| rRNA methylase [Staphylococcus aureus subsp. aureus ST398] >tpg|HAR4218314.1| 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Staphylococcus aureus ADL-227] >tpg|HAR4232287.1| 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Staphylococcus aureus ADL-331] >tpg|HAR4238652.1| 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Staphylococcus aureus ADL-330]\n",
      "\n",
      "\n",
      "gb|ALY18984.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus] >gb|EUQ74542.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus FVRH6001]\n",
      "\n",
      "\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 1d"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "source": [
    "# blast_record has to be read with NCBIXML each time you call the for loop\n",
    "blast_result_handle = open('day5_1b_output.xml')\n",
    "blast_record = NCBIXML.read(blast_result_handle)\n",
    "\n",
    "for alignment in blast_record.alignments:\n",
    "    # use split() on alignment.title\n",
    "    new_alignment_title = alignment.title.split('>')\n",
    "\n",
    "    #select the first item in new_alignment_title, assign it back to new_alignment_title\n",
    "    new_alignment_title = new_alignment_title[0]\n",
    "\n",
    "    # print the new title\n",
    "    print(new_alignment_title)\n",
    "    print('\\n') # for clarity"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "gb|ETJ08129.1| rRNA adenine N-6-methyltransferase, partial [Streptococcus parasanguinis DORA_23_24]\n",
      "\n",
      "\n",
      "gb|EXM57106.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus DAR133]\n",
      "\n",
      "\n",
      "ref|WP_001072201.1| MULTISPECIES: 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Bacilli] \n",
      "\n",
      "\n",
      "gb|ALY18984.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus] \n",
      "\n",
      "\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 1f"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "source": [
    "from Bio.SeqRecord import SeqRecord\n",
    "\n",
    "data_store = []\n",
    "\n",
    "# blast_record has to be read with NCBIXML each time you call the for loop\n",
    "blast_result_handle = open('day5_1b_output.xml')\n",
    "blast_record = NCBIXML.read(blast_result_handle)\n",
    "\n",
    "for alignment in blast_record.alignments:\n",
    "    # use split()\n",
    "    new_alignment_title = alignment.title.split('>')\n",
    "\n",
    "    #select the first item in new_alignment_title, assign it back to new_alignment_title\n",
    "    new_alignment_title = new_alignment_title[0]\n",
    "\n",
    "    # assign hsp.sbjct to sbjct_sequence\n",
    "    for hsp in alignment.hsps:\n",
    "        sbjct_sequence = hsp.sbjct\n",
    "        break\n",
    "    \n",
    "    # Make a Seq object out of sbjct_sequence, Assign it to extracted_sequence \n",
    "    extracted_sequence = Seq(sbjct_sequence)\n",
    "\n",
    "    # make a SeqRecord object called extracted_record. The first item will be the extracted_sequence, id is the new_alignment_title\n",
    "    extracted_record = SeqRecord(extracted_sequence, id = new_alignment_title, description = '')\n",
    "\n",
    "    # store the SeqRecord in a list\n",
    "    data_store.append(extracted_record)\n",
    "data_store"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "[SeqRecord(seq=Seq('MNQKNPKDTQNFITSKKHVKEILNHTNISKQDNVIEIGSGKGHFTKELVKMSRS...LFH'), id='gb|ETJ08129.1| rRNA adenine N-6-methyltransferase, partial [Streptococcus parasanguinis DORA_23_24]', name='<unknown name>', description='', dbxrefs=[]),\n",
       " SeqRecord(seq=Seq('MNQKNPKDTQNFITSKKHVKEILNHTNISKQDNVIEIGSGKGHFTKELVKMSRS...LFH'), id='gb|EXM57106.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus DAR133]', name='<unknown name>', description='', dbxrefs=[]),\n",
       " SeqRecord(seq=Seq('MNQKNPKDTQNFITSKKHVKEILNHTNISKQDNVIEIGSGKGHFTKELVKMSRS...LFH'), id='ref|WP_001072201.1| MULTISPECIES: 23S rRNA (adenine(2058)-N(6))-methyltransferase Erm(A) [Bacilli] ', name='<unknown name>', description='', dbxrefs=[]),\n",
       " SeqRecord(seq=Seq('MNQKNPKDTQNFITSKKHVKEILNHTNISKQDNVIEIGSGKGHFTKELVKMSRS...LFH'), id='gb|ALY18984.1| rRNA adenine N-6-methyltransferase [Staphylococcus aureus] ', name='<unknown name>', description='', dbxrefs=[])]"
      ]
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 1g"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "source": [
    "SeqIO.write(data_store,'day5_1g.fasta','fasta')"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 2a"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "source": [
    "# open short read fasta sequences using open() \n",
    "read_seqs = open('short_fasta_reads.fasta')\n",
    "# place the contents of read_seq into memory with read()\n",
    "read_seqs_memory = read_seqs.read()\n",
    "# supply qblast parameters\n",
    "blast_results = qblast(program='blastn', database='nr', sequence=read_seqs_memory, hitlist_size=3)\n",
    "\n",
    "# write the contents of the blast result to a file called 'day5_2a.xml'\n",
    "output_file = open('day5_2a.xml', 'w')\n",
    "output_file.write(blast_results.read())\n",
    "\n",
    "# close the output file\n",
    "output_file.close()\n",
    "# close the blast_results file\n",
    "blast_results.close()\n",
    "# close the fasta file\n",
    "read_seqs.close()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 2b"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "source": [
    "blast_result_handle = open('day5_2a.xml') # open the blast output file\n",
    "blast_records = NCBIXML.parse(blast_result_handle) # use parse to read a blast result with more than one query\n",
    "\n",
    "for blast_record in blast_records:\n",
    "    print('Read ID: ', blast_record.query)\n",
    "    for alignment in blast_record.alignments:\n",
    "        print('Hit title:', alignment.title)\n",
    "    print('\\n')"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:9154:2416\n",
      "Hit title: gi|982198872|gb|CP013997.1| Pseudomonas monteilii strain USDA-ARS-USMARC-56711, complete genome\n",
      "Hit title: gi|675318909|gb|CP009048.1| Pseudomonas alkylphenolica strain KL28 chromosome, complete genome\n",
      "Hit title: gi|2042583920|gb|CP071007.1| Pseudomonas sp. SORT22 chromosome, complete genome\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:16178:2450\n",
      "Hit title: gi|1594476705|gb|CP037922.1| Hymenobacter radiodurans strain 17J36-26 chromosome, complete genome\n",
      "Hit title: gi|1963973487|emb|LR991037.1| Cosmia trapezina genome assembly, chromosome: 18\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:16201:2387\n",
      "Hit title: gi|982198872|gb|CP013997.1| Pseudomonas monteilii strain USDA-ARS-USMARC-56711, complete genome\n",
      "Hit title: gi|1850326283|gb|CP054128.1| Pseudomonas sp. MPDS chromosome, complete genome\n",
      "Hit title: gi|1373984073|gb|CP024081.1| Pseudomonas sp. s211(2017) chromosome\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:24424:2376\n",
      "Hit title: gi|1524514342|gb|CP027762.1| Pseudomonas sp. LBUM920 chromosome, complete genome\n",
      "Hit title: gi|387159426|gb|CP003041.1| Pseudomonas fluorescens A506, complete genome\n",
      "Hit title: gi|2034361303|gb|CP074078.1| Pseudomonas synxantha strain 27 chromosome, complete genome\n",
      "\n",
      "\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 2c"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "source": [
    "blast_result_handle = open('day5_2a.xml') # open the blast output file\n",
    "blast_records = NCBIXML.parse(blast_result_handle) # use parse to read a blast result with more than one query\n",
    "\n",
    "for blast_record in blast_records:\n",
    "    print('Read ID: ', blast_record.query)\n",
    "    for alignment in blast_record.alignments:\n",
    "        print('Hit title:', alignment.title)\n",
    "        for hsp in alignment.hsps:\n",
    "            print('Expected value: ', hsp.expect)\n",
    "    print('\\n')"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:9154:2416\n",
      "Hit title: gi|982198872|gb|CP013997.1| Pseudomonas monteilii strain USDA-ARS-USMARC-56711, complete genome\n",
      "Expected value:  7.08205e-29\n",
      "Hit title: gi|675318909|gb|CP009048.1| Pseudomonas alkylphenolica strain KL28 chromosome, complete genome\n",
      "Expected value:  1.66701e-11\n",
      "Hit title: gi|2042583920|gb|CP071007.1| Pseudomonas sp. SORT22 chromosome, complete genome\n",
      "Expected value:  2.47406e-09\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:16178:2450\n",
      "Hit title: gi|1594476705|gb|CP037922.1| Hymenobacter radiodurans strain 17J36-26 chromosome, complete genome\n",
      "Expected value:  4.32048\n",
      "Hit title: gi|1963973487|emb|LR991037.1| Cosmia trapezina genome assembly, chromosome: 18\n",
      "Expected value:  4.32048\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:16201:2387\n",
      "Hit title: gi|982198872|gb|CP013997.1| Pseudomonas monteilii strain USDA-ARS-USMARC-56711, complete genome\n",
      "Expected value:  1.62956e-24\n",
      "Hit title: gi|1850326283|gb|CP054128.1| Pseudomonas sp. MPDS chromosome, complete genome\n",
      "Expected value:  7.90608e-16\n",
      "Hit title: gi|1373984073|gb|CP024081.1| Pseudomonas sp. s211(2017) chromosome\n",
      "Expected value:  7.90608e-16\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:24424:2376\n",
      "Hit title: gi|1524514342|gb|CP027762.1| Pseudomonas sp. LBUM920 chromosome, complete genome\n",
      "Expected value:  2.7031e-60\n",
      "Hit title: gi|387159426|gb|CP003041.1| Pseudomonas fluorescens A506, complete genome\n",
      "Expected value:  2.7031e-60\n",
      "Hit title: gi|2034361303|gb|CP074078.1| Pseudomonas synxantha strain 27 chromosome, complete genome\n",
      "Expected value:  1.14939e-58\n",
      "\n",
      "\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 2d"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "source": [
    "blast_result_handle = open('day5_2a.xml') # open the blast output file\n",
    "blast_records = NCBIXML.parse(blast_result_handle) # use parse to read a blast result with more than one query\n",
    "\n",
    "e_value_threshold = 0.0001\n",
    "\n",
    "for blast_record in blast_records:\n",
    "    print('Read ID: ', blast_record.query)\n",
    "    for alignment in blast_record.alignments:\n",
    "        print('Hit title:', alignment.title)\n",
    "        for hsp in alignment.hsps:\n",
    "            if hsp.expect <= e_value_threshold:\n",
    "                print('Expected value: ', hsp.expect)\n",
    "            else:\n",
    "                print('Expected value: did not pass')\n",
    "            break\n",
    "    print('\\n')"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:9154:2416\n",
      "Hit title: gi|982198872|gb|CP013997.1| Pseudomonas monteilii strain USDA-ARS-USMARC-56711, complete genome\n",
      "Expected value:  7.08205e-29\n",
      "Hit title: gi|675318909|gb|CP009048.1| Pseudomonas alkylphenolica strain KL28 chromosome, complete genome\n",
      "Expected value:  1.66701e-11\n",
      "Hit title: gi|2042583920|gb|CP071007.1| Pseudomonas sp. SORT22 chromosome, complete genome\n",
      "Expected value:  2.47406e-09\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:16178:2450\n",
      "Hit title: gi|1594476705|gb|CP037922.1| Hymenobacter radiodurans strain 17J36-26 chromosome, complete genome\n",
      "Expected value: did not pass\n",
      "Hit title: gi|1963973487|emb|LR991037.1| Cosmia trapezina genome assembly, chromosome: 18\n",
      "Expected value: did not pass\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:16201:2387\n",
      "Hit title: gi|982198872|gb|CP013997.1| Pseudomonas monteilii strain USDA-ARS-USMARC-56711, complete genome\n",
      "Expected value:  1.62956e-24\n",
      "Hit title: gi|1850326283|gb|CP054128.1| Pseudomonas sp. MPDS chromosome, complete genome\n",
      "Expected value:  7.90608e-16\n",
      "Hit title: gi|1373984073|gb|CP024081.1| Pseudomonas sp. s211(2017) chromosome\n",
      "Expected value:  7.90608e-16\n",
      "\n",
      "\n",
      "Read ID:  M05164:33:000000000-BGDN6:1:1101:24424:2376\n",
      "Hit title: gi|1524514342|gb|CP027762.1| Pseudomonas sp. LBUM920 chromosome, complete genome\n",
      "Expected value:  2.7031e-60\n",
      "Hit title: gi|387159426|gb|CP003041.1| Pseudomonas fluorescens A506, complete genome\n",
      "Expected value:  2.7031e-60\n",
      "Hit title: gi|2034361303|gb|CP074078.1| Pseudomonas synxantha strain 27 chromosome, complete genome\n",
      "Expected value:  1.14939e-58\n",
      "\n",
      "\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Exercise 2e"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "source": [
    "blast_result_handle = open('day5_2a.xml') # open the blast output file\n",
    "blast_records = NCBIXML.parse(blast_result_handle) # use parse to read a blast result with more than one query\n",
    "\n",
    "store_results = []\n",
    "\n",
    "e_value_threshold = 0.0001\n",
    "\n",
    "for blast_record in blast_records:\n",
    "    # get the read_id\n",
    "    read_id = blast_record.query\n",
    "\n",
    "    for alignment in blast_record.alignments:\n",
    "        # get the name of the hit\n",
    "        hit_title = alignment.title\n",
    "\n",
    "        for hsp in alignment.hsps:\n",
    "            # get the evalue \n",
    "            evalue = hsp.expect\n",
    "\n",
    "            # get a qualitative result for the evalue threshold\n",
    "            if hsp.expect <= e_value_threshold:\n",
    "                evalue_pass = 1\n",
    "            else:\n",
    "                evalue_pass = 0\n",
    "            \n",
    "            # append to list\n",
    "            store_results.append([read_id, hit_title, evalue, evalue_pass])\n",
    "\n",
    "\n",
    "store_results"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "[['M05164:33:000000000-BGDN6:1:1101:9154:2416',\n",
       "  'gi|982198872|gb|CP013997.1| Pseudomonas monteilii strain USDA-ARS-USMARC-56711, complete genome',\n",
       "  7.08205e-29,\n",
       "  1],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:9154:2416',\n",
       "  'gi|675318909|gb|CP009048.1| Pseudomonas alkylphenolica strain KL28 chromosome, complete genome',\n",
       "  1.66701e-11,\n",
       "  1],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:9154:2416',\n",
       "  'gi|2042583920|gb|CP071007.1| Pseudomonas sp. SORT22 chromosome, complete genome',\n",
       "  2.47406e-09,\n",
       "  1],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:16178:2450',\n",
       "  'gi|1594476705|gb|CP037922.1| Hymenobacter radiodurans strain 17J36-26 chromosome, complete genome',\n",
       "  4.32048,\n",
       "  0],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:16178:2450',\n",
       "  'gi|1963973487|emb|LR991037.1| Cosmia trapezina genome assembly, chromosome: 18',\n",
       "  4.32048,\n",
       "  0],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:16201:2387',\n",
       "  'gi|982198872|gb|CP013997.1| Pseudomonas monteilii strain USDA-ARS-USMARC-56711, complete genome',\n",
       "  1.62956e-24,\n",
       "  1],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:16201:2387',\n",
       "  'gi|1850326283|gb|CP054128.1| Pseudomonas sp. MPDS chromosome, complete genome',\n",
       "  7.90608e-16,\n",
       "  1],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:16201:2387',\n",
       "  'gi|1373984073|gb|CP024081.1| Pseudomonas sp. s211(2017) chromosome',\n",
       "  7.90608e-16,\n",
       "  1],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:24424:2376',\n",
       "  'gi|1524514342|gb|CP027762.1| Pseudomonas sp. LBUM920 chromosome, complete genome',\n",
       "  2.7031e-60,\n",
       "  1],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:24424:2376',\n",
       "  'gi|387159426|gb|CP003041.1| Pseudomonas fluorescens A506, complete genome',\n",
       "  2.7031e-60,\n",
       "  1],\n",
       " ['M05164:33:000000000-BGDN6:1:1101:24424:2376',\n",
       "  'gi|2034361303|gb|CP074078.1| Pseudomonas synxantha strain 27 chromosome, complete genome',\n",
       "  1.14939e-58,\n",
       "  1]]"
      ]
     },
     "metadata": {},
     "execution_count": 16
    }
   ],
   "metadata": {}
  }
 ]
}