<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Source of the materials
Step1:  Note that the default settings on the NCBI BLAST website are not quite
the same as the defaults on QBLAST. If you get different results, you’ll
need to check the parameters (e.g., the expectation value threshold and
the gap values).
For example, if you have a nucleotide sequence you want to search
against the nucleotide database (nt) using BLASTN, and you know the GI
number of your query sequence, you can use
Step2:  Alternatively, if we have our query sequence already in a FASTA
formatted file, we just need to open the file and read in this record as
a string, and use that as the query argument
Step3:  We could also have read in the FASTA file as a SeqRecord and then
supplied just the sequence itself
Step4:  Supplying just the sequence means that BLAST will assign an identifier
for your sequence automatically. You might prefer to use the SeqRecord
object’s format method to make a FASTA string (which will include the
existing identifier)
Step5:  This approach makes more sense if you have your sequence(s) in a
non-FASTA file format which you can extract using Bio.SeqIO (see
Chapter 5 - Sequence Input and Output.)
Whatever arguments you give the qblast() function, you should get back
your results in a handle object (by default in XML format). The next
step would be to parse the XML output into Python objects representing
the search results (Section [sec
Step6:  After doing this, the results are in the file my_blast.xml and the
original handle has had all its data extracted (so we closed it).
However, the parse function of the BLAST parser (described
in [sec
Step7:  Now that we’ve got the BLAST results back into a handle again, we are
ready to do something with them, so this leads us right into the parsing
section (see Section [sec
Step8:  In this example there shouldn’t be any output from BLASTX to the
terminal, so stdout and stderr should be empty. You may want to check
the output file opuntia.xml has been created.
As you may recall from earlier examples in the tutorial, the
opuntia.fasta contains seven sequences, so the BLAST XML output should
contain multiple results. Therefore use Bio.Blast.NCBIXML.parse() to
parse it as described below in Section [sec
Step9:  If instead you ran BLAST some other way, and have the BLAST output (in
XML format) in the file my_blast.xml, all you need to do is to open
the file for reading
Step10:  Now that we’ve got a handle, we are ready to parse the output. The code
to parse it is really quite small. If you expect a single BLAST result
(i.e., you used a single query)
Step11:  or, if you have lots of results (i.e., multiple query sequences)
Step12:  Just like Bio.SeqIO and Bio.AlignIO (see
Chapters [chapter
Step13:  Or, you can use a for-loop. Note though that you can step through the BLAST records only once. Usually, from each BLAST record you would save the information that you are interested in. If you want to save all returned BLAST records, you can convert the iterator into a list
Step14:  Now you can access each BLAST record in the list with an index as usual. If your BLAST file is huge though, you may run into memory problems trying to save them all in a list.
Usually, you’ll be running one BLAST search at a time. Then, all you need to do is to pick up the first (and only) BLAST record in blast_records
Step15:  I guess by now you’re wondering what is in a BLAST record.
The BLAST record class
A BLAST Record contains everything you might ever want to extract from
the BLAST output. Right now we’ll just show an example of how to get
some info out of the BLAST report, but if you want something in
particular that is not described here, look at the info on the record
class in detail, and take a gander into the code or automatically
generated documentation – the docstrings have lots of good info about
what is stored in each piece of information.
To continue with our example, let’s just print out some summary info
about all hits in our blast report greater than a particular threshold.
The following code does this
<ASSISTANT_TASK:>
Python Code:
from Bio.Blast import NCBIWWW
help(NCBIWWW.qblast)
Explanation: Source of the materials: Biopython cookbook (Adapted)
<font color='red'>
New status: Draft</font>
BLAST
Running BLAST over the Internet
Saving blast output
Running BLAST locally
Parsing BLAST output
The BLAST record class
Parsing plain-text BLAST output
Hey, everybody loves BLAST right? I mean, geez, how can it get any
easier to do comparisons between one of your sequences and every other
sequence in the known world? But, of course, this section isn’t about
how cool BLAST is, since we already know that. It is about the problem
with BLAST – it can be really difficult to deal with the volume of data
generated by large runs, and to automate BLAST runs in general.
Fortunately, the Biopython folks know this only too well, so they’ve
developed lots of tools for dealing with BLAST and making things much
easier. This section details how to use these tools and do useful things
with them.
Dealing with BLAST can be split up into two steps, both of which can be
done from within Biopython. Firstly, running BLAST for your query
sequence(s), and getting some output. Secondly, parsing the BLAST output
in Python for further analysis.
Your first introduction to running BLAST was probably via the NCBI
web-service. In fact, there are lots of ways you can run BLAST, which
can be categorised in several ways. The most important distinction is
running BLAST locally (on your own machine), and running BLAST remotely
(on another machine, typically the NCBI servers). We’re going to start
this chapter by invoking the NCBI online BLAST service from within a
Python script.
NOTE: The following Chapter [chapter:searchio] describes
Bio.SearchIO, an experimental module in Biopython. We intend this to
ultimately replace the older Bio.Blast module, as it provides a more
general framework handling other related sequence searching tools as
well. However, until that is declared stable, for production code please
continue to use the Bio.Blast module for dealing with NCBI BLAST.
Running BLAST over the Internet
We use the function qblast() in the Bio.Blast.NCBIWWW module to call
the online version of BLAST. This has three non-optional arguments:
The first argument is the blast program to use for the search, as a
    lower case string. The options and descriptions of the programs are
    available at
    https://blast.ncbi.nlm.nih.gov/Blast.cgi.. Currently
    qblast only works with blastn, blastp, blastx, tblast and tblastx.
The second argument specifies the databases to search against.
    Again, the options for this are available on the NCBI web pages at
    http://www.ncbi.nlm.nih.gov/BLAST/blast_databases.shtml.
The third argument is a string containing your query sequence. This
    can either be the sequence itself, the sequence in fasta format, or
    an identifier like a GI number.
The qblast function also take a number of other option arguments which
are basically analogous to the different parameters you can set on the
BLAST web page. We’ll just highlight a few of them here:
The argument url_base sets the base URL for running BLAST over the internet. By default it connects to the NCBI, but one can use this to connect to an instance of NCBI BLAST running in the cloud. Please refer to the documentation for the qblast function for further details
The qblast function can return the BLAST results in various
    formats, which you can choose with the optional format_type
    keyword: "HTML", "Text", "ASN.1", or "XML". The default is
    "XML", as that is the format expected by the parser, described in
    section [sec:parsing-blast] below.
The argument expect sets the expectation or e-value threshold.
For more about the optional BLAST arguments, we refer you to the NCBI’s
own documentation, or that built into Biopython:
End of explanation
from Bio.Blast import NCBIWWW
result_handle = NCBIWWW.qblast("blastn", "nt", "8332116")
Explanation: Note that the default settings on the NCBI BLAST website are not quite
the same as the defaults on QBLAST. If you get different results, you’ll
need to check the parameters (e.g., the expectation value threshold and
the gap values).
For example, if you have a nucleotide sequence you want to search
against the nucleotide database (nt) using BLASTN, and you know the GI
number of your query sequence, you can use:
End of explanation
from Bio.Blast import NCBIWWW
fasta_string = open("data/m_cold.fasta").read()
result_handle = NCBIWWW.qblast("blastn", "nt", fasta_string)
Explanation: Alternatively, if we have our query sequence already in a FASTA
formatted file, we just need to open the file and read in this record as
a string, and use that as the query argument:
End of explanation
from Bio.Blast import NCBIWWW
from Bio import SeqIO
record = SeqIO.read("data/m_cold.fasta", format="fasta")
result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)
Explanation: We could also have read in the FASTA file as a SeqRecord and then
supplied just the sequence itself:
End of explanation
from Bio.Blast import NCBIWWW
from Bio import SeqIO
record = SeqIO.read("data/m_cold.fasta", format="fasta")
result_handle = NCBIWWW.qblast("blastn", "nt", record.format("fasta"))
Explanation: Supplying just the sequence means that BLAST will assign an identifier
for your sequence automatically. You might prefer to use the SeqRecord
object’s format method to make a FASTA string (which will include the
existing identifier):
End of explanation
with open("data/my_blast.xml", "w") as out_handle:
    out_handle.write(result_handle.read())
    result_handle.close()
Explanation: This approach makes more sense if you have your sequence(s) in a
non-FASTA file format which you can extract using Bio.SeqIO (see
Chapter 5 - Sequence Input and Output.)
Whatever arguments you give the qblast() function, you should get back
your results in a handle object (by default in XML format). The next
step would be to parse the XML output into Python objects representing
the search results (Section [sec:parsing-blast]), but you might want
to save a local copy of the output file first. I find this especially
useful when debugging my code that extracts info from the BLAST results
(because re-running the online search is slow and wastes the NCBI
computer time).
Saving blast output
We need to be a bit careful since we can use result_handle.read() to
read the BLAST output only once – calling result_handle.read() again
returns an empty string.
End of explanation
 result_handle = open("data/my_blast.xml")
Explanation: After doing this, the results are in the file my_blast.xml and the
original handle has had all its data extracted (so we closed it).
However, the parse function of the BLAST parser (described
in [sec:parsing-blast]) takes a file-handle-like object, so we can
just open the saved file for input:
End of explanation
from Bio.Blast.Applications import NcbiblastxCommandline
help(NcbiblastxCommandline)
blastx_cline = NcbiblastxCommandline(query="opuntia.fasta", db="nr", evalue=0.001,
outfmt=5, out="opuntia.xml")
blastx_cline
print(blastx_cline)
# stdout, stderr = blastx_cline()
Explanation: Now that we’ve got the BLAST results back into a handle again, we are
ready to do something with them, so this leads us right into the parsing
section (see Section [sec:parsing-blast] below). You may want to jump
ahead to that now ….
Running BLAST locally
Introduction
Running BLAST locally (as opposed to over the internet, see
Section [sec:running-www-blast]) has at least major two advantages:
Local BLAST may be faster than BLAST over the internet;
Local BLAST allows you to make your own database to search for
    sequences against.
Dealing with proprietary or unpublished sequence data can be another
reason to run BLAST locally. You may not be allowed to redistribute the
sequences, so submitting them to the NCBI as a BLAST query would not be
an option.
Unfortunately, there are some major drawbacks too – installing all the
bits and getting it setup right takes some effort:
Local BLAST requires command line tools to be installed.
Local BLAST requires (large) BLAST databases to be setup (and
    potentially kept up to date).
To further confuse matters there are several different BLAST packages
available, and there are also other tools which can produce imitation
BLAST output files, such as BLAT.
Standalone NCBI BLAST+
The “new” NCBI
BLAST+
suite was released in 2009. This replaces the old NCBI “legacy” BLAST
package (see below).
This section will show briefly how to use these tools from within
Python. If you have already read or tried the alignment tool examples in
Section [sec:alignment-tools] this should all seem quite
straightforward. First, we construct a command line string (as you would
type in at the command line prompt if running standalone BLAST by hand).
Then we can execute this command from within Python.
For example, taking a FASTA file of gene nucleotide sequences, you might
want to run a BLASTX (translation) search against the non-redundant (NR)
protein database. Assuming you (or your systems administrator) has
downloaded and installed the NR database, you might run:
```
blastx -query opuntia.fasta -db nr -out opuntia.xml -evalue 0.001 -outfmt 5
```
This should run BLASTX against the NR database, using an expectation
cut-off value of $0.001$ and produce XML output to the specified file
(which we can then parse). On my computer this takes about six minutes -
a good reason to save the output to a file so you can repeat any
analysis as needed.
From within Biopython we can use the NCBI BLASTX wrapper from the
Bio.Blast.Applications module to build the command line string, and
run it:
End of explanation
from Bio.Blast import NCBIWWW
result_handle = NCBIWWW.qblast("blastn", "nt", "8332116")
Explanation: In this example there shouldn’t be any output from BLASTX to the
terminal, so stdout and stderr should be empty. You may want to check
the output file opuntia.xml has been created.
As you may recall from earlier examples in the tutorial, the
opuntia.fasta contains seven sequences, so the BLAST XML output should
contain multiple results. Therefore use Bio.Blast.NCBIXML.parse() to
parse it as described below in Section [sec:parsing-blast].
Other versions of BLAST
NCBI BLAST+ (written in C++) was first released in 2009 as a replacement
for the original NCBI “legacy” BLAST (written in C) which is no longer
being updated. There were a lot of changes – the old version had a
single core command line tool blastall which covered multiple
different BLAST search types (which are now separate commands in
BLAST+), and all the command line options were renamed. Biopython’s
wrappers for the NCBI “legacy” BLAST tools have been deprecated and will
be removed in a future release. To try to avoid confusion, we do not
cover calling these old tools from Biopython in this tutorial.
You may also come across Washington University
BLAST (WU-BLAST), and its successor, Advanced
Biocomputing BLAST (AB-BLAST, released in
2009, not free/open source). These packages include the command line
tools wu-blastall and ab-blastall, which mimicked blastall from
the NCBI “legacy” BLAST suite. Biopython does not currently provide
wrappers for calling these tools, but should be able to parse any NCBI
compatible output from them.
Parsing BLAST output
As mentioned above, BLAST can generate output in various formats, such
as XML, HTML, and plain text. Originally, Biopython had parsers for
BLAST plain text and HTML output, as these were the only output formats
offered at the time. Unfortunately, the BLAST output in these formats
kept changing, each time breaking the Biopython parsers. Our HTML BLAST
parser has been removed, but the plain text BLAST parser is still
available (see Section [sec:parsing-blast-deprecated]). Use it at your
own risk, it may or may not work, depending on which BLAST version
you’re using.
As keeping up with changes in BLAST became a hopeless endeavor,
especially with users running different BLAST versions, we now recommend
to parse the output in XML format, which can be generated by recent
versions of BLAST. Not only is the XML output more stable than the plain
text and HTML output, it is also much easier to parse automatically,
making Biopython a whole lot more stable.
You can get BLAST output in XML format in various ways. For the parser,
it doesn’t matter how the output was generated, as long as it is in the
XML format.
You can use Biopython to run BLAST over the internet, as described
    in section [sec:running-www-blast].
You can use Biopython to run BLAST locally, as described
    in section [sec:running-local-blast].
You can do the BLAST search yourself on the NCBI site through your
    web browser, and then save the results. You need to choose XML as
    the format in which to receive the results, and save the final BLAST
    page you get (you know, the one with all of the
    interesting results!) to a file.
You can also run BLAST locally without using Biopython, and save the
    output in a file. Again, you need to choose XML as the format in
    which to receive the results.
The important point is that you do not have to use Biopython scripts to
fetch the data in order to be able to parse it. Doing things in one of
these ways, you then need to get a handle to the results. In Python, a
handle is just a nice general way of describing input to any info source
so that the info can be retrieved using read() and readline()
functions (see Section <span>sec:appendix-handles</span>).
If you followed the code above for interacting with BLAST through a
script, then you already have result_handle, the handle to the BLAST
results. For example, using a GI number to do an online search:
End of explanation
result_handle = open("data/my_blast.xml")
Explanation: If instead you ran BLAST some other way, and have the BLAST output (in
XML format) in the file my_blast.xml, all you need to do is to open
the file for reading:
End of explanation
from Bio.Blast import NCBIXML
blast_record = NCBIXML.read(result_handle)
Explanation: Now that we’ve got a handle, we are ready to parse the output. The code
to parse it is really quite small. If you expect a single BLAST result
(i.e., you used a single query):
End of explanation
from Bio.Blast import NCBIXML
blast_records = NCBIXML.parse(result_handle)
Explanation: or, if you have lots of results (i.e., multiple query sequences):
End of explanation
from Bio.Blast import NCBIXML
blast_records = NCBIXML.parse(result_handle)
blast_record = next(blast_records)
print(blast_record.database_sequences)
# # ... do something with blast_record
Explanation: Just like Bio.SeqIO and Bio.AlignIO (see
Chapters [chapter:Bio.SeqIO] and [chapter:Bio.AlignIO]), we have a
pair of input functions, read and parse, where read is for when
you have exactly one object, and parse is an iterator for when you can
have lots of objects – but instead of getting SeqRecord or
MultipleSeqAlignment objects, we get BLAST record objects.
To be able to handle the situation where the BLAST file may be huge,
containing thousands of results, NCBIXML.parse() returns an iterator.
In plain English, an iterator allows you to step through the BLAST
output, retrieving BLAST records one by one for each BLAST search
result:
End of explanation
for blast_record in blast_records:
#Do something with blast_records
    blast_records = list(blast_records)
blast_records = list(blast_records)
Explanation: Or, you can use a for-loop. Note though that you can step through the BLAST records only once. Usually, from each BLAST record you would save the information that you are interested in. If you want to save all returned BLAST records, you can convert the iterator into a list:
End of explanation
from Bio.Blast import NCBIXML
blast_records = NCBIXML.parse(result_handle)
Explanation: Now you can access each BLAST record in the list with an index as usual. If your BLAST file is huge though, you may run into memory problems trying to save them all in a list.
Usually, you’ll be running one BLAST search at a time. Then, all you need to do is to pick up the first (and only) BLAST record in blast_records:
End of explanation
E_VALUE_THRESH = 0.04
from Bio.Blast import NCBIXML
result_handle = open("data/my_blast.xml", "r")
blast_records = NCBIXML.parse(result_handle)
for alignment in blast_record.alignments:
    for hsp in alignment.hsps:
        if hsp.expect < E_VALUE_THRESH:
            print("****Alignment****")
            print("sequence:", alignment.title)
            print("length:", alignment.length)
            print("e value:", hsp.expect)
            print(hsp.query[0:75] + "...")
            print(hsp.match[0:75] + "...")
            print(hsp.sbjct[0:75] + "...")
Explanation: I guess by now you’re wondering what is in a BLAST record.
The BLAST record class
A BLAST Record contains everything you might ever want to extract from
the BLAST output. Right now we’ll just show an example of how to get
some info out of the BLAST report, but if you want something in
particular that is not described here, look at the info on the record
class in detail, and take a gander into the code or automatically
generated documentation – the docstrings have lots of good info about
what is stored in each piece of information.
To continue with our example, let’s just print out some summary info
about all hits in our blast report greater than a particular threshold.
The following code does this:
End of explanation
<END_TASK>