4.1 Substitution matrices #
4.1.1 Scores as log-odds #
The toy +1/โ1 scheme treats every substitution as equal. Proteins do not. Leucine replaces isoleucine often among real homologues. Tryptophan rarely changes, and a conserved tryptophan is surprising because tryptophan is rare, so a tryptophan identity scores high. A substitution matrix stores log-odds: if qแตขโฑผ is how often residues i and j are paired in trusted alignments, and pแตข pโฑผ is the chance frequency of that pair, the score is proportional to log(qแตขโฑผ / (pแตข pโฑผ)). Positive means more common among homologues than by chance.
The scale reappears as ฮป in the E-value, so raw scores from different matrices are not comparable until they are bit scores. A matrix counted on close proteins punishes changes that are normal at long distance. The choice of matrix is a claim about that distance.
4.1.2 PAM #
PAM (accepted point mutation) matrices were built by Dayhoff and colleagues from closely related proteins, on the order of 85% identical or higher, where the alignment is secure and multiple hits at one site are rare. PAM1 is one accepted mutation per 100 residues: fixed and still visible, not every mutation that occurred in a cell.
Longer distances are extrapolated by raising the PAM1 probability matrix to the n-th power, then converting back to log-odds. Do not power the integer scores. PAM250 is that 250-step model. The number counts steps, not percent identity. Expected identity at PAM250 is classically around 20%, because of multiple hits and reversions. The extrapolation assumes one process at every site, which real proteins violate. A low PAM, such as PAM30, remains the sensible choice for close sequences and short peptides.
4.1.3 BLOSUM #
BLOSUM matrices, from Henikoff and Henikoff, count substitutions in conserved ungapped blocks, including blocks that are already distant. Sequences more identical than a threshold are clustered and counted once, so near-copies cannot outvote a divergent sequence. BLOSUM62 means that clustering threshold was 62%. It does not mean “use this when the sequences are 62% identical”. Pairs inside a cluster do not contribute, and pairs between clusters can lie far below 62%.
The number runs opposite to PAM. BLOSUM80 is harsher, because only very similar sequences count as separate observations. BLOSUM45 is more tolerant. On the comparisons that chose it, BLOSUM62 recovered distant homologues better than the PAM matrices tried against it. That result, not a theorem, is why it is the protein-BLAST default. Tryptophan matched to itself scores 11 and cysteine matched to itself scores 9, higher than identities of common residues, because those conservations are more surprising.
4.1.4 Choosing one #
| Situation | Usual choice | Why |
|---|---|---|
| First general protein search | BLOSUM62 | The usual default |
| Close homologues, or a short peptide | BLOSUM80 or BLOSUM90, or PAM30 | A mismatch is surprising, and a short alignment cannot absorb noise |
| Remote homologues | BLOSUM45 or PAM250 | Tolerates substitutions that occur at long distance |
A harsh matrix misses a distant homologue. A tolerant matrix on a short peptide reports chance as homology. Several matrices are several looks (Chapter 10). Do not compare a raw PAM score with a raw BLOSUM score.
4.2 BLAST #
The study of BLAST provides fundamental principles governing cellular architecture, molecular interactions, and physiological regulation.
Blast Sample Output
4.2.1 Seeds, then extension #
BLAST answers a local-alignment question without filling a full matrix for every subject. SmithโWaterman remains the right tool for one query against one protein. It is not the default against an archive of hundreds of millions of sequences.
The shortcut is a seed. Protein search classically uses words of length 3; nucleotide search uses longer words, classically 11. Protein words need not match exactly: neighbourhood words that score above a threshold T against a query word are indexed too. Extension grows a hit in both directions, and gapped BLAST often requires two nearby hits on one diagonal before a gapped local alignment. The reported alignment is that extension, not the word. FASTA is the other classical heuristic, and the warnings below apply to both.
4.2.2 A heuristic can miss the optimum #
BLAST is not guaranteed to return the optimal local alignment, nor every homologue an exhaustive search would have ranked as significant. A homologous region with no seed above T is invisible. Distant homologues, short motifs, and unusual residue patterns are the usual misses. Lowering T and shortening the word buys sensitivity and spends time.
“No hit” means this word length, matrix, database, and filter produced nothing above the threshold you set. PSI-BLAST builds a position-specific matrix from one round and searches again. It reaches further, and a non-homologue admitted in one round brings its relatives into the next. Inspect every round. Keeping only the p-value you like, after several searches, is multiple testing.
4.2.3 E-value and bit score #
The raw score S sums substitution scores and gap penalties. The E-value is how many alignments of at least this score a chance model expects for this query in a database of this size, E = K ยท m ยท n ยท e^(โฮปS). Here m and n are effective lengths of query and database. ฮป and K depend on the scoring system and the background, not on one hit. The formula describes the tail of high scores, and it is a poor description of a tiny peptide.
The bit score removes that scoring system: Sโฒ = (ฮปS โ ln K) / ln 2, so E = m ยท n ยท 2^(โSโฒ). Doubling the database doubles E and leaves Sโฒ unchanged. Lengthening the query also changes the search space. For one query, one database, and one scoring system, bit score and E-value rank hits the same way. Change the database and they do not: an alignment with E = 0.01 has E near 1 in a database a hundred times larger, at the same bit score. Quote the database and its date.
E is an expectation, not a probability. The chance of at least one such alignment is about p = 1 โ e^(โE). When E is much less than 1, p โ E. When E is 5, about five chance alignments are expected and p is about 0.99, not 5. A printed E of `0.0` is underflow. A tiny E is evidence against the random model, not the probability that the sequences are homologous.
4.2.4 Three mistakes #
The first is reading E as a p-value, or as the probability that the proteins are unrelated. Convert when E is not small, and call neither number a probability of homology.
The second is ignoring database size. A moderate bit score that would matter against one proteome is ordinary against the non-redundant archive, because m ยท n is huge. The top of a long noise list is still noise. The printed E-value is per query, so a hundred queries are a hundred chances the program was not told about.
The third is low-complexity sequence: a glutamine run, a proline run, a hydrophobic transmembrane peptide, AT-rich DNA. Composition alone matches unrelated sequences, and the random model assumes something nearer independent draws, so E is too small. SEG (protein) and DUST (nucleotide) mask those stretches before they can seed a hit. A hit that vanishes when the filter is on was carried by composition. Later composition adjustments do not replace reading the alignment. Percent identity without a length is not a result, and the product name on the subject is someone else’s annotation.
4.3 Key terms, summary, and self-check #
Key terms. Log-odds; PAM; BLOSUM; BLOSUM62; seed; neighbourhood word; heuristic; PSI-BLAST; E-value; bit score; KarlinโAltschul; low-complexity filter.
Summary. Scores are log-odds against chance. PAM counts accepted mutations at short distance and extrapolates by powering a probability matrix; a higher PAM number is a longer distance. BLOSUM counts substitutions in blocks after clustering; a higher BLOSUM number is a shorter distance, and BLOSUM62 is the default because it worked on the tests that chose it. BLAST seeds and extends, and can miss an alignment that never seeds. E scales with database size and query length. The bit score does not. E is not a probability of homology, and low-complexity sequence makes the null model false unless it is masked.
Self-check.
1. Why is PAM250 not a matrix “for 250% identity”, and BLOSUM62 not a matrix “for sequences that are 62% identical”?
2. For a 20-residue peptide from two strains of one species, which end of the BLOSUM series is the better first choice, and why?
3. BLAST reports no hit. Give two reasons still compatible with a real homologue.
4. Bit score 50 and E = 0.02. The database is taken ten times larger. What happens to each number?
5. A glutamine repeat matches dozens of proteins only when filtering is off. What was wrong with those E-values?
