22111 - User contributions [en]

Exercise:BLAST new

2026-09-22T08:57:59Z

Carol: /* Random protein sequences and BLASTP */

Exercise written by Rasmus Wernersson and modified by Henrik Nielsen and Carolina Barra.

==Introduction==

In this exercise you will learn by doing to BLAST DNA and protein sequences.
BLAST is short for '''B'''asic '''L'''ocal '''A'''lignment '''S'''earch '''T'''ool.

After completion of the exercise you should be able to choose an appropriate BLAST search tool selecting the correct set of parameters, interpret a BLAST result, explain what an E-value means and how database size affects it, and decide what biological conclusion is justified.

For the exercise, we will be using two approaches; the Web interface to '''BLAST''' hosted by the '''NCBI''', and using the '''BioPython library'' in a Notebook using '''Google Colab'''. There are several advantages to each approach:

The Web interface is very intuitive as it provides a graphical summary and links to other databases, while the programmatic way is more flexible allowing us to further customize the search, the output, and do many searches in a single run for a batch of query sequences.

===Useful Links===
* NCBI BLAST main page: http://blast.ncbi.nlm.nih.gov/
* NCBI [http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs BLAST help pages]

 

<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' BLAST is a quite computationally intensive algorithm, and we have in recent years run into issues with overburdening the NCBI server, with 150+ students submitting jobs at the same time from the same DTU IP address. We have therefore implemented a few optimization/workarounds.

:If you still experience long waiting times you can try an alternative internet connection, such as a mobile hotspot. Alternatively, you can copy the pre-saved searches under '''Blast Job ID:XXX''' and use them on the link from BLAST called '''Saved Strategies'''. Please notice that these links will only work for the duration of the exercise session.

[[Image:BLAST_Saved_Strategies.png|center|600px|border]]

 
</div>

==Part 1: Your first BLAST search==

Below is the mRNA sequence for insulin from a close relative, the common macaque (''Macaca fascicularis'').

>J00336.1 Monkey (M. fascicularis) preproinsulin mRNA, complete cds
AACCCTCCGGGACAGGCTGCATCAGAAGAGGTCAGCAAGCAGGTCACTGTCCTTCGGCTATGGCCCTGTG
GATGCGCCTCCTGCCCCTGCTGGCGCTGCTGGCCCTCTGGGGACCTGACCCGGCCCCGGCCTTTGTGAAC
CAGCACCTGTGCGGCTCCCACCTAGTGGAAGCTCTCTACCTGGTGTGCGGGGAACGAGGCTTCTTCTACA
CACCCAAGACCCGCCGGGAGGCAGAGGACCCTCAGGTGGGGCAGGTGGAGCTGGGCGGGGGCCCTGGCGC
AGGCAGCCTGCAGCCCTTGGCGCTGGAGGGGTCCCTGCAGAAGCGCGGCATCGTGGAGCAGTGCTGCACC
AGCATCTGCTCCCTCTACCAGCTGGAGAACTACTGCAACTAA

We will use BLAST to determine whether this sequence looks like the human mRNA for insulin. There are two ways we can do this:

:A. Search the entire database and look for human hits in the results
:B. Specifically search the human part of the database

We will try both possibilities.

=== Search against a broad database ===

Follow the "nucleotide blast" link from the main [https://blast.ncbi.nlm.nih.gov/Blast.cgi BLAST] page.

After pasting the sequence, we are asked to "Choose a search set" or database to perform our search.

In the previous years the students were using a database called "'''nr/nt'''", but today the default is a database called Core Nucleotide database "'''core_nt'''". We want to learn more about these broad databases of nucleotides.
You can click on the question mark of the databases to learn more about them or you can ask your favourite AI as a sparring partner about what is the difference among them.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.1''': What is the main difference of the two databases in regards to content? Why NCBI offers now core_nt as the default. AI loves long and nice explanations that will help you understand, but the TAs have many submissions to review. Please rephrase the answer with your own words in one line.



Based on the answer of Question 1.1, choose the most efficient database from the two.
In the section "Program Selection" select the option "Somewhat similar sequences (blastn)"
Click the BLAST button to launch the search.

After the search has been completed, make yourself familiar with the BLAST output page. Note that you can toggle between hiding and showing each tab by clicking on the title: Descriptions, Graphic Summary, Alignments and Taxonomy (try it!).


First, take a look at the best hit. Since our search sequence (the query) was taken from GenBank which is part of this broad database, we should find an identical sequence in the search. Make sure this is the case!

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.2''':
:Answer the following questions about the best hit:
:* Which database did you use?
:* What is the identifier (Accession)?
:* What is the alignment score ("max score")?
:* What is the percent identity and query coverage?
:* What is the E-value?
:* Are there any gaps in the alignment?

'''Blast Job ID: B2C1HZPV014'''


Then, find the best hit from human (''Homo sapiens'') that is ''not'' a synthetic construct. ('''Tip:''' you can press Ctrl-F in most browsers to search in the page).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.3''':

:Answer the same questions as before about the human hit you found now.


=== Search against Human G+T ===

'''Note:''' In this context, G+T does not mean Gin and Tonic.

Open ''a new window/tab'' with the BLAST home page. Make a new BLASTN search with the same query sequence, this time with Database set to Human genomic + transcript (Human G+T). Remember again to select Somewhat similar sequences (blastn) under Program Selection. Consider the best hit.

'''Note:''' even though you may not have found exactly the same database entry in the two searches, the ''alignment'' should be the same. Make sure this is the case by comparing the actual alignments in the two windows where you made the searches.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.4''':
Answer the same questions as before about the best hit you found in this search.

'''Blast Job ID: B2CCWWRV014'''



===Concerning database size and E-values===

When answering the previous two questions, you may have noticed that the E-value changed, while the alignment score did not. We will now investigate this further.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.5''':
What are the sizes (in basepairs) of the databases we used for the two BLAST searches? ('''Tip:''' Expand the "Search summary" section near the top by clicking it).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.6''':
:*What is the ratio between the database sizes in the two BLAST searches?
:*What is the ratio between the E-values (for the best human hits) in the two BLAST searches?
:*What is the relationship between database size and E-value for hits with identical alignment scores?
:*In conclusion: if the database size is doubled, what will happen to the E-value?

==Part 2: Assessing the statistical significance of BLAST hits==

As discussed in the lecture, there will be a risk of getting false positive results (hits to sequences that are not related to our input sequence) by purely stochastic means. In this first part of the exercise, we will be investigating this further, by examining what happens when we submit randomly generated sequences to BLAST searches.

Rather than giving out a set of pre-generated DNA/Peptide sequences where you only have our word for their randomness, you'll be generating your own random sequences.
It is important to understand that these computer-generated sequences are ''totally random'', just as if you were rolling a die to determine each nucleotide/amino acid in each sequence.

In this part of the exercise, we will run BLAST through a Python Notebook from Google Colab. We will install the package called BioPython that will connect to the NCBI servers where BLAST is hosted.

[[File:Blast exercise in Colab.png|center|600px|border]]

For our purpose, there are several advantages to this approach:

* We don't need to install it on our computers, but we can still perform custom searches which can be fine-grained control of DATA and workflow (everything can be scripted/automated) than running BLAST through a web interface
* NCBI offers direct access to preformatted BLAST databases of all the data that they host:
** GenBank (+ derivates)
** Full Genome database
** Protein databases (Both from translated GenBank and UniProt)

===Links===
A few notes on how to use a [https://colab.research.google.com/notebooks/basic_features_overview.ipynb#scrollTo=WUtu4316QSHL Google Colab]

<div style="background-color: lavender; border: solid thin grey;">
'''IMPORTANT:''' Check that your notebook is always connected to the server.

If your laptop is idle for some time, you may get the error 'runtime disconnected'

[[File:Runtime_disconnected.png|center|300px|border]]
 

To re-connect to the server

[[File:Reconnect_to_a_hosted_runtime.png|center|600px|border]]

Remember that after is disconnected we will need to run all the code cells from that part of the exercise. This could take some minutes for the blast searches, so try to avoid having the screen idle for a long time :)

 
</div>

* Biopython [https://biopython.org/DIST/docs/tutorial/Tutorial.html tutorial]

* Biopython [https://github.com/biopython/biopython/blob/a09f38b468cbedaf25aec6cfed1475b0849011e6/Bio/Blast/NCBIWWW.py examples]
 

==Create your first Python Notebook==

[https://colab.research.google.com/drive/1DpqOH8O_KK0rEohjv-20CCx2THwCq89h#scrollTo=_7tJ_S6wc-Os Click Here]

We are going to make a local copy of this notebook on our own drive.
* First you will need to log in or create a Google account.
* Make a copy of the Colab file in your own drive. Go to File > Save in Drive

[[File:save_a_copy_in_drive.png|center|300px|border]]

Now you can start modifying the file!

* First, we will install BioPython in our notebook. This will only take a few seconds, but note, that this installation is temporary. If your notebook gets idle and disconnected you will need to run this code cell again.

! pip install biopython

* Now import Blast tools from BioPython

from Bio.Blast import NCBIWWW
help(NCBIWWW.qblast)

===Random DNA sequences and BLASTN===

Generate 10 random DNA sequence of length 25bp using Google Colab.

We will now do a BLASTN search using each random sequence as query.

In Google Colab, this will be analogous to the web browser using the option "blastn" in the section "Program Selection". And in this case we will select "Nucleotide Collection (nr/nt)" as the search database. In addition we will have to modify the E-value cut-off ("E_VALUE_THRESH") to '''50'''

'''Note''': For this special situation where we BLAST small artificial sequences, we need to turn off some the automatics NCBI incorporate when short sequences are detected. Otherwise we'll not be able to see the intended results:

In the Google Colab we do this inadvertently setting the parameter "short_query=False". In the web browser version this will be analogous as going into the "Algorithm parameters" section (see the screen shot below) in order to gain access to fine-tuning the options, and '''Deselect''' the "Automatically adjust parameters for short input sequences" option.

[[file:Blastn_cropped+circle.png‎|center|frame|'''Remember to adjust the BLAST settings''']]

Inspect the results.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 2.1''':
:Answer the following small questions, and '''document your findings''' by pasting in examples of alignments / text snippets from the overview plots generated on google colab:

