Suppose that we are walking along a street and see something that we have never seen before. As good scientists, we decide that we would like to know what this unkown object is and what it does. We could do several things. We could dismantle it and try and see what happens, or we could begin by looking at it and asking does it look like anything that we already know what the function is. The second option is obviously easier, and quicker in the long run.
Note this was inspired in part by The Salvation of Doug. I encourage everyone to read it.
That is one of the most important aspects of the BLAST search. If we come across a DNA or protein sequence and we don't know what it is, we can compare it to all the sequences that we do know and try and identify something similar where the function has already been described. This would be a good starting point to track down the actual function of the DNA or protein that we are interested in.
BLAST searches take DNA sequences or protein sequences and use some very fast computer tricks to find optimal alignments of these sequences to a database. The alignment process involves some very complicated statistics for identifying High-scoring Segment Pairs, regions of sequence (segments) that align together (pairs) with a probability that is greater than a random chance (high-scoring)
The results are then displayed in several forms. There is a graphical overview, a list of the most significant HSPs with their score and probability, and a view of the alignments of those high-scoring segment pairs
For this course, we are not going to consider the details of the BLAST system. We are going to focus on understanding how to use the system, understanding the results, and what the results actually mean.
Statistics
There is a lot of information about the BLAST search scoring system from NCBI, the Nucleic Acids Research paper describing BLAST, linked from the BLAST home page, and in the references.
For this course, we are not going to place a lot of emphasis on theory behind how the statistics are generated, however it is important that you know what they mean and understand the mechanisms used to generate them as this will affect how you view the results and the results that you get.
The score is the basis of the BLAST search. Each nucleotide or amino acid in an alignment is given a score depending on whether it matches or not. If there is no match the score is negative, and if there is a match the score is positive. For each region of alignment all the scores are added up, though the score can never go below zero. So, for example, consider the following alignment:
position 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
A T T T C C G C C C G A T G G C G C A G T C C A C
T T C G C C C A G T
C G T T C A G C C C A T A A A A C T A G T T G T G
score -1 -1 1 1 1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1
running score 0 0 1 2 3 2 3 4 5 6 5 4 3 2 1 0 0 0 1 2 3 2 1 0 0
In this simple alignment the first positive running score value is at position 3. This is therefore defined as the start of an high scoring region. The maximum running score is at position 10 (maximum score = 6) and so this is defined as the potential end of the high scoring region. The value would turn negative at position 17, and so this confirms the end of the high scoring region. The next region that is positive is at position 19 and runs through position 21. Again, the result would turn negative at position 25, which is also the end of the sequence, and so 19-21 is also a high scoring region.Although I used a simple +1/-1 system to demonstrate the scoring principle, in reality the BLAST program running a nucleotide search by default will use a +1 score for a positive match and a -3 score for a mismatch. This score is then normalized based on the scoring system used and whether gaps were included.
DNA alignments are simple because there are only limited changes that can take place, and so we can score a match with a positive value and a mismatch with a negative value. However with protein sequences the changes that can occur are more limited. For example, consider the following amino acids: (note the full genetic code is available in the Glossary)
Serine: AGT or AGC
Arginine: AGA or AGG
Phenylalanine: TTT or TTC
To change from a serine to an arginine will only take one base change (the last base) for either of the serine residues. To change from a serine to a phenylalanine will take two base changes (the first two bases). Therefore a Ser->Phe change is much less likely to happen than a Ser->Arg change. Of course, there are several other factors like function and protein structure that have to be taken into account. If you are aligning two protein sequences you should give a Ser->Phe change more penalty than a Ser->Arg change because there will have been more changes at the DNA level (therefore less chance of an amino acid substitution happening).
Dayhoff and others studied how frequently amino acids change to other observed aino acids in closely related proteins. From these results they constructed the PAM (for "point accepted mutation") model for amino acid substitutions to give each amino acid a score for the likelihood that it will change to another amino acid.
For example, a PAM matrix may look like the table below. Positive scores are highlighted for easy identification, and the scores for Ser->Arg and Ser-Phe are shown in red
PAM100 matrix
| Amino Acid | A | R | N | D | C | Q | E | G | H | I | L | K | M | F | P | S | T | W | Y | V | B | Z | X | * |
| A | 4 | -3 | -1 | -1 | -3 | -2 | 0 | 1 | -3 | -2 | -3 | -3 | -2 | -5 | 1 | 1 | 1 | -7 | -4 | 0 | -1 | -1 | -1 | -9 |
| R | -3 | 7 | -2 | -4 | -5 | 1 | -3 | -5 | 1 | -3 | -5 | 2 | -1 | -6 | -1 | -1 | -3 | 1 | -6 | -4 | -3 | -1 | -2 | -9 |
| N | -1 | -2 | 5 | 3 | -5 | -1 | 1 | -1 | 2 | -3 | -4 | 1 | -4 | -5 | -2 | 1 | 0 | -5 | -2 | -3 | 4 | 0 | -1 | -9 |
| D | -1 | -4 | 3 | 5 | -7 | 0 | 4 | -1 | -1 | -4 | -6 | -1 | -5 | -8 | -3 | -1 | -2 | -9 | -6 | -4 | 4 | 3 | -2 | -9 |
| C | -3 | -5 | -5 | -7 | 9 | -8 | -8 | -5 | -4 | -3 | -8 | -8 | -7 | -7 | -4 | -1 | -4 | -9 | -1 | -3 | -6 | -8 | -5 | -9 |
| Q | -2 | 1 | -1 | 0 | -8 | 6 | 2 | -3 | 3 | -4 | -2 | 0 | -2 | -7 | -1 | -2 | -2 | -7 | -6 | -3 | 0 | 5 | -2 | -9 |
| E | 0 | -3 | 1 | 4 | -8 | 2 | 5 | -1 | -1 | -3 | -5 | -1 | -4 | -8 | -2 | -1 | -2 | -9 | -5 | -3 | 3 | 4 | -2 | -9 |
| G | 1 | -5 | -1 | -1 | -5 | -3 | -1 | 5 | -4 | -5 | -6 | -3 | -4 | -6 | -2 | 0 | -2 | -9 | -7 | -3 | -1 | -2 | -2 | -9 |
| H | -3 | 1 | 2 | -1 | -4 | 3 | -1 | -4 | 7 | -4 | -3 | -2 | -4 | -3 | -1 | -2 | -3 | -4 | -1 | -3 | 1 | 1 | -2 | -9 |
| I | -2 | -3 | -3 | -4 | -3 | -4 | -3 | -5 | -4 | 6 | 1 | -3 | 1 | 0 | -4 | -3 | 0 | -7 | -3 | 3 | -3 | -3 | -2 | -9 |
| L | -3 | -5 | -4 | -6 | -8 | -2 | -5 | -6 | -3 | 1 | 6 | -4 | 3 | 0 | -4 | -4 | -3 | -3 | -3 | 0 | -5 | -4 | -3 | -9 |
| K | -3 | 2 | 1 | -1 | -8 | 0 | -1 | -3 | -2 | -3 | -4 | 5 | 0 | -7 | -3 | -1 | -1 | -6 | -6 | -4 | 0 | -1 | -2 | -9 |
| M | -2 | -1 | -4 | -5 | -7 | -2 | -4 | -4 | -4 | 1 | 3 | 0 | 9 | -1 | -4 | -3 | -1 | -6 | -5 | 1 | -4 | -2 | -2 | -9 |
| F | -5 | -6 | -5 | -8 | -7 | -7 | -8 | -6 | -3 | 0 | 0 | -7 | -1 | 8 | -6 | -4 | -5 | -1 | 4 | -3 | -6 | -7 | -4 | -9 |
| P | 1 | -1 | -2 | -3 | -4 | -1 | -2 | -2 | -1 | -4 | -4 | -3 | -4 | -6 | 7 | 0 | -1 | -7 | -7 | -3 | -3 | -1 | -2 | -9 |
| S | 1 | -1 | 1 | -1 | -1 | -2 | -1 | 0 | -2 | -3 | -4 | -1 | -3 | -4 | 0 | 4 | 2 | -3 | -4 | -2 | 0 | -2 | -1 | -9 |
| T | 1 | -3 | 0 | -2 | -4 | -2 | -2 | -2 | -3 | 0 | -3 | -1 | -1 | -5 | -1 | 2 | 5 | -7 | -4 | 0 | -1 | -2 | -1 | -9 |
| W | -7 | 1 | -5 | -9 | -9 | -7 | -9 | -9 | -4 | -7 | -3 | -6 | -6 | -1 | -7 | -3 | -7 | 12 | -2 | -9 | -6 | -8 | -6 | -9 |
| Y | -4 | -6 | -2 | -6 | -1 | -6 | -5 | -7 | -1 | -3 | -3 | -6 | -5 | 4 | -7 | -4 | -4 | -2 | 9 | -4 | -4 | -6 | -4 | -9 |
| V | 0 | -4 | -3 | -4 | -3 | -3 | -3 | -3 | -3 | 3 | 0 | -4 | 1 | -3 | -3 | -2 | 0 | -9 | -4 | 5 | -4 | -3 | -2 | -9 |
| B | -1 | -3 | 4 | 4 | -6 | 0 | 3 | -1 | 1 | -3 | -5 | 0 | -4 | -6 | -3 | 0 | -1 | -6 | -4 | -4 | 4 | 2 | -2 | -9 |
| Z | -1 | -1 | 0 | 3 | -8 | 5 | 4 | -2 | 1 | -3 | -4 | -1 | -2 | -7 | -1 | -2 | -2 | -8 | -6 | -3 | 2 | 5 | -2 | -9 |
| X | -1 | -2 | -1 | -2 | -5 | -2 | -2 | -2 | -2 | -2 | -3 | -2 | -2 | -4 | -2 | -1 | -1 | -6 | -4 | -2 | -2 | -2 | -2 | -9 |
| * | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | -9 | 1 |
Using the Phe/Ser/Arg example above, what are the relative scores for Ser->Arg and Ser->Phe substitutions? Does this make sense?
There are a lot of different matrices that have been optimized for different alignment parameters. An alternative to PAM is the BLOSUM matrix which was derived using blocks of similar protein sequences.
In essence, the critical parameter for protein searches is the length of the query sequence. As a rough guide the following should be used as each has been optimized for these conditions
| Query length | Substitution matrix |
| <35 amino acids | PAM-30 |
| 35-50 amino acids | PAM-70 |
| 50-85 amino acids | BLOSUM-80 |
| >85 amino acids | BLOSUM-62 |
The latest BLAST searches can introduce gaps into a sequence, and these are included as negative scores in the match. The default values for a gap is -5 to start a gap and -2 to continue a gap once it has been started. Gap scores are calculated based on the length of the gap, so a gap score would be:
-(opening cost + (length of gap x continuing cost)),
-(5 + (8 x 2)) = -21
Therefore there is a heavy penalty for opening and continuing a gap in a run of bases.
When you get a BLAST result you also get a probability value [reported as P (probability) or E (expected occurence of random sequences)]. This value is the probability that the high-scoring segment pair occurs by random. For example if you take a completely random sequence the same length as your query sequence there is a chance that you could get exactly the same sequence by random. This probability is based on the length of the query sequence and the total length of the database. If you are comparing a short sequence it is more likely that a random sequence could give you the same result. If you are comparing a sequence to a large database of sequences there is more chance that you have a random sequence in there that matches your query sequence than if you are comparing your sequence to a small database.
A probability of 0 means that there is essential no chance your match was random. A probability of 10-50 (reported as 1 E -50) means that 1 result in 1050 results reported will be random and not significant. This is not very likely, and so most hits returned with a score of 10-50 are probably real. In contrast a result of 0.1 means that 1 in 10 of the results that you have got back occurred by random and are not significant. It is up to you to figure out which is the random result!
The National Centers used to maintain a single database that contained all the sequence data. A lot of effort went into maintaining that database (imaginatively called the non-redundant dtabase) and removing redundant sequences from the database. A submission to any of the centers would result in the permeation of the data into all the databases. Therefore people in Europe could submit to EMBL; people in Japan to the DNA databank of Japan; and people in the Americas to the NCBI. This is still true to some extent, and data is shared among all these systems, however because of the overwhelming amount of DNA sequence being made available (for example see this page on NCBI's growth statistics it is now not feasible to maintain separate databases. Instead there are a large, and ever growing, number of databases that you can search against. These include:
Query sequences can be entered in two formats. Most often you will use the fasta format. In this case it should have an ID line that begins with a greater than sign (>) and then the sequence should start on the next line like this:
>sequence ID 1 tatatatattatttcttcgcgatcgctagctagctagcatgcatgcatgcatcgataggaggatcgcgagatcgatcagNote that the sequence can be in upper or lower case (or both mixed together!). If the sequence contains numbers or spaces these are ignored.
Alternatively you can enter an accession number of GI (genbank identifier) number. For example, for the S. Enteritidis SefR sequence the accession number would be AF239978. This is the sequence that has been used in all of the following examples.
These are the basics behind the BLAST system. We are not going to discuss the actual algorithims, a lot more information is available on this from the appropriate papers (see references below) than I can go into here. Instead we are going to focus on how to use BLAST and what the results mean.
A lot of information was taken from Gary Olsen's bioinformatics course taught at the University of Illinois, Urbana-Champaign, (UIUC) and I am indebted to him for his help, advice, and suggestions.
BLAST searches
Matrices
More references are available from the NCBI site.