<jupyter_start><jupyter_text>Starting with our mystery fasta file containing an unknown, interesting sequences we found in our studies, let's use Biopython's http BLAST methods to retrieve a BLAST query result.  We can save the result to an XML file the default output method of the BLAST search.<jupyter_code>from Bio.Blast import NCBIWWW
fasta_string = open("./files/yakuba.fa").read()
result_handle = NCBIWWW.qblast("blastn", "nt", fasta_string)

with open("blast_output.xml", "w") as output_xml:
    output_xml.write(result_handle.read())
result_handle.close()<jupyter_output><empty_output><jupyter_text>Now that the result is saved to file, we do not need to re-run this cell, and thus avoid having to re-run the BLAST query.  With the query result saved to file, we can load the result into a BLAST record object (a type of [python class](https://docs.python.org/3/tutorial/classes.html)) and use Biopython's NCBIXML to return the BLAST record object.  The BLAST record object essentially contains all the information you might need in a BLAST record.  You can view the raw [BLAST record code](http://biopython.org/DIST/docs/api/Bio.Blast.Record-module.html) to see what data are stored in the object.  We can view the alignments, alignment scores, expectation values, and other parameters that are stored in the hsps objects within the BLAST record class.  Because there are typically many alignments for each hit from BLAST, we will just examine the alignments from our top hit.<jupyter_code>result_handle = open("blast_output.xml")
from Bio.Blast import NCBIXML
blast_records = NCBIXML.read(result_handle)
for alignment in blast_records.alignments:
    for hsp in alignment.hsps:
        print(hsp)<jupyter_output>Score 4022 (3627 bits), expectation 0.0e+00, alignment length 2011
Query:    3364 CAGCAGAGGACTGACCGAATGAGTCCAATTCCTTTGGTGATAGAT...GAC 5374
               |||||||||||||||||||||||||||||||||||||||||||||...|||
Sbjct:    2910 CAGCAGAGGACTGACCGAATGAGTCCAATTCCTTTGGTGATAGAT...GAC 900
Score 3902 (3519 bits), expectation 0.0e+00, alignment length 1951
Query:     602 AAGAATTTACAGTTCTCTTAATATTGATTACAATAAAGTTGAGAG...GCT 2552
               |||||||||||||||||||||||||||||||||||||||||||||...|||
Sbjct:    5056 AAGAATTTACAGTTCTCTTAATATTGATTACAATAAAGTTGAGAG...GCT 3106
Score 580 (524 bits), expectation 2.5e-143, alignment length 300
Query:    8819 CTTCCTCGTCCTTCACACTACATTTCTTTTTCTTCGATTGCTCGT...TAG 9118
               |||||||||||||||||||||||||||||||||||||||||||||...|||
Sbjct:     690 CTTCCTCGTCCTTCACACTACATTTCTTTTTCTTCGATTGCTCGT...TAG 391
Score 476 (430 bits), expectation 4.3e-115, alignment length 241
Query:    9167 CCGGCGACATAGTAGAAAAGAGTTCGTTTTTAAGGTTCTTTGGAG...CTA 9407
               ||||||||||||||||||[...]<jupyter_text>Let's now look at what the top 10 hits we get from NCBI's Nucleotide (nt) database.  We will view the decriptions.<jupyter_code>result_handle = open("blast_output.xml")
blast_records = NCBIXML.read(result_handle)
counter = 0
for description in blast_records.descriptions:
    print(description)
    counter += 1
    if counter == 10:
        break<jupyter_output>gi|969472224|ref|XM_002099563.2| Drosophila yakuba uncharacterized protein, transcript variant A (Dyak\GE14515), mRNA 4022.0  0.0
gi|969471588|ref|XM_015191338.1| Drosophila yakuba uncharacterized protein, transcript variant B (Dyak\GE14515), mRNA 4022.0  0.0
gi|968004926|ref|XM_001982693.2| Drosophila erecta uncharacterized protein (Dere\GG16448), mRNA 3276.0  0.0
gi|667678241|gb|AE014135.4| Drosophila melanogaster chromosome 4   3276.0  0.0
gi|16874825|gb|AC012390.6| Drosophila melanogaster, chromosome 4, region 101F-102F, BAC clone BACR30L17, complete sequence 3276.0  0.0
gi|1013933709|ref|XM_016180582.1| Drosophila simulans uncharacterized protein (Dsim\GD24381), mRNA 3086.0  0.0
gi|442614379|ref|NM_143665.4| Drosophila melanogaster legless (lgs), mRNA 3062.0  0.0
gi|19550444|gb|AF457205.1| Drosophila melanogaster legless mRNA, complete cds 3062.0  0.0
gi|15291612|gb|AY051651.1| Drosophila melanogaster LD02061 full length cDNA 3062.0  0.0
gi|195354374|ref|XM_002043637.1| Drosophila[...]<jupyter_text>Our top hit is a mRNA transcript record from an organism named <NAME>.  While outside the scope of this tutorial, one could get more information on this record is by using efetch from the nucleotide database of NCBI via Biopython's Entrez module, as shown below.  [EFetch](https://www.ncbi.nlm.nih.gov/books/NBK25499/chapter4.EFetch) is part of [NCBI's E-Utilities Suite](https://www.ncbi.nlm.nih.gov/books/NBK25497/), which provides programmatic access to records and information from NCBI's Entrez databases.<jupyter_code>from Bio import Entrez
Entrez.email = "<EMAIL>"     # Always tell NCBI who you are
top_hit_record = Entrez.efetch(db="nucleotide", id="969472224", rettype="gb", retmode="text")
print(top_hit_record.read())<jupyter_output>LOCUS       XM_002099563            5056 bp    mRNA    linear   INV 22-DEC-2015
DEFINITION  Drosophila yakuba uncharacterized protein, transcript variant A
            (Dyak\GE14515), mRNA.
ACCESSION   XM_002099563
VERSION     XM_002099563.2
DBLINK      BioProject: PRJNA29999
            BioSample: SAMN02952910
KEYWORDS    RefSeq.
SOURCE      Drosophila yakuba
  ORGANISM  Drosophila yakuba
            Eukaryota; Metazoa; Ecdysozoa; Arthropoda; Hexapoda; Insecta;
            Pterygota; Neoptera; Holometabola; Diptera; Brachycera;
            Muscomorpha; Ephydroidea; Drosophilidae; Drosophila; Sophophora.
REFERENCE   1  (bases 1 to 5056)
  AUTHORS   Clark,A.G., Eisen,M.B., Smith,D.R., Bergman,C.M., Oliver,B.,
            Markow,T.A., Kaufman,T.C., Kellis,M., Gelbart,W., Iyer,V.N.,
            Pollard,D.A., Sackton,T.B., Larracuente,A.M., Singh,N.D.,
            Abad,J.P., Abt,D.N., Adryan,B., Aguade,M., Akashi,H.,
            Anderson,W.W., Aquadro,C.F., Ardell,D.H., Arguello,R.,
      [...]<jupyter_text>Returning to our BLAST query, let's explore a little further different parameters for BLAST.  We can first print out the parameters for the qblast function using inspect.signature(function_name).  One could also print the more verbose help(function_name).<jupyter_code>import inspect
print(inspect.signature(NCBIWWW.qblast))<jupyter_output>(program, database, sequence, url_base='https://blast.ncbi.nlm.nih.gov/Blast.cgi', auto_format=None, composition_based_statistics=None, db_genetic_code=None, endpoints=None, entrez_query='(none)', expect=10.0, filter=None, gapcosts=None, genetic_code=None, hitlist_size=50, i_thresh=None, layout=None, lcase_mask=None, matrix_name=None, nucl_penalty=None, nucl_reward=None, other_advanced=None, perc_ident=None, phi_pattern=None, query_file=None, query_believe_defline=None, query_from=None, query_to=None, searchsp_eff=None, service=None, threshold=None, ungapped_alignment=None, word_size=None, alignments=500, alignment_view=None, descriptions=500, entrez_links_new_window=None, expect_low=None, expect_high=None, format_entrez_query=None, format_object=None, format_type='XML', ncbi_gi=None, results_file=None, show_overview=None, megablast=None)<jupyter_text>There are many parameters for this function.  Many of them need not be adjusted in most cases.  For simplicity, let's focus on the four parameters below.  First, we create a python dictionary for a few different parameters we would like to adjust to ultimately compare outputs.  We can start by adjusting the database we want to query against.  We will store the output of each blast XML separately should we choose to access the data at a later point.  As we loop through each blast query for the different programs, we will collect<jupyter_code>%matplotlib inline
import pandas as pd
import matplotlib
blast_params = {'program': 'blastn', 'database': 'nt', 'sequence': fasta_string, 'expect': 10.0}
blast_params['database'] = ['nt', 'refseq_representative_genomes']
print_data = pd.DataFrame()
for database in blast_params['database']:
    db_values = {}
    #result = NCBIWWW.qblast(blast_params['program'], database, blast_params['sequence'], expect=blast_params['expect'])
    file_name = "blast_output_" + database + ".xml"
    #with open(file_name, "w") as output_xml:
    #    output_xml.write(result.read())
    result.close()
    result_input = open(file_name)
    blast_records = NCBIXML.read(result_input)
    for description in blast_records.descriptions:
        if 'score' in db_values:
            db_values['score'].append(description.score)
        else:
            db_values['score'] = [description.score]
        if 'e-value' in db_values:
            db_values['e-value'].append(description.e)
        else:
            db_values['e-value'] = [description.e]
    df = pd.DataFrame.from_dict(db_values)
    df['database'] = database[0:6] # we simply limit the name to the first 6 characters for easier viewing
    frames = [print_data, df]
    print_data = pd.concat(frames, ignore_index=True)
print_data.boxplot('score', by='database')
print_data.boxplot('e-value', by='database')
%matplotlib inline
import pandas as pd
import matplotlib
blast_params = {'program': 'blastn', 'database': 'nt', 'sequence': fasta_string, 'expect': 10.0}
blast_params['database'] = ['nt', 'refseq_representative_genomes']
print_data = pd.DataFrame()
for database in blast_params['database']:
    db_values = {}
    #result = NCBIWWW.qblast(blast_params['program'], database, blast_params['sequence'], expect=blast_params['expect'])
    file_name = "blast_output_" + database + ".xml"
    #with open(file_name, "w") as output_xml:
    #    output_xml.write(result.read())
    #result.close()
    result_input = open(file_name)
    blast_records = NCBIXML.read(result_input)
    for description in blast_records.descriptions:
        if 'score' in db_values:
            db_values['score'].append(description.score)
        else:
            db_values['score'] = [description.score]
        if 'e-value' in db_values:
            db_values['e-value'].append(description.e)
        else:
            db_values['e-value'] = [description.e]
    df = pd.DataFrame.from_dict(db_values)
    df['database'] = database[0:6] # we simply limit the name to the first 6 characters for easier viewing
    frames = [print_data, df]
    print_data = pd.concat(frames, ignore_index=True)
print_data.boxplot('score', by='database')
print_data.boxplot('e-value', by='database')<jupyter_output><empty_output><jupyter_text>Based our initial results, our top hit for a genome in NCBI Entrez has an accession number of AE014135.4.  We can use this entry to explore genes that are present on this model organism to see what potential biological activities our DNA sequence might be involved with.  Before we move on to notebook 2, we need to find out where on this model genomic sequence our test sequence lies (we would not want to look at ALL of the genes on this chromosome, as our sequence only applies to a small portion of that chromosome).<jupyter_code>from Bio import SearchIO
qresults = next(SearchIO.parse('blast_output.xml', 'blast-xml'))
print(qresults[3]) # the 4th target is out first genomic hit<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11207129354577369, 'fraction_numerical': 0.04563867134755604, 'mean_word_length': 4.685772773797339, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 13}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '258762', 'license': 'MIT', 'n_tokens_mistral': 4079, 'n_tokens_neox': 3541, 'n_words': 1030, 'path': 'notebooks/draft_nb1.ipynb', 'repo_name': 'NCBI-Hackathons/EducationalJupyter', 'size': '11129'}