:* What is the typical length of the hits (the alignment length)?
:* What is the typical % identity?
:* In what range is the bit-scores ("max score")?
:** ''Notice: This is conceptually the same as the "alignment score" we have already met in the pairwise alignment exercise''.
:* What is the range of the E-values?

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 2.2''':
:What is the '''biological''' significance of these hits / is there any biological meaning?

===Random protein sequences and BLASTP===
Now it's time to work with a set of '''protein sequences''': Generate 10 peptide sequences of length 25aa using Google colab.

In our Google colab we will now select blastp as the algorithm to use, and the "NR" database (this is the protein version, consisting of translated CDS'es, UniProt etc).

:'''Note 1:''' In our random generated sequences, the distribution of amino acids will be equal (5% prob) and this is different from true biological sequences - however this is not important for this first part of the exercise.
:'''Note 2:''' Please recall from the lecture that the way <tt>BLASTP</tt> selects candidate sequences for full Smith-Waterman alignment is different from <tt>BLASTN</tt>. (<tt>BLASTN</tt> - a single short (11 bp +) perfect match hit is needed. <tt>BLASTP</tt> - a pair of "near match" hits of 3 aa within a 40 aa window is needed). We also need to tweak the parameters this time. We are doing this by "short_query =False", which will be analogous to set in the Web browser version "Algorithm Parameters" section select BLOSUM62 as the alignment matrix to use and set the "Expect threshold" to 1000 (default: 10).

Perform the BLAST search.
Inspect the results.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 2.3''':
:''(Please '''document your answers''' in the same manner as Q2.1)''

:* What is the typical length of the alignment and do they contain gaps?
:* What is the range of E-values?
:* Try to inspect a few of the alignments in details ("+" means similar sequences) - do you find any that look plausible, if we for a moment ignore the length/E-value?
:* If we had used the default E-value cut-off of 10 would any hits have been found?

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 2.4''':
:If we compare the result from BLAST'ing random DNA sequences to random Peptide sequences - which kind of search has the higher risk of returning false positives (results that appear plausible, maybe even significant, but are truly unrelated)? Remember to take E-values into your consideration.

 

==Part 3: using BLAST to transfer functional information by finding homologs==
<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' The following parts can either be done on the google colab or in the webserver. If you choose the web server version remember to limit your search to "bacteria" (taxid: 2) in ALL of this section (PART 3) to make the BLAST searches run quicker. (The organisms we're looking for all belong to the "Bacteria" domain of life, so this restriction is OK).
</div>

===Homo-, Ortho- and Paralogs===

One of the most common ways to use BLAST as a tool is in the situation where you have a sequence of '''unknown function''', and want to '''find out which function it has'''. Since a large amount of sequence data has been gathered over the years, chances are that an '''evolutionarily related''' sequence with known function has already been identified. In general, such a related sequence is known as a "'''homolog'''".

Homo-, Ortho- and Paralogs:
* A '''Homolog''' is a general term that describes a sequence that is related by any evolutionary means.
* An '''Ortholog''' ("Ortho" = True) is a sequence that is "the same gene" in a different organism: The sequences shared a single common ancestor sequence, and have now diverged through speciation (e.g. the Alpha-globin gene in Human and Mouse).
* A '''Paralog''' arises due to a gene duplication within a species. For example, Alpha- and Beta-globin are paralogs.

[[File:Homo_Ortho_Para-log.gif|center|frame|''Image source: [http://www.thegreatgoodplace.com/tt/gwlee/126 gwLee's blog]'' ]]

Notice that in both cases it's possible to transfer information, for example, about gene family / protein domains.
We have already touched upon the comparison of (potentially) evolutionarily related sequences in the pairwise alignment exercise. However, this time we do not start with two sequences we assume are related, but instead, we start with a single sequence ("query sequence") which we will use to search the databases for homologs (we often informally speak of "BLAST hits", when discussing the sequences found).

 



===BLAST example 2===

The sequence is a '''DNA fragment''' from an unknown non-cultivatable microorganism. It was cloned and sequenced directly from DNA extracted from a soil-sample, and it goes by the poetic name "CLONE12". It was amplified using degenerated PCR primers that target the middle ("core cloning") of the sequence of '''a group of known enzymes'''. (I can guarantee this particular sequence is not in the BLAST databases, since I have cloned and sequenced it myself, and it has never been submitted to GenBank).

LOCUS CLONE12.DNA 609 BP DS-DNA UPDATED 06/14/98
DEFINITION UWGCG file capture
ACCESSION -
KEYWORDS -
SOURCE -
COMMENT Non-sequence data from original file:
BASE COUNT 174 A 116 C 162 G 157 T 0 OTHER
ORIGIN ?
clone12.dna Length: 609 Jun 13, 1998 - 03:39 PM Check: 6014 ..
1 AACGGGCACG GGACGCATGT AGCTGGAACA GTGGCAGCCG TAAATAATAA TGGTATCGGA
61 GTTGCCGGGG TTGCAGGAGG AAACGGCTCT ACCAATAGTG GAGCAAGGTT AATGTCCACA
121 CAAATTTTTA ATAGTGATGG GGATTATACA AATAGCGAAA CTCTTGTGTA CAGAGCCATT
181 GTTTATGGTG CAGATAACGG AGCTGTGATC TCGCAAAATA GCTGGGGTAG TCAGTCTCTG
241 ACTATTAAGG AGTTGCAGAA AGCTGCGATC GACTATTTCA TTGATTATGC AGGAATGGAC
301 GAAACAGGAG AAATACAGAC AGGCCCTATG AGGGGAGGTA TATTTATAGC TGCCGCCGGA
361 AACGATAACG TTTCCACTCC AAATATGCCT TCAGCTTATG AACGGGTTTT AGCTGTGGCC
421 TCAATGGGAC CAGATTTTAC TAAGGCAAGC TATAGCACTT TTGGAACATG GACTGATATT
481 ACTGCTCCTG GCGGAGATAT TGACAAATTT GATTTGTCAG AATACGGAGT TCTCAGCACT
541 TATGCCGATA ATTATTATGC TTATGGAGAG GGAACATCCA TGGCTTGTCC ACATGTCGCC
601 GGCGCCGCC
//

[[Image:Office-notes-line_drawing.png|30px|left]][[Image:Cogs_brain.png|right|150px]]
:'''QUESTION 3.1 (Long question - read all)''':
: ''Your task is now to find out '''what kind of enzyme''' this sequence is likely to encode, '''using the methods''' you have learned''.

'''INSTRUCTIONS''': You are free to write the combined answer to this question in a free-style essay-like fashion - just be sure to include the subquestions in your answers. In an exam situation you will need to put all the clues together yourself, reason about the tools/databases to use, and document your findings.

'''STEP 1 - cleaning up the sequence''':

The sequence is (more or less) in GenBank format and the NCBI BLAST server expects the input to be in FASTA format, or to be "raw" unformatted sequence.

* There are two solutions to this:
** Copy the sequence into a text-editor and manually create a FASTA file ("search and replace" and/or "rectangular selection" is useful for the reformatting). This is the most robust solution: it will always work. (Look at [[Plain text files and Geany|the Geany exercise]] for a reminder of how to do this).
** Hope the creators of the web-server you're using were kind enough to automatically remove non-DNA letters (paste in ONLY the DNA lines) - this turns out to be the case for both NCBI BLAST and VirtualRibosome, but it ''cannot be universally relied upon''.

'''Subquestion''': convert the sequence to FASTA format (manually, in Geany) and quote it in your report.

'''STEP 2 - thinking about the task''':

Consider the following before you start on solving this task:
* Based on the information given: is the sequence protein-coding?
* If it is, can you trust it will contain both a START and STOP codon?
* Do we know if the sequence is sense or anti-sense?
and think which consequences the answers to these questions should have for your choice of methods and parameters.

'''Subquestion''': Give a summary of your considerations.

'''STEP 3 - Performing the database search''':

''Significance'': We will put the criteria for significance at 1e-10 (remember: the higher the E-value, the worse the significance).

'''Subquestion''':

Cover the following in your answer:
* What tool(s) and database(s) will be relevant to use?
* Document the results from the different BLAST searches - what works and what does not work?
* You need to copy in small snippets of the BLAST results to document what you observe.
* '''In conclusion''': What kind of enzyme is CLONE12? Gather as much evidence as possible.

'''Blast Job ID (blastn): B4C91VW4016'''

'''Blast Job ID (blastp): B4CMGVXT016'''

 

==Part 4: BLAST'ing Genomes==
<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' do '''NOT''' limit your search to "bacteria" here - now we are actively looking at organism specific queries.
</div>

So far we have been using BLAST to search in the big broad databases that covers at huge set of sequence from a large range of organisms. In this final part of the exercise we will be doing some more focused searches in smaller databases by targeting specific genomes.

Typically this will be useful if you have a gene of known function from one organism (say a cell-cycle controlling gene from Yeast, ''Saccharomyces cerevisiae'') and want to find the human homolog/ortholog to this gene (genes that control cell division are often involved in cancer).

When you have been performing the BLAST searches, you have probably already noticed, that's it possible to search specifically in the Human and Mouse genomes (these database only contains sequences from Human/Mouse). It's also possible to restrict the output from searches in the large databases (e.g. NR) to specific organisms.

A growing number of organisms have been fully sequenced, and the research teams responsible for a large scale genome project typically put up their own Web resources for accessing the data. For example the Yeast genome is principally hosted in the Saccharomyces Genome Database (SGD - www.yeastgenome.org) - it should be noted that SGD also offers BLAST as a means to search the database.



===Genome specific analysis of histones===
====SGD====
Let's do a small study of the relationship between the histones found in Yeast and in Human (evolutionary distance: ~1-1.5 billion years).

Look up the '''HTA2''' gene in SGD (http://www.yeastgenome.org - use the search box at the top of the page). Notice that a brief description about the function of the gene and its protein product is displayed (a huge amount of additional information can be found further down the page - much of it Yeast specific).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.1''':
:What information is given about the relationship between this gene and the gene "HTA1"?

Browse the page and locate the link to the protein sequence. Save the sequence as a file, '''we'll need it in a moment'''.

====NCBI====


Now return to the NCBI blastp page. Set Database to "Reference proteins (refseq_protein)", and enter <tt>Saccharomyces cerevisiae</tt> in the Organism field (and accept the suggestion with taxid:4932).

