5 Session 5: Advanced Database Search & Profile Methods
Sources for this page: Petras Kundrotas’s KB8029 lecture “Advanced Database Search” (Feb 2026) — this page follows its structure and adapts several of its worked examples (independently recomputed and verified, not copied), not copied verbatim; Altschul, Gish, Miller, Myers & Lipman (1990), J. Mol. Biol. 215:403-410 (BLAST); Altschul et al. (1997), Nucleic Acids Research 25:3389-3402 (PSI-BLAST and gapped BLAST); Karlin & Altschul (1990), PNAS 87:2264-2268 (the statistical foundation BLAST’s E-values are built on); Steinegger & Söding (2017), Nature Biotechnology 35:1026-1028 (MMseqs2); Korf, Yandell & Bedell, BLAST (O’Reilly, 2003), pp. 79-87 (the drop-off extension example, adapted); Wikipedia, BLAST (biotechnology), PSI-BLAST, Position weight matrix, Pseudocount, Hidden Markov model, Pfam, InterPro, E-values — CC BY-SA 4.0; NCBI, BLAST Help — public domain (US government work); EMBL-EBI’s InterPro — CC0. The HBB Pfam domain hit and all worked-example numbers below were computed or fetched live while writing this page (2026-09-18) — see the notebook.
5.1 Learning goals
By the end of this session you should be able to:
- explain why exact dynamic programming does not scale to searching a whole database;
- describe BLAST’s seed-extend-evaluate heuristic, and how MMseqs2 scales the same ideas;
- interpret the E-value and bit score of a real BLAST result;
- construct a small position-specific scoring matrix (PSSM) by hand, and explain why pseudocounts matter;
- describe PSI-BLAST as BLAST iterated with a PSSM instead of a single sequence, and explain why iteration can drift;
- explain, at a conceptual level, what a profile Hidden Markov Model adds beyond a PSSM.
5.2 The scaling problem, and BLAST’s heuristic answer
Session 4’s dynamic programming finds the provably optimal alignment between two sequences, but its cost is proportional to the product of their lengths. For two sequences a few hundred residues long, this is no problem. For searching one query against a database, it is. GenBank and UniProt hold well over a hundred million sequences between them, and running full dynamic programming against every entry would take far too long to be useful. A database search therefore needs a faster method.
BLAST (Basic Local Alignment Search Tool) is such a method, and it trades the guarantee of optimality for speed. Instead of filling a DP matrix for every database sequence, it first looks for short exact or near-exact matches. Only the promising regions that these matches point to get the expensive alignment step. Concretely, BLAST runs in three stages:
- Seeding. Break the query into short “words”, list every word that scores well against each of them, and scan the database for occurrences of those words.
- Extension. Grow each seed hit outwards, in both directions, into a full local alignment.
- Evaluation. Ask whether the score of the resulting alignment is better than what pure chance would produce.
Each stage is covered in its own section below. This three-stage picture is the core idea, and it is deliberately simplified. Production BLAST adds refinements such as the two-hit rule (only extend when two seeds lie close together on the same diagonal), masking of low-complexity regions, gapped extension, and corrections for unusual amino-acid composition. None of these refinements changes the logic of seed, extend, evaluate.
5.2.1 Further reading
- Wikipedia, BLAST (biotechnology)
- NCBI, BLAST Help — the practical reference for choosing blastp/blastn/blastx and reading real output
5.3 Stage 1: seeding — turning the query into fast lookups
The first stage turns the query into a set of short words. For protein sequences, a “word” is usually only 3 residues long, so a query of length 200 has 198 overlapping 3-letter words. BLAST does not search the database only for exact matches to these words, because a biologically real match may substitute a residue or two without losing significance. Instead, it uses a substitution matrix (BLOSUM62, from Session 4) to precompute every 3-letter combination that scores above a threshold \(T\) against each query word.
Worked example. Scored with the real BLOSUM62 matrix, several other words besides the query word LWG itself clear a typical threshold of \(T=12\). The notebook computes them directly from BLOSUM62, and 49 words clear the threshold in total. The top ones are:
LWG.
| Neighbour word | Score vs. LWG |
|---|---|
LWG |
21 |
IWG |
19 |
MWG |
19 |
VWG |
18 |
FWG |
17 |
AWG |
16 |
CWG |
16 |
TWG |
16 |
YWG |
16 |
KWG |
15 |
These query words and their high-scoring neighbours go into a lookup table, which is built once per query and records, for every 3-letter word, which query positions it would seed. BLAST then scans each database sequence once, from start to end. For every 3-letter word it reads, it asks the table “does this word seed anything in the query, and where?” Because that check is a single table lookup, the scan is fast and linear in the size of the database. Only the few places that produce a seed go on to the expensive extension step.
5.4 Stage 2: extension — growing a seed into a real alignment
A seed hit is only a 3-letter coincidence and not yet evidence of anything. BLAST therefore extends it outwards in both directions and recomputes a running score as it goes. The scoring uses the same substitution-matrix logic as Session 4’s dynamic programming, but without filling a full matrix. Instead, extension stops as soon as the running score has fallen too far below its own best point so far.
The X-drop rule. BLAST tracks the running score as the alignment extends, and separately the best score seen so far. If the running score ever falls more than a threshold \(X\) below that best score, the alignment is very unlikely to recover, so extension stops. The reported alignment is then trimmed back to the point where the best score occurred.
Worked example, adapted from Korf, Yandell & Bedell (2003). We treat two short strings as if they were sequences, score \(+1\) for a matching character and \(-1\) for a mismatch (ignoring spaces), and use \(X=5\):
query: Thequickbrownfoxjumps...
database: Thequietbrownncatpurrs...The running score climbs steadily through Thequ, dips slightly at the ick/iet mismatch, and then recovers and climbs further through brown to a peak of 9 right after ...brown/...brown. After that it falls steadily, as foxjump diverges from catpurr. Once the running score has dropped 5 points below the peak of 9, extension stops. The reported alignment is trimmed back to the peak itself, so it pairs Thequickbrown with Thequietbrown and leaves out the lower-scoring tail beyond it. The notebook reproduces this computation in code.
5.4.1 Further watching
BLAST 2: seed and extend (Rob Edwards, San Diego State University):
5.5 Stage 3: is this hit significant? E-values and bit scores
Extending every seed produces a raw alignment score. A raw score alone, however, does not say whether the alignment is meaningful or simply what we would expect from two random sequences of that length. To decide, we need to know how scores are distributed by chance. Karlin and Altschul proved in 1990 that the scores from aligning a query against random sequence data follow an extreme value (Gumbel) distribution, not a normal distribution, and this result underpins BLAST’s statistics.
From that distribution follows the expected number of alignments that score at least \(S\) purely by chance, the E-value:
\[ E = K m n\, e^{-\lambda S} \]
Here \(m\) and \(n\) are the lengths of the query and the database, and \(K\) and \(\lambda\) are constants of the particular scoring system (substitution matrix and gap penalties), which real BLAST output reports as “Karlin-Altschul statistics”. This is the simplified textbook form. Real BLAST replaces \(m\) and \(n\) by slightly smaller effective lengths (an alignment cannot start right at the end of a sequence) and adjusts for composition, but the behaviour is the same.
Two consequences of the formula matter. First, a higher score gives an exponentially smaller \(E\). Second, a bigger database (\(n\)) inflates \(E\) for the same alignment, because chance has more opportunities to produce a false hit. An E-value is therefore only meaningful together with the database it was computed against, and the same alignment can be “significant” against a small database but not against all of GenBank.
Raw scores and E-values both depend on the scoring system, which makes them awkward to compare across different searches. The bit score solves this by normalising away the scoring system:
\[ S' = \frac{\lambda S - \ln K}{\ln 2} \]
Bit scores from completely different searches, even ones that use different substitution matrices, can therefore be compared directly. A bit score is not a significance test on its own, though. There is no universal bit-score cutoff, because whether a given score is surprising still depends on the size of the search space it was found in. For significance, read the E-value.
Reading real output. A real local BLASTP search of HBB (beta-globin, P68871) against its five closest human globin relatives gives the following hits:
| Hit | % identity | E-value | Bit score |
|---|---|---|---|
Delta-globin (HBD) |
93.2% | 9.3×10⁻¹⁰⁵ | 284 |
Epsilon-globin (HBE) |
75.5% | 2.1×10⁻⁸⁷ | 240 |
Alpha-globin (HBA) |
43.4% | 9.0×10⁻³⁸ | 114 |
Zeta-globin (HBAZ) |
35.9% | 3.8×10⁻³² | 100 |
Myoglobin (MYG) |
16.0% | 63 | 10.8 |
The first four are all confident, high-quality hits, with high percent identity, astronomically small E-values and comfortably high bit scores. Myoglobin is the interesting case. It shares the globin fold and a real evolutionary origin with haemoglobin, yet BLAST reports only a weak alignment with an E-value of 63, which taken at face value is no statistical evidence at all. That number, however, only evaluates the alignment BLAST found. If we run the same search with the two-hit rule switched off (so that every single seed gets extended; -window_size 0 in BLAST+), BLAST finds a full-length alignment, 145 residues at 25.5% identity, with an E-value of 5×10⁻¹¹. With the default settings, no pair of seeds happened to fall close enough together on the right diagonal, so the good alignment was never extended. This is the price of a heuristic. BLAST is fast because it skips most of the search, and occasionally it skips the answer, which an optimal Smith–Waterman search would not have missed.
Database size still matters as well. These E-values are relative to a tiny five-sequence database. Against a database the size of UniProtKB (about \(9\times10^{10}\) residues), the same full-length myoglobin alignment would get an E-value of about 0.007, which is significant, but only just. Relatives that have diverged further than myoglobin, such as the bacterial truncated haemoglobins, give no detectable pairwise signal at all, even with the optimal alignment (Session 6 shows that they still share the globin fold).
5.5.1 Further reading
- Wikipedia, E-values
- Original source: Karlin & Altschul (1990), PNAS 87:2264-2268 — the full derivation behind the formula above
5.6 BLAST’s ideas at modern database scale: MMseqs2
BLAST was designed in 1990, when the protein databases held tens of thousands of sequences. Today’s metagenomic collections hold billions, and at that scale even BLAST becomes too slow. MMseqs2 (Steinegger & Söding, 2017) is a widely used modern search tool that keeps BLAST’s logic but turns one part of it around. In simplified form, it works as follows:
- Index the target database by short words (k-mers, typically 6–7 residues for proteins). The index is a table that maps each k-mer to the database sequences and positions where it occurs, and it is built once and reused.
- For each k-mer in the query, generate the similar k-mers (scored with a substitution matrix, exactly like BLAST’s neighbour words) and look them up in the index.
- Keep only database sequences in which two matching k-mers fall on the same diagonal, that is, at the same offset between query and target position, just as in a dot plot from Session 4. Two consistent hits are much stronger evidence of a real alignment than one.
- Score the surviving candidates quickly with an ungapped alignment along that diagonal, and discard the weak ones.
- Run the expensive gapped Smith–Waterman alignment, heavily optimised to use the processor’s vector instructions, only on the small set that is left.
The difference from classic BLAST lies in which side gets indexed. BLAST builds a lookup table from the query and its neighbour words, and then reads through every database sequence. MMseqs2 instead indexes the database once, so each query only has to look up its own k-mers. Steps 3–5 form a cascade of ever more expensive filters, each of which passes on only a small fraction of the candidates. As a result, MMseqs2 is often hundreds of times faster than BLAST at similar sensitivity. No machine learning is involved, because MMseqs2 simply engineers the same seed-extend-evaluate ideas for scale.
Like PSI-BLAST below, MMseqs2 can also iterate, turning its hits into a profile and searching again, and ColabFold uses it to build the alignments that AlphaFold needs (Session 12). Session 6 meets its structural counterpart, Foldseek, from the same group.
5.7 From one sequence to many: PSI-BLAST and the PSSM
Myoglobin sits near the edge of what a single sequence can detect, and more distant relatives lie beyond it. The reason is that ordinary BLAST scores everything against one fixed substitution matrix, whatever family the query belongs to, so a faint but real family signal is easily lost against the background of a large database. PSI-BLAST (Position-Specific Iterated BLAST) does better, because it builds a query-specific scoring matrix and refines it over several rounds:
- Run an ordinary BLAST search with the query sequence.
- Take the hits that pass the inclusion threshold (an E-value cutoff, by default stricter than the one used for reporting hits), align them to the query, and turn that alignment into a position-specific scoring matrix (PSSM), described below. Hits that are displayed but fall short of the inclusion threshold do not enter the profile.
- Search the database again, this time scoring with the PSSM instead of a generic substitution matrix.
- Repeat, so that each round’s included hits refine the PSSM further, until no new sequences pass the threshold (convergence) or a set number of rounds is reached.
This makes PSI-BLAST more sensitive than a single BLAST pass, particularly for remote homologs whose similarity to the original query is too faint for a fixed substitution matrix to catch.
Iteration is not automatically better: profile drift. Suppose that one unrelated sequence passes the inclusion threshold. It then enters the alignment and shifts the next PSSM a little towards itself. The next round finds more sequences like the intruder, which shift the profile further. After a few rounds, the “family” can be dominated by proteins that have nothing to do with the original query, often through a shared low-complexity or coiled-coil region rather than real homology. This is a general lesson that returns in the machine-learning half of the course. A model is only as good as the examples it was built from, and a model that chooses its own training data can reinforce its own mistakes. In practice, keep the inclusion threshold strict, look at what entered in each round, and check surprising late-round hits independently.
5.7.1 Building a PSSM, step by step
Session 4 ended with alignment columns, their entropy, and sequence logos. A logo describes an alignment column, whereas a PSSM turns the same column statistics into numbers that can score a new sequence. The recipe below is a teaching model. The real PSI-BLAST computation also weights the sequences and uses more sophisticated pseudocounts, as noted along the way, but every step serves the same purpose in the real thing.
1. Count. From a multiple sequence alignment, count how often each residue appears at each position. Take this small toy alignment of 5 sequences:
Position: 1 2 3 4 5
Protein 1: A G L S P
Protein 2: A G L T P
Protein 3: R G I S P
Protein 4: A A L S Q
Protein 5: A A L S PAt position 1, 4 of the 5 sequences have A and 1 has R. At position 5, 4 have P and 1 has Q, and so on.
Real alignments are rarely this balanced. If a search returns fifty nearly identical mammalian sequences and three distant ones, raw counts would make the profile describe mammals rather than the family. Real tools therefore reduce redundancy or weight sequences, removing or down-weighting near-duplicates so that an over-represented subfamily does not dominate the counts. The exact recipe varies between programs, but the principle is what matters.
2. Add pseudocounts. Turning raw counts straight into frequencies causes two problems. First, a residue that never appears at a position gets a frequency of exactly zero, and its log-odds score becomes \(-\infty\). Second, with only 5 sequences, the absence of a residue might be bad luck rather than a real biological signal. Pseudocounts, small constants added to every count before the frequencies are computed, solve both problems. In our teaching model we add \(+1\) uniformly to the count of every residue type at every position. Conceptually, a pseudocount is a prior, that is, what we assume about a column before looking at the alignment (here, “every residue is possible”). With only five sequences the prior has a large influence, while with hundreds of sequences the observed counts swamp it. Real tools use smarter priors, for example by giving larger pseudocounts to residues that substitute easily for the observed ones (according to a substitution matrix).
3. Compare to background. Compute each residue’s background frequency \(q_j\), that is, how often it occurs in proteins generally (e.g. alanine ≈7.4%, tryptophan ≈1.3%), independently of this alignment.
4. Compute log-odds scores. For each position \(i\) and residue \(j\),
\[ S_{i,j} = \log_2 \frac{f_{i,j}}{q_j} \]
Here \(f_{i,j}\) is the (pseudocount-adjusted) observed frequency. With base 2, a score of \(+1\) means “twice as likely here as by chance”, \(+3\) means “8× as likely”, and so on. Arranged into a matrix, with one row per alignment position and one column per residue, the scores \(S_{i,j}\) are the PSSM. (Some tools store the matrix the other way round, but the orientation is only a storage convention.) Scoring a candidate sequence against the PSSM means looking up one value per position and summing them. This is the weighted-sum, or dot-product, idea from Session 1, generalised to sequence positions.
The notebook builds this PSSM from the alignment above and confirms that the highest-scoring sequence the matrix can produce is AGLSP, the alignment’s own consensus, as it should be.
Why pseudocounts really matter. With only five sequences, every frequency is a small-sample estimate, and the raw log-odds score takes it at face value. Take a position where 4 of 5 sequences show alanine (background 7.4%) and 1 shows the much rarer tryptophan (1.3%). The raw scores are \(\log_2(0.8/0.074) \approx 3.43\) for alanine and \(\log_2(0.2/0.013) \approx 3.94\) for tryptophan, so a single observation outscores four. Yet one tryptophan in five sequences could easily be chance, and the score treats it as a reliable frequency of 20%. With \(+1\) pseudocounts the frequencies become 5/25 and 2/25, and the scores fall to 1.43 and 2.62. The pseudocounts make every score more cautious, but they do not change the ranking: tryptophan still wins. A pseudocount cannot tell a lucky single observation from a real one; only more sequences can.
Gaps. A PSSM contains only substitution scores, one row per position. To align a new sequence to it, PSI-BLAST therefore uses ordinary gapped alignment, the same dynamic programming as in Session 4, with the PSSM row for position \(i\) taking the place of the substitution-matrix row for query residue \(i\). The gap-open and gap-extension penalties are the ordinary ones, the same at every position.
A PSSM is already a linear model. For an ungapped sequence \(x = x_1 x_2 \dots x_L\), the PSSM score is
\[ S(x) = \sum_{i=1}^{L} w[i, x_i], \]
where \(w[i, x_i]\) is the matrix entry for residue \(x_i\) at position \(i\). The score is thus one weight looked up per position and added up. Session 7 keeps this weighted-sum idea but learns the weights from labelled examples instead of deriving them from alignment counts.
5.7.2 Further reading
- Wikipedia, PSI-BLAST, Position weight matrix, and Pseudocount
5.9 Worked example: a real profile-HMM hit on HBB
The same idea works on a real sequence at production scale. To tie this back to the HBB (beta-globin) locus from Session 3, we search Pfam’s library of profile HMMs against the real UniProt entry P68871 (HBB_HUMAN), fetched live via EBI’s InterPro API. The search finds a match to PF00042 (“Globin”), spanning residues 27-142 of the 147-residue protein, with an E-value of around \(1.3\times10^{-29}\). In a database of this size, essentially no match that good is expected by chance. The profile HMM thus does on a real sequence what the toy PSSM above did on paper: it scores how well a sequence matches a family, position by position.
5.10 The course notebook
notebooks/session05-blast-and-profiles-in-code.ipynb reproduces every worked example on this page in code. It traces the X-drop extension, builds the PSSM from the toy alignment (with and without pseudocounts, including the rare-residue-dominance check), and fetches HBB’s real Pfam domain hit live from the InterPro API. Open it in Google Colab via the badge at the top of the notebook page to edit and run it yourself.
5.11 Check your understanding
The questions below give immediate feedback and cover the sections above. They are in the same style as the course’s live in-class peer-instruction quizzes, and the full question bank is in quizzes/session05.md in the repository.
5.12 The lab
The full lab covers BLAST and profile methods as one continuous exercise. (Session 4 now has its own short lab for alignment practice, so this one no longer needs to cover alignment as well.) You first identify a real, unknown gene from a raw sequence through a real NCBI BLAST search. It turns out to be a sixth real human globin-family gene, continuing this course’s running thread. Three code exercises follow: a real Clustal Omega alignment and PSSM on five globin-family proteins, a real iterated profile search (jackhmmer, the profile-HMM counterpart of PSI-BLAST) that compares its first round with its third against every sequence in the PDB, and a live Pfam domain search for HIV-1’s Gag-Pol polyprotein via the same InterPro API that this page’s HBB worked example uses.
Every search in notebooks/session05-lab.ipynb runs on a server at EMBL-EBI, reached through its web API: the notebook sends the sequences, the server runs the program, and the notebook reads the result back. There is nothing to install, so the lab runs on Colab, in the browser or from VS Code. The same pattern (submit a job, wait for it, fetch the result) works for most large bioinformatics services, and it is how you will usually run big searches from your own code.
5.12.1 Submitting this lab
This lab is administered through Canvas, not through this book. See the course’s Canvas page for the current quiz, the deadline and how to submit.
5.13 What’s next
A remote homolog can become impossible to recognise from its sequence while its fold stays conserved. Session 6 therefore turns to structure: the four levels of protein structure, explored hands-on in PyMOL, structure alignment, and structural search with Foldseek. Session 7 then returns to the weighted sum \(S(x) = \sum_i w[i, x_i]\) that a PSSM computes. It formalises the sum as a perceptron and shows that its weights can be learned from labelled examples by gradient descent, instead of being derived from alignment statistics as this session’s PSSM was.