chp7-blast.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BLAST\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remote blast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIWWW\n",
    "result_handle = NCBIWWW.qblast(\"blastn\", \"nt\", \"8332116\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "remote_result = result_handle.readlines()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1442"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(remote_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = NCBIWWW.qblast(\"blastn\", \"nt\", \"8332116\")\n",
    "blast_records = NCBIXML.parse(result_handle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "blast_rec_list = list(blast_records)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Score 482 (435 bits), expectation 1.2e-117, alignment length 624\n",
      "Query:      59 ACAGAAAATGGGGAGAGAAATGAAGTACTTGGCCATGAAAACTGA...GTA 678\n",
      "               || ||||||||| |||| | |||| ||  |||| |||| | ||||... ||\n",
      "Sbjct:     278 ACCGAAAATGGGCAGAGGAGTGAATTATATGGCAATGACACCTGA...TTA 901\n",
      "Score 468 (423 bits), expectation 7.5e-114, alignment length 590\n",
      "Query:      63 AAAATGGGGAGAGAAATGAAGTACTTGGCCATGAAAACTGATCAA...CCC 649\n",
      "               |||||||| |||  |||| | || ||||| |||||||| || |||...|||\n",
      "Sbjct:      11 AAAATGGGTAGACGAATGGATTATTTGGCGATGAAAACCGAGCAA...CCC 600\n",
      "Score 448 (405 bits), expectation 2.0e-108, alignment length 597\n",
      "Query:      87 TTGGCCATGAAAACTGATCAATTGGCCGTGGCTAATATGATCGAT...TAG 679\n",
      "               ||||||||||||||||| ||| ||||  |||||||| |||| |||...|||\n",
      "Sbjct:      81 TTGGCCATGAAAACTGAGCAAATGGCGTTGGCTAATTTGATAGAT...TAG 677\n",
      "Score 439 (397 bits), expectation 1.0e-105, alignment length 596\n",
      "Query:      63 AAAATGGGGAGAGAAATGAAGTACTTGGCCATGAAAACTGATCAA...CTG 655\n",
      "               |||||||||||    ||| |||  ||||| |||| |||||||| |...|||\n",
      "Sbjct:     170 AAAATGGGGAGG---ATGGAGTTTTTGGCTATGAGAACTGATCCA...CTG 756\n",
      "Score 430 (389 bits), expectation 1.6e-103, alignment length 602\n",
      "Query:      65 AATGGGGAGAGAAATGAAGTACTTGGCCATGAAAACTGATCAATT...CTG 655\n",
      "               ||||||||| |||  | | || |||||||||||||||||||    ... ||\n",
      "Sbjct:     191 AATGGGGAG-GAA--GGATTATTTGGCCATGAAAACTGATCCGGC...ATG 787\n"
     ]
    }
   ],
   "source": [
    "for align in blast_rec_list[0].alignments[0:5]:\n",
    "    for hsp in align.hsps:\n",
    "        print(hsp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parsing BLAST XML file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3056429"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"files/my_blast.xml\")\n",
    "blast_record = NCBIXML.read(result_handle)\n",
    "blast_record.database_sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "refseq_rna\n"
     ]
    }
   ],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"files/my_blast.xml\")\n",
    "blast_records = NCBIXML.parse(result_handle)\n",
    "\n",
    "for rec in blast_records:\n",
    "    print(rec.database)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BLAST record class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****Alignment****\n",
      "sequence: gi|301171311|ref|NR_035856.1| Pan troglodytes microRNA mir-520b (MIR520B), microRNA\n",
      "length: 60\n",
      "e value: 1.71483e-22\n",
      "CCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG...\n",
      "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||...\n",
      "CCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG...\n",
      "****Alignment****\n",
      "sequence: gi|270133242|ref|NR_032573.1| Macaca mulatta microRNA mir-519a (MIR519A), microRNA\n",
      "length: 85\n",
      "e value: 2.54503e-20\n",
      "CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG...\n",
      "||||||| ||||||||||||||||| |||||||||||||||||||||||||||||||||||...\n",
      "CCCTCTAGAGGGAAGCGCTTTCTGTGGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG...\n",
      "****Alignment****\n",
      "sequence: gi|301171322|ref|NR_035857.1| Pan troglodytes microRNA mir-520c (MIR520C), microRNA\n",
      "length: 86\n",
      "e value: 8.88303e-20\n",
      "CCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG...\n",
      "|||||| |||||||| ||||||||||||||||||||||||||||||||||||||||||||...\n",
      "CCTCTAGAGGGAAGCACTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG...\n",
      "****Alignment****\n",
      "sequence: gi|301171322|ref|NR_035857.1| Pan troglodytes microRNA mir-520c (MIR520C), microRNA\n",
      "length: 86\n",
      "e value: 3.31332e-06\n",
      "CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGG...\n",
      "||||||| | |||||| |||||| || | | | | || ||||||||||||| | ||||||...\n",
      "CCCTCTAAAAGGAAGCACTTTCTTTTCTTTCAGACAACAGAAAGTGCTTCCCTCTAGAGG...\n",
      "****Alignment****\n",
      "sequence: gi|301171267|ref|NR_035851.1| Pan troglodytes microRNA mir-519b (MIR519B), microRNA\n",
      "length: 80\n",
      "e value: 8.88303e-20\n",
      "CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGG...\n",
      "||||||| ||||||||||||||||||||||||||||||||||||||| ||||||||||||...\n",
      "CCCTCTAGAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCATCCTTTTAGAGG...\n"
     ]
    }
   ],
   "source": [
    "from Bio.Blast import NCBIXML\n",
    "result_handle = open(\"files/my_blast.xml\")\n",
    "blast_record = NCBIXML.read(result_handle)\n",
    "\n",
    "E_VALUE_THRESH = 0.04\n",
    "for alignment in blast_record.alignments[1:5]:\n",
    "    for hsp in alignment.hsps:\n",
    "        if hsp.expect < E_VALUE_THRESH:\n",
    "            print(\"****Alignment****\")\n",
    "            print(\"sequence:\", alignment.title)\n",
    "            print(\"length:\", alignment.length)\n",
    "            print(\"e value:\", hsp.expect)\n",
    "            print(hsp.query[0:75] + \"...\")\n",
    "            print(hsp.match[0:75] + \"...\")\n",
    "            print(hsp.sbjct[0:75] + \"...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local BLAST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Building a new DB, current time: 12/21/2018 14:05:06\n",
      "New DB name:   /home/alper/biopython/files/orchid_blastdb\n",
      "New DB title:  files/ls_orchid.fasta\n",
      "Sequence type: Nucleotide\n",
      "Keep MBits: T\n",
      "Maximum file size: 1000000000B\n",
      "Adding sequences from FASTA; added 94 sequences in 0.00462008 seconds.\n"
     ]
    }
   ],
   "source": [
    "! makeblastdb -in files/ls_orchid.fasta -dbtype nucl -out files/orchid_blastdb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# BLASTN 2.6.0+\n",
      "# Query: gi|6273291|gb|AF191665.1|AF191665 Opuntia marenae rpl16 gene; chloroplast gene for chloroplast product, partial intron sequence\n",
      "# Database: files/orchid_blastdb\n",
      "# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score\n",
      "# 2 hits found\n",
      "gi|6273291|gb|AF191665.1|AF191665\tgi|2765597|emb|Z78472.1|PLZ78472\t100.000\t13\t0\t0\t139\t151\t208\t220\t2.1\t24.7\n",
      "gi|6273291|gb|AF191665.1|AF191665\tgi|2765594|emb|Z78469.1|PHZ78469\t100.000\t12\t0\t0\t140\t151\t210\t221\t7.3\t22.9\n",
      "# BLASTN 2.6.0+\n",
      "# Query: gi|6273290|gb|AF191664.1|AF191664 Opuntia clavata rpl16 gene; chloroplast gene for chloroplast product, partial intron sequence\n",
      "# Database: files/orchid_blastdb\n",
      "# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score\n",
      "# 2 hits found\n",
      "gi|6273290|gb|AF191664.1|AF191664\tgi|2765597|emb|Z78472.1|PLZ78472\t100.000\t13\t0\t0\t139\t151\t208\t220\t2.1\t24.7\n",
      "gi|6273290|gb|AF191664.1|AF191664\tgi|2765594|emb|Z78469.1|PHZ78469\t100.000\t12\t0\t0\t140\t151\t210\t221\t7.3\t22.9\n",
      "# BLASTN 2.6.0+\n",
      "# Query: gi|6273289|gb|AF191663.1|AF191663 Opuntia bradtiana rpl16 gene; chloroplast gene for chloroplast product, partial intron sequence\n",
      "# Database: files/orchid_blastdb\n",
      "# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score\n",
      "# 2 hits found\n",
      "gi|6273289|gb|AF191663.1|AF191663\tgi|2765597|emb|Z78472.1|PLZ78472\t100.000\t13\t0\t0\t139\t151\t208\t220\t2.1\t24.7\n",
      "gi|6273289|gb|AF191663.1|AF191663\tgi|2765594|emb|Z78469.1|PHZ78469\t100.000\t12\t0\t0\t140\t151\t210\t221\t7.3\t22.9\n",
      "# BLASTN 2.6.0+\n",
      "# Query: gi|6273287|gb|AF191661.1|AF191661 Opuntia kuehnrichiana rpl16 gene; chloroplast gene for chloroplast product, partial intron sequence\n",
      "# Database: files/orchid_blastdb\n",
      "# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score\n",
      "# 2 hits found\n",
      "gi|6273287|gb|AF191661.1|AF191661\tgi|2765597|emb|Z78472.1|PLZ78472\t100.000\t13\t0\t0\t139\t151\t208\t220\t2.1\t24.7\n",
      "gi|6273287|gb|AF191661.1|AF191661\tgi|2765594|emb|Z78469.1|PHZ78469\t100.000\t12\t0\t0\t140\t151\t210\t221\t7.2\t22.9\n",
      "# BLASTN 2.6.0+\n",
      "# Query: gi|6273286|gb|AF191660.1|AF191660 Opuntia echinacea rpl16 gene; chloroplast gene for chloroplast product, partial intron sequence\n",
      "# Database: files/orchid_blastdb\n",
      "# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score\n",
      "# 2 hits found\n",
      "gi|6273286|gb|AF191660.1|AF191660\tgi|2765597|emb|Z78472.1|PLZ78472\t100.000\t13\t0\t0\t139\t151\t208\t220\t2.1\t24.7\n",
      "gi|6273286|gb|AF191660.1|AF191660\tgi|2765594|emb|Z78469.1|PHZ78469\t100.000\t12\t0\t0\t140\t151\t210\t221\t7.2\t22.9\n",
      "# BLASTN 2.6.0+\n",
      "# Query: gi|6273285|gb|AF191659.1|AF191659 Opuntia pachypus rpl16 gene; chloroplast gene for chloroplast product, partial intron sequence\n",
      "# Database: files/orchid_blastdb\n",
      "# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score\n",
      "# 2 hits found\n",
      "gi|6273285|gb|AF191659.1|AF191659\tgi|2765597|emb|Z78472.1|PLZ78472\t100.000\t13\t0\t0\t139\t151\t208\t220\t2.1\t24.7\n",
      "gi|6273285|gb|AF191659.1|AF191659\tgi|2765594|emb|Z78469.1|PHZ78469\t100.000\t12\t0\t0\t140\t151\t210\t221\t7.2\t22.9\n",
      "# BLASTN 2.6.0+\n",
      "# Query: gi|6273284|gb|AF191658.1|AF191658 Opuntia subulata rpl16 gene; chloroplast gene for chloroplast product, partial intron sequence\n",
      "# Database: files/orchid_blastdb\n",
      "# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score\n",
      "# 2 hits found\n",
      "gi|6273284|gb|AF191658.1|AF191658\tgi|2765597|emb|Z78472.1|PLZ78472\t100.000\t13\t0\t0\t139\t151\t208\t220\t2.1\t24.7\n",
      "gi|6273284|gb|AF191658.1|AF191658\tgi|2765594|emb|Z78469.1|PHZ78469\t100.000\t12\t0\t0\t140\t151\t210\t221\t7.2\t22.9\n",
      "# BLAST processed 7 queries\n"
     ]
    }
   ],
   "source": [
    "! blastn -db files/orchid_blastdb -query files/opuntia.fasta -task blastn -outfmt 7   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Blast.Applications import NcbiblastnCommandline\n",
    "blastn_cline = NcbiblastnCommandline(query=\"files/opuntia.fasta\", \n",
    "                                     db=\"files/orchid_blastdb\", \n",
    "                                     evalue=10,\n",
    "                                     outfmt=5, \n",
    "                                     task=\"blastn\",\n",
    "                                     out=\"files/opuntia.xml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "NcbiblastnCommandline(cmd='blastn', out='files/opuntia.xml', outfmt=5, query='files/opuntia.fasta', db='files/orchid_blastdb', evalue=10, task='blastn')"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blastn_cline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blastn -out files/opuntia.xml -outfmt 5 -query files/opuntia.fasta -db files/orchid_blastdb -evalue 10 -task blastn\n"
     ]
    }
   ],
   "source": [
    "print(blastn_cline)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "stdout, stderr = blastn_cline()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stderr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.Blast.Applications import NcbiblastnCommandline\n",
    "blastn_cline = NcbiblastnCommandline(query=\"files/opuntia.fasta\", \n",
    "                                     db=\"files/orchid_blastdb\", \n",
    "                                     evalue=10,\n",
    "                                     outfmt=7, \n",
    "                                     task=\"blastn\",\n",
    "                                     out=\"files/opuntia_blast.out\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "stdout, stderr = blastn_cline()"
   ]
  }
 ],
 "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
}