'''Blast Job ID: B4D0NJAJ016'''

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.2''':
:''(Remember to document your answers)''
:* How many high-confidence hits do we get?
:* Do the hits make sense, from what you have read about HTA2 at the SGD webpage?
:'''Tip:''' click on the Gene links under Related Information (to the right of the alignments) to see the gene names for the protein hits.

The next step is to search the translated version of the human genome.

Do as before, still with Database set to "Reference proteins (refseq_protein)", just enter <tt>Human</tt> in the Organism field.

'''Blast Job ID: B4D73S06016'''

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.3''':
:* How many high-confidence hits (with E-value better than 10-10) are found? (Approximately)
:* What are all the high-confidence hits called?




==Concluding remarks==
Today we have been using BLAST to find a number of homologous genes (and protein-products). If we want to go even deeper into the analysis of the homologs, the next logical step would be to build a dataset of the full-length versions of the sequences we have found (not just the part found by the local alignment in BLAST).

A further analysis could consist of a series of pairwise alignments (for finding out what is similar/different between pairs of sequences) or a multiple alignment which could form the basis of establishing the evolutionary relationship between the entire set of sequences.

BLAST can also be used as way to build a dataset of sequences base on a known "seed" sequence. As we saw in the GenBank exercise, free-text searching in the GenBank can be difficult, and if we for instance wanted to build a dataset of variants of the insulin gene, an easiy way to go around this would be to BLAST the normal version of the insulin against the sequence database of choice, and pick the best matching hits from here.

ExBlast-Answers new

2026-09-22T08:48:16Z

Carol: /* QUESTION 2.3 */

Answers to the BLAST exercise. Values for database sizes etc. retrieved
September 21, 2026

==Part 1: Your first BLAST search==

===QUESTION 1.1===
The main difference is that core_nt excludes most eukaryotic chromosome sequences, making it smaller and faster than nr/nt, so NCBI uses it as the default for more efficient searches.

===QUESTION 1.2===
* ''Which database did you use?''
:core_nt
* ''what is the identifier (Accession)?''
:J00336.1 (Note that the latter was also part of the sequence name for your query sequence!)
* ''what is the alignment score ("Max score in bits")?''
:The max score is 708 bits
* ''what is the percent identity and query coverage?''
:100% and 100%
* ''what is the E-value?''
:0.0 (actually, a number so small that it is rounded off to 0.0)
* ''are there any gaps in the alignment?''
:No, of course not, since the sequences are identical

===QUESTION 1.3===
* ''Which database did you use?''
:core_nt
* ''what is the identifier (Accession)?''
:NM_000207.3
* ''what is the alignment score ("max score")?''
:603
* ''what is the percent identity and query coverage?''
:identity: 94.34% and query coverage: 99%
* ''what is the E-value?''
:7e-168 (meaning 7×10-168)
* ''are there any gaps in the alignment?''
:No

===QUESTION 1.4===
* ''Which database did you use?''
:Human G+T (2 databases)
* ''what is the identifier (Accession)?''
:NM_000207.3, they have the same score and are therefore equally good. Note that these are among the equally good hits found in the previous question.
* ''what is the alignment score ("max score")?''
:603
* ''what is the percent identity and query coverage?''
:identity: 94.34% and query coverage: 99%
* ''what is the E-value?''
:3e-170 (meaning 3×10-170)
* ''are there any gaps in the alignment?''
:No

===QUESTION 1.5===
''What are the sizes (in basepairs) of the databases we used for the two BLAST searches?''

nt(Human): 1,081,933,534,067 letters (= basepairs), Human G+T: 4,017,001,775 letters (= basepairs).

===QUESTION 1.6===

*''What is the ratio between the database sizes in the two BLAST searches?''
: 1,081,933,534,067 / 4,017,001,775 = ~269
*''What is the ratio between the E-values (for the best human hits) in the two BLAST searches? ''
:7e-168 / 3e-170 = ~233
:Note: since the E-values have only three significant digits, you cannot expect to get the exact same result.
:Also note, you can google "7e-168 / 3e-170" directly and the answer will show up in the results.
*''What is the relationship between database size and E-value for hits with identical alignment score?''
:The E-value is directly proportional to the database size.
:Note: Conceptually this means that an alignment within the given score (603 bits) is more SIGNIFICANT in the smaller database. In larger database there is a larger chance of randomly picking up matches.
*''In conclusion: if the database size is doubled, what will happen to the E-value?''
:Each time the database size doubles, the E-value doubles as well.

==Part 2: Assessing the statistical significance of BLAST hits==

===QUESTION 2.1===

****Alignment**** 1
Title: gi|2440392781|emb|OX421481.1| Eilema caniola genome assembly, chromosome: 20
Accession: OX421481
Length: 22119023
Max Score: 44.0
Bits: 40.9604
Identities: 22
Align_length: 22
Gaps: 0
%Ident: 100.00 %
Query Cover: 88 %
E value: 1.88e+00
TTCTGAAAGGTCCTCTCGATAC
||||||||||||||||||||||
TTCTGAAAGGTCCTCTCGATAC

*''What is the typical length of the hits (the alignment length)?'':
:Typically around 17-24 base pairs.

*''What is the typical % identity?'':
:90% - 100%

*''In what range are the bit-scores ("max score)?'':
:typically 30-45 bits.

*''What is the range of the E-values?'':
:1-15
:usually varying from 1 to 50 (occasionally, you might find hits as "good" as 0.1).
:'''Note''' that here we choose to use an E-value threshold of 50.0 as a cut-off. That means that some sequences do not give an output if the E-value of the best alignment is higher than that. Also note that when running BLAST to identify homologous sequences the default cut-off is 0.05.

===QUESTION 2.2===
''What is the biological significance of the hits you found / is there any biological meaning?'':

This makes absolutely NO biological sense(!) The hits are real enough as such, they represent sequences that actually are in the database. But we know that our query sequences are completely random and therefore have no evolutionary relationship with the hits. The only reason we found our hits is that the database is so vast that we for for purely stochastic reasons happen upon sequences that are similar.

The E-values tell us precisely this: As described in the BLAST lecture, the alignment score will follow an extreme value distribution for those sequences that are not related to our query sequences, and the E-value is ''the expected number'' of spurious (unrelated) hits with the given alignment score or better, given the database size.

'''Note:''' Don't be confused by the difference between alignment score and bit score; bit score is simply the alignment score normalized by a constant factor which gives a result expressible in bits.

===QUESTION 2.3===


*''What is the typical length of the alignment and do they contain gaps?'':
:Typically 15-23. Rarely gaps, but several mismatches.

*''What is the range of E-values?'':
:Typically 100-1000

*''Try to inspect a few of the alignments in details ("+" means similar) - do you find any that look plausible, if we for a moment ignore the length/E-value?''
:Yes, maybe. See ''e.g.'' the alignment below, it has 56.5% identities (but it is way too short to be significant, as the E-value tells us).

Generating Random Sequence 1 for BLAST...
DTGVHDAQEIGITCQFQPVCFLFSP
****Alignment**** 1
Title: ref|WP_243336483.1| class III lanthionine synthetase LanKC [Streptomyces sp. AP-93] >gb|MCJ0871799.1| class III lanthionine synthetase LanKC [Streptomyces sp. AP-93]
Accession: WP_243336483
Length: 827
Max Score: 69.0
Bits: 31.187
Identities: 13
Align_length: 23
Gaps: 0
%Ident: 56.52 %
Query Cover: 92 %
E value: 1.12e+02
Query: TGVHDAQEIGITCQFQPVCFLFS
Match: T VHDA I +T QF P C LF+
Subject: TAVHDALLIDLTSQFHPACGLFT

*''If we had used the default E-value cutoff of 10 would any hits have been found?'':
:No (note: the default is actually 0.05 now). Note the difference from the nucleotide database searches (whose E-values were typically in the range 1-50): if we had run BLASTN with an E-value threshold of 1000, we would have had many pages of hits for each query sequence.

===QUESTION 2.4===

*''If we compare the result from BLAST'ing random DNA sequences to random Peptide sequences - which kind of search has the higher risk of returning false positives (results that appear plausible, maybe even significant, but are truly unrelated)?'':
:The risk of getting a false hit (an unrelated sequence with a "decent" E-value) is much larger when working with DNA sequences. Remember than we used 50 as E-value cut-off for BLASTN, while we used 1000 with BLASTP in order to see any hits at all.

==Part 3: using BLAST to transfer functional information by finding homologs==



*''Do we find any significant hits? (E-value?)'':
:Yes, a lot. The first many hits have an E-value of 0.0, and hit #100 is still very significant (3e-98) — note that by default, only the top 100 hits are shown!

*''Are all the best hits the same category of enzymes?'':
:Yes, they are alkaline proteases (except a few that are hypothetical proteins).
:Note that you can click the Accession code for a hit and go directly to the corresponding entry in the database.

*''From what you have seen, what is best for identifying intermediate quality hits - DNA or Protein BLAST?'':
:Protein BLAST (BLASTP). If you have very high quality hits, they can be identified by both methods, but if the evolutionary distance is larger, BLASTP is clearly better.
:Note: Recall from the PyMOL exercises that information between distant genes/proteins are conserved from: Structure > Peptide Sequence > Nucleotide sequence. So when the evolutionary distance is larger, blastp would generally give better hits than blastn.

-->

===QUESTION 3.1===
'''STEP 1 - cleaning up the sequence: '''

*''Subquestion: convert the sequence to FASTA format (manually, in Geany) and quote it in your report.''

>CLONE12
AACGGGCACGGGACGCATGTAGCTGGAACAGTGGCAGCCGTAAATAATAATGGTATCGGA
GTTGCCGGGGTTGCAGGAGGAAACGGCTCTACCAATAGTGGAGCAAGGTTAATGTCCACA
CAAATTTTTAATAGTGATGGGGATTATACAAATAGCGAAACTCTTGTGTACAGAGCCATT
GTTTATGGTGCAGATAACGGAGCTGTGATCTCGCAAAATAGCTGGGGTAGTCAGTCTCTG
ACTATTAAGGAGTTGCAGAAAGCTGCGATCGACTATTTCATTGATTATGCAGGAATGGAC
GAAACAGGAGAAATACAGACAGGCCCTATGAGGGGAGGTATATTTATAGCTGCCGCCGGA
AACGATAACGTTTCCACTCCAAATATGCCTTCAGCTTATGAACGGGTTTTAGCTGTGGCC
TCAATGGGACCAGATTTTACTAAGGCAAGCTATAGCACTTTTGGAACATGGACTGATATT
ACTGCTCCTGGCGGAGATATTGACAAATTTGATTTGTCAGAATACGGAGTTCTCAGCACT
TATGCCGATAATTATTATGCTTATGGAGAGGGAACATCCATGGCTTGTCCACATGTCGCC
GGCGCCGCC

'''STEP 2 - thinking about the task: '''

*''Subquestion: Give a summary of your considerations.''
**''Based on the information given: is the sequence protein-coding? ''
::Yes — we know this because the PCR primers used to clone the sequence target '''known enzymes'''. Therefore, it will make sense to try to translate the sequence using Virtual Ribosome.
:*''If it is, can you trust it will contain both a START and STOP codon? ''
::No — the PCR primers used to clone the sequence target '''the middle of the sequence''', in other words we must assume that our sequence is a fragment. Therefore, the ORF finder in Virtual Ribosome should be set to Start codon: None.
:*''Do we know if the sequence is sense or anti-sense? ''
::No — the PCR process amplifies a stretch of double-stranded DNA. Therefore, we should let Virtual Ribosome search in '''all 6 reading frames'''.

'''STEP 3 - Performing the database search''':
We want to use BLAST to search the large databases. Let's therefore try the following:
# BLASTN
# Translate to protein (using Virtual Ribosome).
# BLASTP
Both when doing BLASTN and BLASTP we will use the NR database in order to search as broadly as possible. It would not make sense to use an organism-specific database when we don't know which organism our sequence stems from.

1) BLASTN. When trying BLASTN against NR we get some borderline significant results, but observe how small the query coverage percentages are (check also the Graphic Summary tab!).

