<jupyter_start><jupyter_text>BLAST an unknown sequence  As stated in the introduction, we have an sequence from *D.yakuba*, but we don't know much about it. First, let's examine the [sequence](./files/yakuba.fa), which is saved in the same directory as this notebook.  We will use Linux's `head` command to to preview the first few line of the file. > Tip: To execute a bash command we can place a `!` in front of the command to launch within this Python Jupyter notebook.<jupyter_code>!head ./files/yakuba.fa<jupyter_output><empty_output><jupyter_text>Starting with BiopythonIn these notebooks, we will be using [Biopython](http://biopython.org/) a set of free software tools for a variety of bioinformatics applications. While this tutorial will not teach Biopython comprehensively, you will learn some useful features and we will refer you to the [Biopython documentation](http://biopython.org/wiki/Documentation) to learn more.   Load Biopython and check versionFirst, let's check that Biopython is installed and check the version.<jupyter_code>import Bio
print("Biopython version is " + Bio.__version__)<jupyter_output><empty_output><jupyter_text>> tip: If you did not have Biopython installed, see their [installation instructions](http://biopython.org/wiki/Download)  Load a fasta file for use in BiopythonIn this step, we want to load the yakuba.fa sequence into a variable that can be used in our blast search. To to this we create a variable called `fasta_file` and use Python's `open()` function to read the file. As shown above, the yakuba file is in a folder called `files` at `./files/yakuba.fa`<jupyter_code># Complete this code by entering the name of your file. The filename and 
# filepath should be in quotes

fasta_file = open().read()
fasta_file = open('./files/yakuba.fa').read()<jupyter_output><empty_output><jupyter_text>We can preview what was read into the fasta file by printing it:<jupyter_code>print(fasta_file)<jupyter_output><empty_output><jupyter_text>Preform a BLAST search using BiopythonAs mentioned in the introduction, BLAST is a tool for similarity searching. This is done by taking your **query** sequence (the sequence you want to find matches for), as well as **search parameters** (some optional adjustments to the way you wish to limit or expand your search) and searching a **database** (a repository of known DNA sequences). First, we will load the appropriate Biopython module for doing a BLAST search over the Internet. The [NCBIWWW module](http://biopython.org/DIST/docs/api/Bio.Blast.NCBIWWW-module.html) has a variety of features we will explore in a moment.<jupyter_code>from Bio.Blast import NCBIWWW<jupyter_output><empty_output><jupyter_text>We will do our first BLAST using this piece of Biopython code. > tip: Since this is a real BLAST search, you will get an 'In [\*]' in the cell below for up to several minutes as the search is executed. Don't proceed in the notebook until the '\*' turns into a number.<jupyter_code>blast_result_1 = NCBIWWW.qblast("blastn", "nt", fasta_file)<jupyter_output><empty_output><jupyter_text>The blast result returned by the NCBIWWW.qblast function is not easy to read as it is an [XML file](https://en.wikipedia.org/wiki/XML). We will use some additional code to examine. First, let's save the blast result as its own file. This<jupyter_code>with open("./files/blast_output.xml", "w") as output_xml:
    output_xml.write(blast_result_1.read())
blast_result_1.close()<jupyter_output><empty_output><jupyter_text>We can preview the first few lines of the `blast_output.xml` file and then go on to extract the information we need.<jupyter_code># Use the `!head` command (using the -n argument to specify the 
# number of lines) to preview the first 50 lines of the blast_output.xml file

### your code here
!head -n 50 ./files/blast_output.xml<jupyter_output><empty_output><jupyter_text>Examining your BLAST resultNext, we will use some additional Biopython tools to view the results of our BLAST search stored in the XML file. We will start by importing Biopython's [SearchIO module](http://biopython.org/DIST/docs/api/Bio.SearchIO-module.html):<jupyter_code>from Bio import SearchIO
# you may get a warning that this feature is exprimental, we can ignore for now.<jupyter_output><empty_output><jupyter_text>Next, we will use SearchIO.read to read in the file to a variable and take a look.<jupyter_code>blast_result_1_xml = SearchIO.read('./files/blast_output.xml', 'blast-xml')
print(blast_result_1_xml)<jupyter_output><empty_output><jupyter_text>First let's interpret what we are seeing in this output:The first three lines are giving us some information about BLAST and our search:- **Program: blastn (2.6.1+)**: This is the BLAST tool we are using (blastn) and the version of the software. - **Query: unknown_yakuba_sequence (11001)**: This is the name and length of our sequence- **Target: nt**: This is the database we are searching called - 'nt' (more on this later). The next section (hits) is useful information on which sequences in the 'nt' database were close matches to our query sequence. - **Column 1 ()**: This is the hit number 1..n- **Column 2 ( HSPs)**: The number "[High-scoring Pairs](https://www.ncbi.nlm.nih.gov/books/NBK62051/); these are the number of places where there was a potentially valid match on a given sequence from the target database - **Column 3 (ID + description)**: [Genbank identifiers](https://www.ncbi.nlm.nih.gov/Sitemap/sequenceIDs.html) and description of the matching sequence In this form we are still looking at a lot of information, so let's look at just a single record. On our list of hits, hit 3 is the first D.melanogaster sequence, so let's look at that one.<jupyter_code>blast_hit = next(SearchIO.parse('./files/blast_output.xml', 'blast-xml'))

print(blast_hit[3])<jupyter_output><empty_output><jupyter_text>We can also view the HSP alignments for the D.melaogaster sequence:<jupyter_code>blast_hit = next(SearchIO.parse('./files/blast_output.xml', 'blast-xml'))

for hsps in blast_hit[3].hsps:
    print(hsps)<jupyter_output><empty_output><jupyter_text>Each HSP record gives us some additional information including the location on our D.yakuba sequence (Query) and the D.melanogaster sequence (Hit). For example, the first HSP matches the coordinates [436867:439261] on chromosome 4 of D.melanogaster.   Critical assessment of BLAST resultsTo understand the meaning of our BLAST results, we have to define a few terms. At this point, you should review the [recommended reading](./files/reading1_nihms519883.pdf) to help you answer the following questions.  **Question**: What is a e-value? **Answer**: e value: Expect value is the number of matches by chance to the provided sequence one can expect in a database of a given size. Lower e values indicate more “significant” or better alignments. **Question**: What is a bitscore **Answer**: is a normalized score expressed in bits that lets you estimate the search space you would have to look through before you would expect to find an score as good as or better than this one by chance.  **Question**: What is the difference between blast record IDs have an 'XM' identifier prefixes and others have NM?  What about parametersNow that we have run BLAST using some defaults, what happens when we want to adjust parameters? Let's examine the options:<jupyter_code>import inspect
print(inspect.signature(NCBIWWW.qblast))<jupyter_output><empty_output><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_file, '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>
{'fraction_non_alphanumeric': 0.08873974645786727, 'fraction_numerical': 0.005859166932992437, 'mean_word_length': 5.056774193548387, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 60, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 17}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '505698', 'license': 'MIT', 'n_tokens_mistral': 2780, 'n_tokens_neox': 2643, 'n_words': 1132, 'path': 'notebooks/Notebook_1_BLAST.ipynb', 'repo_name': 'NCBI-Hackathons/EducationalJupyter', 'size': '9387'}