[[image:NCBI BlastN_CLONE12 new_version.png]]
[[image:NCBI BlastN_CLONE12 Graphic Summary.png]]

There is simply nothing in the entire NR database that has enough similarity to our whole query sequence. A search on the DNA level is only suited for finding very close hits.

2) Translate using Virtual Ribosome with the settings we chose under Step 2 above.

The result from the ORF finder:

VIRTUAL RIBOSOME
----------------
Translation table: Standard SGC0

>CLONE12_rframe1_ORF
Reading frame: 1

N G H G T H V A G T V A A V N N N G I G V A G V A G G N G S
5' AACGGGCACGGGACGCATGTAGCTGGAACAGTGGCAGCCGTAAATAATAATGGTATCGGAGTTGCCGGGGTTGCAGGAGGAAACGGCTCT 90
..........................................................................................

T N S G A R L M S T Q I F N S D G D Y T N S E T L V Y R A I
5' ACCAATAGTGGAGCAAGGTTAATGTCCACACAAATTTTTAATAGTGATGGGGATTATACAAATAGCGAAACTCTTGTGTACAGAGCCATT 180
.....................>>>..................................................................

V Y G A D N G A V I S Q N S W G S Q S L T I K E L Q K A A I
5' GTTTATGGTGCAGATAACGGAGCTGTGATCTCGCAAAATAGCTGGGGTAGTCAGTCTCTGACTATTAAGGAGTTGCAGAAAGCTGCGATC 270
.........................................................)))............)))...............

D Y F I D Y A G M D E T G E I Q T G P M R G G I F I A A A G
5' GACTATTTCATTGATTATGCAGGAATGGACGAAACAGGAGAAATACAGACAGGCCCTATGAGGGGAGGTATATTTATAGCTGCCGCCGGA 360
........................>>>..............................>>>..............................

N D N V S T P N M P S A Y E R V L A V A S M G P D F T K A S
5' AACGATAACGTTTCCACTCCAAATATGCCTTCAGCTTATGAACGGGTTTTAGCTGTGGCCTCAATGGGACCAGATTTTACTAAGGCAAGC 450
........................>>>....................................>>>........................

Y S T F G T W T D I T A P G G D I D K F D L S E Y G V L S T
5' TATAGCACTTTTGGAACATGGACTGATATTACTGCTCCTGGCGGAGATATTGACAAATTTGATTTGTCAGAATACGGAGTTCTCAGCACT 540
...............................................................)))........................

Y A D N Y Y A Y G E G T S M A C P H V A G A A
5' TATGCCGATAATTATTATGCTTATGGAGAGGGAACATCCATGGCTTGTCCACATGTCGCCGGCGCCGCC 609
.......................................>>>...........................

('''Tip:''' Remember that you can get the sequence in FASTA format via the FASTA link on the result page):

>CLONE12_rframe1_ORF
NGHGTHVAGTVAAVNNNGIGVAGVAGGNGSTNSGARLMSTQIFNSDGDYTNSETLVYRAI
VYGADNGAVISQNSWGSQSLTIKELQKAAIDYFIDYAGMDETGEIQTGPMRGGIFIAAAG
NDNVSTPNMPSAYERVLAVASMGPDFTKASYSTFGTWTDITAPGGDIDKFDLSEYGVLST
YADNYYAYGEGTSMACPHVAGAA

3) BLASTP



We get several very significant hits. When looking at the top hits and disregarding "hypothetical" and "uncharacterized" proteins, we can see that the rest are almost all serine proteases. Some of them are described as belonging to the of the S8 family.

[[image:NCBI BlastP_CLONE12_rframe1_ORF new version.png]]

Let's take a closer look at the first hit that is not "uncharacterized":
[[image:NCBI_BlastP_CLONE12_best_hit.png]]

Note that although it is not a perfect hit (our query sequence not existing in the database) it looks reasonable: the alignment covers a large part of the query with Identity of 54% and Similarity (Positives) of 69%.

Taken together with the fact that almost all the best non-hypothetical hits are serine proteases, we have a very strong indication that our mystery sequence, CLONE12, is a peptidase or protease of the S8 family.

==Part 4: BLAST'ing Genomes==

===QUESTION 4.1===
''What information is given about the relationship between this gene and the gene "HTA1"?''

They are nearly identical ("one of two nearly identical (see also HTA1) subtypes").

Protein sequence:
>YBL003C
MSGGKGGKAGSAAKASQSRSAKAGLTFPVGRVHRLLRRGNYAQRIGSGAPVYLTAVLEYL
AAEILELAGNAARDNKKTRIIPRHLQLAIRNDDELNKLLGNVTIAQGGVLPNIHQNLLPK
KSAKTAKASQEL*

===QUESTION 4.2===
*''How many high-confidence hits do we get?'':
:3 — HTA1, HTA2 and HTZ1.
:Note: If you click on the Gene links for the two top hits, you will see that one is HTA1 and the other is HTA2.

*''Do the hits make sense, from what you have read about HTA2 at the SGD webpage?'':
:Yes; HTA1 and HTA2 are indeed nearly identical (only 2 amino acids differ).

===QUESTION 4.3===

*''How many high-confidence hits (with E-value better than 10-10) are found?''
:Answer: 29.

2026-09-21T13:40:09Z

Carol: /* Useful Links */

Exercise written by Rasmus Wernersson and modified by Henrik Nielsen and Carolina Barra.

==Introduction==

In this exercise you will learn by doing to BLAST DNA and protein sequences.
BLAST is short for '''B'''asic '''L'''ocal '''A'''lignment '''S'''earch '''T'''ool.

After completion of the exercise you should be able to choose an appropriate BLAST search tool selecting the correct set of parameters, interpret a BLAST result, explain what an E-value means and how database size affects it, and decide what biological conclusion is justified.

For the exercise, we will be using two approaches; the Web interface to '''BLAST''' hosted by the '''NCBI''', and using the '''BioPython library'' in a Notebook using '''Google Colab'''. There are several advantages to each approach:

The Web interface is very intuitive as it provides a graphical summary and links to other databases, while the programmatic way is more flexible allowing us to further customize the search, the output, and do many searches in a single run for a batch of query sequences.

===Useful Links===
* NCBI BLAST main page: http://blast.ncbi.nlm.nih.gov/
* NCBI [http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs BLAST help pages]

 

<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' BLAST is a quite computationally intensive algorithm, and we have in recent years run into issues with overburdening the NCBI server, with 150+ students submitting jobs at the same time from the same DTU IP address. We have therefore implemented a few optimization/workarounds.

:If you still experience long waiting times you can try an alternative internet connection, such as a mobile hotspot. Alternatively, you can copy the pre-saved searches and use the link from BLAST called '''Saved Strategies'''. Please notice that these links will only work for the duration of the exercise session.

[[Image:BLAST_Saved_Strategies.png|center|600px|border]]

 
</div>

==Part 1: Your first BLAST search==

Below is the mRNA sequence for insulin from a close relative, the common macaque (''Macaca fascicularis'').

>J00336.1 Monkey (M. fascicularis) preproinsulin mRNA, complete cds
AACCCTCCGGGACAGGCTGCATCAGAAGAGGTCAGCAAGCAGGTCACTGTCCTTCGGCTATGGCCCTGTG
GATGCGCCTCCTGCCCCTGCTGGCGCTGCTGGCCCTCTGGGGACCTGACCCGGCCCCGGCCTTTGTGAAC
CAGCACCTGTGCGGCTCCCACCTAGTGGAAGCTCTCTACCTGGTGTGCGGGGAACGAGGCTTCTTCTACA
CACCCAAGACCCGCCGGGAGGCAGAGGACCCTCAGGTGGGGCAGGTGGAGCTGGGCGGGGGCCCTGGCGC
AGGCAGCCTGCAGCCCTTGGCGCTGGAGGGGTCCCTGCAGAAGCGCGGCATCGTGGAGCAGTGCTGCACC
AGCATCTGCTCCCTCTACCAGCTGGAGAACTACTGCAACTAA

We will use BLAST to determine whether this sequence looks like the human mRNA for insulin. There are two ways we can do this:

:A. Search the entire database and look for human hits in the results
:B. Specifically search the human part of the database

We will try both possibilities.

=== Search against a broad database ===

Follow the "nucleotide blast" link from the main [https://blast.ncbi.nlm.nih.gov/Blast.cgi BLAST] page.

After pasting the sequence, we are asked to "Choose a search set" or database to perform our search.

In the previous years the students were using a database called "'''nr/nt'''", but today the default is a database called Core Nucleotide database "'''core_nt'''". We want to learn more about these broad databases of nucleotides.
You can click on the question mark of the databases to learn more about them or you can ask your favourite AI as a sparring partner about what is the difference among them.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.1''': What is the main difference of the two databases in regards to content? Why NCBI offers now core_nt as the default. AI loves long and nice explanations that will help you understand, but the TAs have many submissions to review. Please rephrase the answer with your own words in one line.



Based on the answer of Question 1.1, choose the most efficient database from the two.
In the section "Program Selection" select the option "Somewhat similar sequences (blastn)"
Click the BLAST button to launch the search.

After the search has been completed, make yourself familiar with the BLAST output page. Note that you can toggle between hiding and showing each tab by clicking on the title: Descriptions, Graphic Summary, Alignments and Taxonomy (try it!).


First, take a look at the best hit. Since our search sequence (the query) was taken from GenBank which is part of this broad database, we should find an identical sequence in the search. Make sure this is the case!

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.2''':
:Answer the following questions about the best hit:
:* Which database did you use?
:* What is the identifier (Accession)?
:* What is the alignment score ("max score")?
:* What is the percent identity and query coverage?
:* What is the E-value?
:* Are there any gaps in the alignment?

'''Blast Job ID: B2C1HZPV014'''


Then, find the best hit from human (''Homo sapiens'') that is ''not'' a synthetic construct. ('''Tip:''' you can press Ctrl-F in most browsers to search in the page).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.3''':

:Answer the same questions as before about the human hit you found now.


=== Search against Human G+T ===

'''Note:''' In this context, G+T does not mean Gin and Tonic.

Open ''a new window/tab'' with the BLAST home page. Make a new BLASTN search with the same query sequence, this time with Database set to Human genomic + transcript (Human G+T). Remember again to select Somewhat similar sequences (blastn) under Program Selection. Consider the best hit.

'''Note:''' even though you may not have found exactly the same database entry in the two searches, the ''alignment'' should be the same. Make sure this is the case by comparing the actual alignments in the two windows where you made the searches.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.4''':
Answer the same questions as before about the best hit you found in this search.

'''Blast Job ID: B2CCWWRV014'''



===Concerning database size and E-values===

When answering the previous two questions, you may have noticed that the E-value changed, while the alignment score did not. We will now investigate this further.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.5''':
What are the sizes (in basepairs) of the databases we used for the two BLAST searches? ('''Tip:''' Expand the "Search summary" section near the top by clicking it).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.6''':
:*What is the ratio between the database sizes in the two BLAST searches?
:*What is the ratio between the E-values (for the best human hits) in the two BLAST searches?
:*What is the relationship between database size and E-value for hits with identical alignment scores?
:*In conclusion: if the database size is doubled, what will happen to the E-value?

==Part 2: Assessing the statistical significance of BLAST hits==
<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' Limit your search to "bacteria" (taxid: 2) in ALL of this section (PART 2) to make the BLAST searches run quicker.
</div>

As discussed in the lecture, there will be a risk of getting false positive results (hits to sequences that are not related to our input sequence) by purely stochastic means. In this first part of the exercise, we will be investigating this further, by examining what happens when we submit randomly generated sequences to BLAST searches.

Rather than giving out a set of pre-generated DNA/Peptide sequences where you only have our word for their randomness, you'll be generating your own random sequences.
It is important to understand that these computer-generated sequences are ''totally random'', just as if you were rolling a die to determine each nucleotide/amino acid in each sequence.

In this part of the exercise, we will run BLAST through a Python Notebook from Google Colab. We will install the package called BioPython that will connect to the NCBI servers where BLAST is hosted.

[[File:Blast exercise in Colab.png|center|600px|border]]

For our purpose, there are several advantages to this approach:

* We don't need to install it on our computers, but we can still perform custom searches which can be fine-grained control of DATA and workflow (everything can be scripted/automated) than running BLAST through a web interface
* NCBI offers direct access to preformatted BLAST databases of all the data that they host:
** GenBank (+ derivates)
** Full Genome database
** Protein databases (Both from translated GenBank and UniProt)

===Links===
* A few notes on how to use a [https://colab.research.google.com/notebooks/basic_features_overview.ipynb#scrollTo=WUtu4316QSHL Google Colab]

<div style="background-color: lavender; border: solid thin grey;">
'''IMPORTANT:''' Check that your notebook is always connected to the server.

If your laptop is idle for some time, you may get the error 'runtime disconnected'

[[File:Runtime_disconnected.png|center|300px|border]]
 

To re-connect to the server

[[File:Reconnect_to_a_hosted_runtime.png|center|600px|border]]

Remember that after is disconnected we will need to run all the code cells from that part of the exercise. This could take some minutes for the blast searches, so try to avoid having the screen idle for a long time :)

 
</div>

* Biopython [https://biopython.org/DIST/docs/tutorial/Tutorial.html tutorial]

* Biopython [https://github.com/biopython/biopython/blob/a09f38b468cbedaf25aec6cfed1475b0849011e6/Bio/Blast/NCBIWWW.py examples]
 

==Part 2.0: Create your first Python Notebook==

[https://colab.research.google.com/drive/1DpqOH8O_KK0rEohjv-20CCx2THwCq89h#scrollTo=_7tJ_S6wc-Os Click Here]

We are going to make a local copy of this notebook on our own drive.
* First you will need to log in or create a Google account.
* Make a copy of the Colab file in your own drive. Go to File > Save in Drive

[[File:save_a_copy_in_drive.png|center|300px|border]]

Now you can start modifying the file!

* First, we will install BioPython in our notebook. This will only take a few seconds, but note, that this installation is temporary. If your notebook gets idle and disconnected you will need to run this code cell again.

! pip install biopython

* Now import Blast tools from BioPython

from Bio.Blast import NCBIWWW
help(NCBIWWW.qblast)

===Random DNA sequences and BLASTN===

*Generate 10 random DNA sequence of length 25bp using Google Colab

[[Image:Office-notes-line_drawing.png|30px|left]]

We will now do a BLASTN search using the random sequence as query. Select the option "blastn" in the section "Program Selection". Choose "Nucleotide Collection (nr/nt)" as the search database.

'''VERY IMPORTANT''':
For this special situation where we BLAST small artificial sequences, we need to turn off some the automatics NCBI incorporate when short sequences are detected. Otherwise we'll not be able to see the intended results:

* Extend the "Algorithm parameters" section (see the screen shot below) in order to gain access to fine-tuning the options.
*# '''Deselect''' the "Automatically adjust parameters for short input sequences" option.
*# Set the E-value cut-off ("Expect threshold") to '''50'''

[[file:Blastn_cropped+circle.png‎|center|frame|'''Remember to adjust the BLAST settings''']]

* Inspect the results.

:'''QUESTION 2.1''':
:Answer the following small questions, and '''document your findings''' by pasting in examples of alignments / text snippets from the overview plots generated on google colab:

:* What is the typical length of the hits (the alignment length)?
:* What is the typical % identity?
:* In what range is the bit-scores ("max score")?
:** ''Notice: This is conceptually the same as the "alignment score" we have already met in the pairwise alignment exercise''.
:* What is the range of the E-values?

:'''QUESTION 2.2''':
:*What is the '''biological''' significance of these hits / is there any biological meaning?

===Random protein sequences and BLASTP===
Now it's time to work with a set of '''protein sequences''': Generate 10 peptide sequences of length 25aa using Google colab.

* '''Notice 1:''' The distribution of amino acids will be equal (5% prob) and this is different from true biological sequences - however this is not important for this first part of the exercise.
* '''Notice 2:''' Please recall from the lecture that the way <tt>BLASTP</tt> selects candidate sequences for full Smith-Waterman alignment is different from <tt>BLASTN</tt>. (<tt>BLASTN</tt> - a single short (11 bp +) perfect match hit is needed. <tt>BLASTP</tt> - a pair of "near match" hits of 3 aa within a 40 aa window is needed).

[[Image:Office-notes-line_drawing.png|30px|left]]

Choose blastp as the algorithm to use, the "NR" database (this is the protein version, consisting of translated CDS'es, UniProt etc).

'''VERY IMPORTANT''': We also need to tweak the parameters this time - in the "Algorithm Parameters" section select BLOSUM62 as the alignment matrix to use and set the "Expect threshold" to 1000 (default: 10).

* Perform the BLAST search.
* Inspect the results.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 2.3''':
:''(Remember to '''document your answers''' in the same manner as Q2.2)''

:* What is the typical length of the alignment and do they contain gaps?
:* What is the range of E-values?
:* Try to inspect a few of the alignments in details ("+" means similar sequences) - do you find any that look plausible, if we for a moment ignore the length/E-value?
:* If we had used the default E-value cut-off of 10 would any hits have been found?

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 2.4''':
:* If we compare the result from BLAST'ing random DNA sequences to random Peptide sequences - which kind of search has the higher risk of returning false positives (results that appear plausible, maybe even significant, but are truly unrelated)?
:** Remember to take E-values into your consideration.

 

==Part 3: using BLAST to transfer functional information by finding homologs==
<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' The following parts can either be done on the google colab or in the webserver. If you choose the web server version remember to limit your search to "bacteria" (taxid: 2) in ALL of this section (PART 3) to make the BLAST searches run quicker. (The organisms we're looking for all belong to the "Bacteria" domain of life, so this restriction is OK).
</div>

===Homo-, Ortho- and Paralogs===

One of the most common ways to use BLAST as a tool is in the situation where you have a sequence of '''unknown function''', and want to '''find out which function it has'''. Since a large amount of sequence data has been gathered over the years, chances are that an '''evolutionarily related''' sequence with known function has already been identified. In general, such a related sequence is known as a "'''homolog'''".

Homo-, Ortho- and Paralogs:
* A '''Homolog''' is a general term that describes a sequence that is related by any evolutionary means.
* An '''Ortholog''' ("Ortho" = True) is a sequence that is "the same gene" in a different organism: The sequences shared a single common ancestor sequence, and have now diverged through speciation (e.g. the Alpha-globin gene in Human and Mouse).
* A '''Paralog''' arises due to a gene duplication within a species. For example, Alpha- and Beta-globin are paralogs.

[[File:Homo_Ortho_Para-log.gif|center|frame|''Image source: [http://www.thegreatgoodplace.com/tt/gwlee/126 gwLee's blog]'' ]]

Notice that in both cases it's possible to transfer information, for example, about gene family / protein domains.
We have already touched upon the comparison of (potentially) evolutionarily related sequences in the pairwise alignment exercise. However, this time we do not start with two sequences we assume are related, but instead, we start with a single sequence ("query sequence") which we will use to search the databases for homologs (we often informally speak of "BLAST hits", when discussing the sequences found).

 



===BLAST example 2===

The sequence is a '''DNA fragment''' from an unknown non-cultivatable microorganism. It was cloned and sequenced directly from DNA extracted from a soil-sample, and it goes by the poetic name "CLONE12". It was amplified using degenerated PCR primers that target the middle ("core cloning") of the sequence of '''a group of known enzymes'''. (I can guarantee this particular sequence is not in the BLAST databases, since I have cloned and sequenced it myself, and it has never been submitted to GenBank).

LOCUS CLONE12.DNA 609 BP DS-DNA UPDATED 06/14/98
DEFINITION UWGCG file capture
ACCESSION -
KEYWORDS -
SOURCE -
COMMENT Non-sequence data from original file:
BASE COUNT 174 A 116 C 162 G 157 T 0 OTHER
ORIGIN ?
clone12.dna Length: 609 Jun 13, 1998 - 03:39 PM Check: 6014 ..
1 AACGGGCACG GGACGCATGT AGCTGGAACA GTGGCAGCCG TAAATAATAA TGGTATCGGA
61 GTTGCCGGGG TTGCAGGAGG AAACGGCTCT ACCAATAGTG GAGCAAGGTT AATGTCCACA
121 CAAATTTTTA ATAGTGATGG GGATTATACA AATAGCGAAA CTCTTGTGTA CAGAGCCATT
181 GTTTATGGTG CAGATAACGG AGCTGTGATC TCGCAAAATA GCTGGGGTAG TCAGTCTCTG
241 ACTATTAAGG AGTTGCAGAA AGCTGCGATC GACTATTTCA TTGATTATGC AGGAATGGAC
301 GAAACAGGAG AAATACAGAC AGGCCCTATG AGGGGAGGTA TATTTATAGC TGCCGCCGGA
361 AACGATAACG TTTCCACTCC AAATATGCCT TCAGCTTATG AACGGGTTTT AGCTGTGGCC
421 TCAATGGGAC CAGATTTTAC TAAGGCAAGC TATAGCACTT TTGGAACATG GACTGATATT
481 ACTGCTCCTG GCGGAGATAT TGACAAATTT GATTTGTCAG AATACGGAGT TCTCAGCACT
541 TATGCCGATA ATTATTATGC TTATGGAGAG GGAACATCCA TGGCTTGTCC ACATGTCGCC
601 GGCGCCGCC
//

[[Image:Office-notes-line_drawing.png|30px|left]][[Image:Cogs_brain.png|right|150px]]
:'''QUESTION 3.1 (Long question - read all)''':
: ''Your task is now to find out '''what kind of enzyme''' this sequence is likely to encode, '''using the methods''' you have learned''.

'''INSTRUCTIONS''': You are free to write the combined answer to this question in a free-style essay-like fashion - just be sure to include the subquestions in your answers. In an exam situation you will need to put all the clues together yourself, reason about the tools/databases to use, and document your findings.

'''STEP 1 - cleaning up the sequence''':

The sequence is (more or less) in GenBank format and the NCBI BLAST server expects the input to be in FASTA format, or to be "raw" unformatted sequence.

* There are two solutions to this:
** Copy the sequence into a text-editor and manually create a FASTA file ("search and replace" and/or "rectangular selection" is useful for the reformatting). This is the most robust solution: it will always work. (Look at [[Plain text files and Geany|the Geany exercise]] for a reminder of how to do this).
** Hope the creators of the web-server you're using were kind enough to automatically remove non-DNA letters (paste in ONLY the DNA lines) - this turns out to be the case for both NCBI BLAST and VirtualRibosome, but it ''cannot be universally relied upon''.

'''Subquestion''': convert the sequence to FASTA format (manually, in Geany) and quote it in your report.

'''STEP 2 - thinking about the task''':

Consider the following before you start on solving this task:
* Based on the information given: is the sequence protein-coding?
* If it is, can you trust it will contain both a START and STOP codon?
* Do we know if the sequence is sense or anti-sense?
and think which consequences the answers to these questions should have for your choice of methods and parameters.

'''Subquestion''': Give a summary of your considerations.

'''STEP 3 - Performing the database search''':

''Significance'': We will put the criteria for significance at 1e-10 (remember: the higher the E-value, the worse the significance).

'''Subquestion''':

Cover the following in your answer:
* What tool(s) and database(s) will be relevant to use?
* Document the results from the different BLAST searches - what works and what does not work?
* You need to copy in small snippets of the BLAST results to document what you observe.
* '''In conclusion''': What kind of enzyme is CLONE12? Gather as much evidence as possible.

Blast Job ID (blastn): DU9XHV8H014

Blast Job ID (blastp): DU9T4K44014

 

==Part 4: BLAST'ing Genomes==
<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' do '''NOT''' limit your search to "bacteria" here - now we are actively looking at organism specific queries.
</div>

So far we have been using BLAST to search in the big broad databases that covers at huge set of sequence from a large range of organisms. In this final part of the exercise we will be doing some more focused searches in smaller databases by targeting specific genomes.

Typically this will be useful if you have a gene of known function from one organism (say a cell-cycle controlling gene from Yeast, ''Saccharomyces cerevisiae'') and want to find the human homolog/ortholog to this gene (genes that control cell division are often involved in cancer).

When you have been performing the BLAST searches, you have probably already noticed, that's it possible to search specifically in the Human and Mouse genomes (these database only contains sequences from Human/Mouse). It's also possible to restrict the output from searches in the large databases (e.g. NR) to specific organisms.

A growing number of organisms have been fully sequenced, and the research teams responsible for a large scale genome project typically put up their own Web resources for accessing the data. For example the Yeast genome is principally hosted in the Saccharomyces Genome Database (SGD - www.yeastgenome.org) - it should be noted that SGD also offers BLAST as a means to search the database.



===Genome specific analysis of histones===
====SGD====
Let's do a small study of the relationship between the histones found in Yeast and in Human (evolutionary distance: ~1-1.5 billion years).

Look up the '''HTA2''' gene in SGD (http://www.yeastgenome.org - use the search box at the top of the page). Notice that a brief description about the function of the gene and its protein product is displayed (a huge amount of additional information can be found further down the page - much of it Yeast specific).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.1''':
:What information is given about the relationship between this gene and the gene "HTA1"?

Browse the page and locate the link to the protein sequence. Save the sequence as a file, '''we'll need it in a moment'''.

====NCBI====


Now return to the NCBI blastp page. Set Database to "Reference proteins (refseq_protein)", and enter <tt>Saccharomyces cerevisiae</tt> in the Organism field (and accept the suggestion with taxid:4932).

Blast Job ID: FR6NDFMJ013

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.2''':
:''(Remember to document your answers)''
:* How many high-confidence hits do we get?
:* Do the hits make sense, from what you have read about HTA2 at the SGD webpage?
:'''Tip:''' click on the Gene links under Related Information (to the right of the alignments) to see the gene names for the protein hits.

The next step is to search the translated version of the human genome.

Do as before, still with Database set to "Reference proteins (refseq_protein)", just enter <tt>Human</tt> in the Organism field.

Blast Job ID: FR65XVA0016

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.3''':
:* How many high-confidence hits (with E-value better than 10-10) are found? (Approximately)
:* What are all the high-confidence hits called?




==Concluding remarks==
Today we have been using BLAST to find a number of homologous genes (and protein-products). If we want to go even deeper into the analysis of the homologs, the next logical step would be to build a dataset of the full-length versions of the sequences we have found (not just the part found by the local alignment in BLAST).

A further analysis could consist of a series of pairwise alignments (for finding out what is similar/different between pairs of sequences) or a multiple alignment which could form the basis of establishing the evolutionary relationship between the entire set of sequences.

BLAST can also be used as way to build a dataset of sequences base on a known "seed" sequence. As we saw in the GenBank exercise, free-text searching in the GenBank can be difficult, and if we for instance wanted to build a dataset of variants of the insulin gene, an easiy way to go around this would be to BLAST the normal version of the insulin against the sequence database of choice, and pick the best matching hits from here.

2026-09-18T16:36:35Z

Carol: /* Search against a broad database */

Exercise written by Rasmus Wernersson and modified by Henrik Nielsen and Carolina Barra.

==Introduction==

In this exercise you will learn by doing to BLAST ('''B'''asic '''L'''ocal '''A'''lignment '''S'''earch '''T'''ool) DNA and protein sequences.
After completion of the exercise you should be able to choose an appropriate BLAST search tool selecting the correct set of parameters, interpret a BLAST result, explain what an E-value means and how database size affects it, and decide what biological conclusion is justified

For the exercise, we will be using the Web interface to '''BLAST hosted by the NCBI''' and from a Python Notebook from Google Colab. There are several advantages to each approach:

The Web interface is very intuitive as it provides a graphical summary and links to other databases, while the programmatic way is more flexible allowing us to further customize the search, the output, and do many searches in a single run for a batch of query sequences.

===Useful Links===
* NCBI BLAST main page: http://blast.ncbi.nlm.nih.gov/
* NCBI [http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs BLAST help pages]

 

<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' BLAST is a quite computationally intensive algorithm, and we have in recent years run into issues with overburdening the NCBI server, with 150+ students submitting jobs at the same time from the same DTU IP address. We have therefore implemented a few optimization/workarounds, that it is '''important you remember to follow'''. In some of the sections below, you will be asked to limit your search to a certain subset of the BLAST database (e.g. only search in the "bacterial" part of the NR database). This will limit the amount of data to search through and will make the search finish faster.

[[Image:BLAST_limit_search.png|center|600px|border]]

:'''NOTE: if you still experience long waiting times you can try an alternative internet connection, such as a mobile hotspot. Alternatively, you can click on the provided links to get the expected output. Please notice that these links will only work for the duration of the exercise session.'''

 
</div>

==Part 1: Your first BLAST search==

Below is the mRNA sequence for insulin from a South American rodent, the Degu (''Octodon degus'').

>gi|202471|gb|M57671.1|OCOINS Octodon degus insulin mRNA, complete cds
GCATTCTGAGGCATTCTCTAACAGGTTCTCGACCCTCCGCCATGGCCCCGTGGATGCATCTCCTCACCGT
GCTGGCCCTGCTGGCCCTCTGGGGACCCAACTCTGTTCAGGCCTATTCCAGCCAGCACCTGTGCGGCTCC
AACCTAGTGGAGGCACTGTACATGACATGTGGACGGAGTGGCTTCTATAGACCCCACGACCGCCGAGAGC
TGGAGGACCTCCAGGTGGAGCAGGCAGAACTGGGTCTGGAGGCAGGCGGCCTGCAGCCTTCGGCCCTGGA
GATGATTCTGCAGAAGCGCGGCATTGTGGATCAGTGCTGTAATAACATTTGCACATTTAACCAGCTGCAG
AACTACTGCAATGTCCCTTAGACACCTGCCTTGGGCCTGGCCTGCTGCTCTGCCCTGGCAACCAATAAAC
CCCTTGAATGAG

We will use BLAST to determine whether this sequence looks like the human mRNA for insulin. There are two ways we can do this:

:A. Search the entire database and look for human hits in the results
:B. Specifically search the human part of the database

We will try both possibilities.

=== Search against a broad database ===

Follow the "nucleotide blast" link from the main [https://blast.ncbi.nlm.nih.gov/Blast.cgi BLAST] page.

After pasting the sequence, we are asked to "Choose a search set" or database to perform our search.

In the previous years the students were using a database called "'''nr/nt'''", but today the default is a database called Core Nucleotide database "'''core_nt'''". We want to learn more about these broad databases of nucleotides.
You can click on the question mark of the databases to learn more about them or you can ask your favourite AI as a sparring partner about what is the difference among them.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.1''': What is the main difference of the two databases in regards to content? Why NCBI offers now core_nt as the default. AI loves long and nice explanations that will help you understand, but the TAs have many submissions to review. Please rephrase the answer with your own words in one line.



Based on the answer of Question 1, choose the most efficient database from the two.
In the section "Program Selection" select the option "Somewhat similar sequences (blastn)"
Click the BLAST button to launch the search.

After the search has been completed, make yourself familiar with the BLAST output page. Note that you can toggle between hiding and showing each tab by clicking on the title (try it!).

'''Descriptions'''
:A table with a one-line description of each hit with some alignment statistics.

'''Graphic Summary'''
:Each hit is represented by a line showing which part of the query sequence the alignment covers. The lines are coloured according to the alignment score.

'''Alignments'''
:The actual alignments between the query and the database hits.

The columns in the '''Descriptions''' tab are:

* Scientific Name — The species name of the matched sequence
* Description — the description line from the database
* Max score — the alignment score of the best match (local alignment) between the query and the database hit
* Total score — the sum of alignment scores for all matches (alignments) between the query and the database hit (if there is only one match per hit, these two scores are identical)
* Query cover — the percentage of the query sequence that is covered by the alignment(s)
* E value — the Expect value calculated from the Max score (''i.e.'' the number of ''unrelated'' hits with that score or better you would expect to find for random reasons)
* Per. Ident — the percent identity in the alignment(s)
* Acc. Len — the accession length
* Accession - the accession ID of the matched sequence

First, take a look at the best hit. Since our search sequence (the query) was taken from GenBank which is part of this broad database, we should find an identical sequence in the search. Make sure this is the case!

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.2''':
:Answer the following questions about the best hit:
:* Which database did you use?
:* what is the identifier (Accession)?
:* what is the alignment score ("max score")?
:* what is the percent identity and query coverage?
:* what is the E-value?
:* are there any gaps in the alignment?

'''Blast Job ID: XXX'''

Then, find the best hit from human (''Homo sapiens'') that is ''not'' a synthetic construct. ('''Tip:''' you can press Ctrl-F in most browsers to search in the page).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.2''':

:Answer the same questions as before about the hit you found now.

=== Search against Human G+T ===

'''Note:''' In this context, G+T does not mean Gin and Tonic.

Open ''a new window/tab'' with the BLAST home page. Make a new BLASTN search with the same query sequence, this time with Database set to Human genomic + transcript (Human G+T). Remember again to select Somewhat similar sequences (blastn) under Program Selection. Consider the best hit.

'''Note:''' even though you may not have found exactly the same database entry in the two searches, the ''alignment'' should be the same. Make sure this is the case by comparing the actual alignments in the two windows where you made the searches.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.3''':
Answer the same questions as before about the best hit you found in this search.

Blast Job ID: DU93YZB1014

===Concerning database size and E-values===

When answering the previous two questions, you may have noticed that the E-value changed, while the alignment score did not. We will now investigate this further.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.4''':
What are the sizes (in basepairs) of the databases we used for the two BLAST searches? ('''Tip:''' Expand the "Search summary" section near the top by clicking it).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 1.5''':
:*What is the ratio between the database sizes in the two BLAST searches?
:*What is the ratio between the E-values (for the best human hits) in the two BLAST searches?
:*What is the relationship between database size and E-value for hits with identical alignment scores?
:*In conclusion: if the database size is doubled, what will happen to the E-value?

==Part 2: Assessing the statistical significance of BLAST hits==
<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' Limit your search to "bacteria" (taxid: 2) in ALL of this section (PART 2) to make the BLAST searches run quicker.
</div>

As discussed in the lecture, there will be a risk of getting false positive results (hits to sequences that are not related to our input sequence) by purely stochastic means. In this first part of the exercise, we will be investigating this further, by examining what happens when we submit randomly generated sequences to BLAST searches.

Rather than giving out a set of pre-generated DNA/Peptide sequences where you only have our word for their randomness, you'll be generating your own random sequences.
It is important to understand that these computer-generated sequences are ''totally random'', just as if you were rolling a die to determine each nucleotide/amino acid in each sequence.

In this part of the exercise, we will run BLAST through a Python Notebook from Google Colab. We will install the package called BioPython that will connect to the NCBI servers where BLAST is hosted.

[[File:Blast exercise in Colab.png|center|600px|border]]

For our purpose, there are several advantages to this approach:

* We don't need to install it on our computers, but we can still perform custom searches which can be fine-grained control of DATA and workflow (everything can be scripted/automated) than running BLAST through a web interface
* NCBI offers direct access to preformatted BLAST databases of all the data that they host:
** GenBank (+ derivates)
** Full Genome database
** Protein databases (Both from translated GenBank and UniProt)

===Links===
* A few notes on how to use a [https://colab.research.google.com/notebooks/basic_features_overview.ipynb#scrollTo=WUtu4316QSHL Google Colab]

<div style="background-color: lavender; border: solid thin grey;">
'''IMPORTANT:''' Check that your notebook is always connected to the server.

If your laptop is idle for some time, you may get the error 'runtime disconnected'

[[File:Runtime_disconnected.png|center|300px|border]]
 

To re-connect to the server

[[File:Reconnect_to_a_hosted_runtime.png|center|600px|border]]

Remember that after is disconnected we will need to run all the code cells from that part of the exercise. This could take some minutes for the blast searches, so try to avoid having the screen idle for a long time :)

 
</div>

* Biopython [https://biopython.org/DIST/docs/tutorial/Tutorial.html tutorial]

* Biopython [https://github.com/biopython/biopython/blob/a09f38b468cbedaf25aec6cfed1475b0849011e6/Bio/Blast/NCBIWWW.py examples]
 

==Part 2.0: Create your first Python Notebook==

[https://colab.research.google.com/drive/1DpqOH8O_KK0rEohjv-20CCx2THwCq89h#scrollTo=_7tJ_S6wc-Os Click Here]

We are going to make a local copy of this notebook on our own drive.
* First you will need to log in or create a Google account.
* Make a copy of the Colab file in your own drive. Go to File > Save in Drive

[[File:save_a_copy_in_drive.png|center|300px|border]]

Now you can start modifying the file!

* First, we will install BioPython in our notebook. This will only take a few seconds, but note, that this installation is temporary. If your notebook gets idle and disconnected you will need to run this code cell again.

! pip install biopython

* Now import Blast tools from BioPython

from Bio.Blast import NCBIWWW
help(NCBIWWW.qblast)

===Random DNA sequences and BLASTN===

*Generate 10 random DNA sequence of length 25bp using Google Colab

[[Image:Office-notes-line_drawing.png|30px|left]]

We will now do a BLASTN search using the random sequence as query. Select the option "blastn" in the section "Program Selection". Choose "Nucleotide Collection (nr/nt)" as the search database.

'''VERY IMPORTANT''':
For this special situation where we BLAST small artificial sequences, we need to turn off some the automatics NCBI incorporate when short sequences are detected. Otherwise we'll not be able to see the intended results:

* Extend the "Algorithm parameters" section (see the screen shot below) in order to gain access to fine-tuning the options.
*# '''Deselect''' the "Automatically adjust parameters for short input sequences" option.
*# Set the E-value cut-off ("Expect threshold") to '''50'''

[[file:Blastn_cropped+circle.png‎|center|frame|'''Remember to adjust the BLAST settings''']]

* Inspect the results.

:'''QUESTION 2.1''':
:Answer the following small questions, and '''document your findings''' by pasting in examples of alignments / text snippets from the overview plots generated on google colab:

:* What is the typical length of the hits (the alignment length)?
:* What is the typical % identity?
:* In what range is the bit-scores ("max score")?
:** ''Notice: This is conceptually the same as the "alignment score" we have already met in the pairwise alignment exercise''.
:* What is the range of the E-values?

:'''QUESTION 2.2''':
:*What is the '''biological''' significance of these hits / is there any biological meaning?

===Random protein sequences and BLASTP===
Now it's time to work with a set of '''protein sequences''': Generate 10 peptide sequences of length 25aa using Google colab.

* '''Notice 1:''' The distribution of amino acids will be equal (5% prob) and this is different from true biological sequences - however this is not important for this first part of the exercise.
* '''Notice 2:''' Please recall from the lecture that the way <tt>BLASTP</tt> selects candidate sequences for full Smith-Waterman alignment is different from <tt>BLASTN</tt>. (<tt>BLASTN</tt> - a single short (11 bp +) perfect match hit is needed. <tt>BLASTP</tt> - a pair of "near match" hits of 3 aa within a 40 aa window is needed).

[[Image:Office-notes-line_drawing.png|30px|left]]

Choose blastp as the algorithm to use, the "NR" database (this is the protein version, consisting of translated CDS'es, UniProt etc).

'''VERY IMPORTANT''': We also need to tweak the parameters this time - in the "Algorithm Parameters" section select BLOSUM62 as the alignment matrix to use and set the "Expect threshold" to 1000 (default: 10).

* Perform the BLAST search.
* Inspect the results.

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 2.3''':
:''(Remember to '''document your answers''' in the same manner as Q2.2)''

:* What is the typical length of the alignment and do they contain gaps?
:* What is the range of E-values?
:* Try to inspect a few of the alignments in details ("+" means similar sequences) - do you find any that look plausible, if we for a moment ignore the length/E-value?
:* If we had used the default E-value cut-off of 10 would any hits have been found?

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 2.4''':
:* If we compare the result from BLAST'ing random DNA sequences to random Peptide sequences - which kind of search has the higher risk of returning false positives (results that appear plausible, maybe even significant, but are truly unrelated)?
:** Remember to take E-values into your consideration.

 

==Part 3: using BLAST to transfer functional information by finding homologs==
<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' The following parts can either be done on the google colab or in the webserver. If you choose the web server version remember to limit your search to "bacteria" (taxid: 2) in ALL of this section (PART 3) to make the BLAST searches run quicker. (The organisms we're looking for all belong to the "Bacteria" domain of life, so this restriction is OK).
</div>

===Homo-, Ortho- and Paralogs===

One of the most common ways to use BLAST as a tool is in the situation where you have a sequence of '''unknown function''', and want to '''find out which function it has'''. Since a large amount of sequence data has been gathered over the years, chances are that an '''evolutionarily related''' sequence with known function has already been identified. In general, such a related sequence is known as a "'''homolog'''".

Homo-, Ortho- and Paralogs:
* A '''Homolog''' is a general term that describes a sequence that is related by any evolutionary means.
* An '''Ortholog''' ("Ortho" = True) is a sequence that is "the same gene" in a different organism: The sequences shared a single common ancestor sequence, and have now diverged through speciation (e.g. the Alpha-globin gene in Human and Mouse).
* A '''Paralog''' arises due to a gene duplication within a species. For example, Alpha- and Beta-globin are paralogs.

[[File:Homo_Ortho_Para-log.gif|center|frame|''Image source: [http://www.thegreatgoodplace.com/tt/gwlee/126 gwLee's blog]'' ]]

Notice that in both cases it's possible to transfer information, for example, about gene family / protein domains.
We have already touched upon the comparison of (potentially) evolutionarily related sequences in the pairwise alignment exercise. However, this time we do not start with two sequences we assume are related, but instead, we start with a single sequence ("query sequence") which we will use to search the databases for homologs (we often informally speak of "BLAST hits", when discussing the sequences found).

 



===BLAST example 2===

The sequence is a '''DNA fragment''' from an unknown non-cultivatable microorganism. It was cloned and sequenced directly from DNA extracted from a soil-sample, and it goes by the poetic name "CLONE12". It was amplified using degenerated PCR primers that target the middle ("core cloning") of the sequence of '''a group of known enzymes'''. (I can guarantee this particular sequence is not in the BLAST databases, since I have cloned and sequenced it myself, and it has never been submitted to GenBank).

LOCUS CLONE12.DNA 609 BP DS-DNA UPDATED 06/14/98
DEFINITION UWGCG file capture
ACCESSION -
KEYWORDS -
SOURCE -
COMMENT Non-sequence data from original file:
BASE COUNT 174 A 116 C 162 G 157 T 0 OTHER
ORIGIN ?
clone12.dna Length: 609 Jun 13, 1998 - 03:39 PM Check: 6014 ..
1 AACGGGCACG GGACGCATGT AGCTGGAACA GTGGCAGCCG TAAATAATAA TGGTATCGGA
61 GTTGCCGGGG TTGCAGGAGG AAACGGCTCT ACCAATAGTG GAGCAAGGTT AATGTCCACA
121 CAAATTTTTA ATAGTGATGG GGATTATACA AATAGCGAAA CTCTTGTGTA CAGAGCCATT
181 GTTTATGGTG CAGATAACGG AGCTGTGATC TCGCAAAATA GCTGGGGTAG TCAGTCTCTG
241 ACTATTAAGG AGTTGCAGAA AGCTGCGATC GACTATTTCA TTGATTATGC AGGAATGGAC
301 GAAACAGGAG AAATACAGAC AGGCCCTATG AGGGGAGGTA TATTTATAGC TGCCGCCGGA
361 AACGATAACG TTTCCACTCC AAATATGCCT TCAGCTTATG AACGGGTTTT AGCTGTGGCC
421 TCAATGGGAC CAGATTTTAC TAAGGCAAGC TATAGCACTT TTGGAACATG GACTGATATT
481 ACTGCTCCTG GCGGAGATAT TGACAAATTT GATTTGTCAG AATACGGAGT TCTCAGCACT
541 TATGCCGATA ATTATTATGC TTATGGAGAG GGAACATCCA TGGCTTGTCC ACATGTCGCC
601 GGCGCCGCC
//

[[Image:Office-notes-line_drawing.png|30px|left]][[Image:Cogs_brain.png|right|150px]]
:'''QUESTION 3.1 (Long question - read all)''':
: ''Your task is now to find out '''what kind of enzyme''' this sequence is likely to encode, '''using the methods''' you have learned''.

'''INSTRUCTIONS''': You are free to write the combined answer to this question in a free-style essay-like fashion - just be sure to include the subquestions in your answers. In an exam situation you will need to put all the clues together yourself, reason about the tools/databases to use, and document your findings.

'''STEP 1 - cleaning up the sequence''':

The sequence is (more or less) in GenBank format and the NCBI BLAST server expects the input to be in FASTA format, or to be "raw" unformatted sequence.

* There are two solutions to this:
** Copy the sequence into a text-editor and manually create a FASTA file ("search and replace" and/or "rectangular selection" is useful for the reformatting). This is the most robust solution: it will always work. (Look at [[Plain text files and Geany|the Geany exercise]] for a reminder of how to do this).
** Hope the creators of the web-server you're using were kind enough to automatically remove non-DNA letters (paste in ONLY the DNA lines) - this turns out to be the case for both NCBI BLAST and VirtualRibosome, but it ''cannot be universally relied upon''.

'''Subquestion''': convert the sequence to FASTA format (manually, in Geany) and quote it in your report.

'''STEP 2 - thinking about the task''':

Consider the following before you start on solving this task:
* Based on the information given: is the sequence protein-coding?
* If it is, can you trust it will contain both a START and STOP codon?
* Do we know if the sequence is sense or anti-sense?
and think which consequences the answers to these questions should have for your choice of methods and parameters.

'''Subquestion''': Give a summary of your considerations.

'''STEP 3 - Performing the database search''':

''Significance'': We will put the criteria for significance at 1e-10 (remember: the higher the E-value, the worse the significance).

'''Subquestion''':

Cover the following in your answer:
* What tool(s) and database(s) will be relevant to use?
* Document the results from the different BLAST searches - what works and what does not work?
* You need to copy in small snippets of the BLAST results to document what you observe.
* '''In conclusion''': What kind of enzyme is CLONE12? Gather as much evidence as possible.

Blast Job ID (blastn): DU9XHV8H014

Blast Job ID (blastp): DU9T4K44014

 

==Part 4: BLAST'ing Genomes==
<div style="background-color: lavender; border: solid thin grey;">
:'''IMPORTANT:''' do '''NOT''' limit your search to "bacteria" here - now we are actively looking at organism specific queries.
</div>

So far we have been using BLAST to search in the big broad databases that covers at huge set of sequence from a large range of organisms. In this final part of the exercise we will be doing some more focused searches in smaller databases by targeting specific genomes.

Typically this will be useful if you have a gene of known function from one organism (say a cell-cycle controlling gene from Yeast, ''Saccharomyces cerevisiae'') and want to find the human homolog/ortholog to this gene (genes that control cell division are often involved in cancer).

When you have been performing the BLAST searches, you have probably already noticed, that's it possible to search specifically in the Human and Mouse genomes (these database only contains sequences from Human/Mouse). It's also possible to restrict the output from searches in the large databases (e.g. NR) to specific organisms.

A growing number of organisms have been fully sequenced, and the research teams responsible for a large scale genome project typically put up their own Web resources for accessing the data. For example the Yeast genome is principally hosted in the Saccharomyces Genome Database (SGD - www.yeastgenome.org) - it should be noted that SGD also offers BLAST as a means to search the database.



===Genome specific analysis of histones===
====SGD====
Let's do a small study of the relationship between the histones found in Yeast and in Human (evolutionary distance: ~1-1.5 billion years).

Look up the '''HTA2''' gene in SGD (http://www.yeastgenome.org - use the search box at the top of the page). Notice that a brief description about the function of the gene and its protein product is displayed (a huge amount of additional information can be found further down the page - much of it Yeast specific).

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.1''':
:What information is given about the relationship between this gene and the gene "HTA1"?

Browse the page and locate the link to the protein sequence. Save the sequence as a file, '''we'll need it in a moment'''.

====NCBI====


Now return to the NCBI blastp page. Set Database to "Reference proteins (refseq_protein)", and enter <tt>Saccharomyces cerevisiae</tt> in the Organism field (and accept the suggestion with taxid:4932).

Blast Job ID: FR6NDFMJ013

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.2''':
:''(Remember to document your answers)''
:* How many high-confidence hits do we get?
:* Do the hits make sense, from what you have read about HTA2 at the SGD webpage?
:'''Tip:''' click on the Gene links under Related Information (to the right of the alignments) to see the gene names for the protein hits.

The next step is to search the translated version of the human genome.

Do as before, still with Database set to "Reference proteins (refseq_protein)", just enter <tt>Human</tt> in the Organism field.

Blast Job ID: FR65XVA0016

[[Image:Office-notes-line_drawing.png|30px|left]]
:'''QUESTION 4.3''':
:* How many high-confidence hits (with E-value better than 10-10) are found? (Approximately)
:* What are all the high-confidence hits called?




==Concluding remarks==
Today we have been using BLAST to find a number of homologous genes (and protein-products). If we want to go even deeper into the analysis of the homologs, the next logical step would be to build a dataset of the full-length versions of the sequences we have found (not just the part found by the local alignment in BLAST).

A further analysis could consist of a series of pairwise alignments (for finding out what is similar/different between pairs of sequences) or a multiple alignment which could form the basis of establishing the evolutionary relationship between the entire set of sequences.

BLAST can also be used as way to build a dataset of sequences base on a known "seed" sequence. As we saw in the GenBank exercise, free-text searching in the GenBank can be difficult, and if we for instance wanted to build a dataset of variants of the insulin gene, an easiy way to go around this would be to BLAST the normal version of the insulin against the sequence database of choice, and pick the best matching hits from here.