Past exam of the mathematical and Computational Biology course of the University of Cambridge / 2024 / Paper 1 / C4

A BLOSUM substitution matrix is estimated from conserved, ungapped blocks of related protein sequences. Closely related sequences are clustered to reduce overrepresentation, and aligned amino acid pairs from different clusters are counted. Let $q_{ij}$ be their normalized observed pair frequency, and $e_{ij}$ the frequency expected from independent draws at the background amino acid frequencies. With unordered pairs, $e_{ii}=p_i^2$ and $e_{ij}=2p_i{p}_j$ for $i\ne j$.

For BLOSUM62, scores are rounded half-bit logarithms of these odds:A positive entry means that the aligned pair occurs more often in related sequences than the independent-background model predicts; a negative entry means less often. Adding these scores accumulates evidence across aligned positions. A score is neither a percentage identity nor a direct substitution probability.^[ref]

A diagonal substitution matrix score measures enrichment over a background match, rather than simply rewarding an identical letter. tryptophan residues are comparatively rare, so a W–W match has a small chance background frequency $p_W^2$. Its observed conservation produces a much larger ratio $q_{WW}/p_W^2$ than the corresponding ratio for the more common leucine. The BLOSUM62 diagonal therefore awards W–W more than L–L. Both the background frequency and the observed pair frequency matter; rarity alone does not determine a diagonal score.^[ref]

The number specifies the sequence-identity threshold used when clustering related sequences before deriving a BLOSUM substitution matrix. At 62%, sufficiently similar sequences are grouped and their within-cluster comparisons are excluded. Lowering the threshold merges more sequences, emphasizing comparisons across greater divergence.

Consequently BLOSUM80 is suited to closer relatives, BLOSUM45 to more divergent sequences, and BLOSUM62 is intermediate. The numbers do not specify a required identity of the query and target, nor a number of evolutionary time steps.^[ref]

Read the ten residue–residue columns from the BLOSUM62 substitution matrix. In their order in the supplied sequence alignment, their contributions areThere are three gap positions, each incurring the stated linear gap penalty of four, including the leading gap. Therefore the alignment score is $\boxed{24-3(4)=12}$.

An affine gap penalty charges a larger opening cost than extension cost, reflecting that a single insertion or deletion event can span several residues. For a gap of length $\ell \ge 1$, explicitly chooseThe supplied sequence alignment has two gap runs, of lengths one and two. Its residue–residue score remains $24$, so this choice givesMore generally, this affine gap penalty convention gives $S=24-2o-e$. There is no unique numerical rescore until opening and extension costs are chosen. A concave length-dependent gap penalty can make very long gaps progressively cheaper per added residue; distinct terminal gap penalties can accommodate incomplete sequence ends. Each choice encodes a different assumption about insertions, deletions or sequence completeness.

Use global sequence alignment, for example the Needleman-Wunsch algorithm with an appropriate substitution matrix and gap penalty. The two full-length globin proteins are related by orthology and are from closely related species, so the expected correspondence extends across both complete sequences. Unmatched ends should normally contribute to the score rather than being discarded.

With a linear gap penalty of magnitude $g$, the dynamic programming entry $F_{ij}$ is the best score for the first $i$ and $j$ residues. It satisfieswith $F_{i0}=-ig$, $F_{0j}=-jg$, and traceback from $F_{mn}$. An affine gap penalty uses separate match and gap states.

Use local sequence alignment, such as the Smith-Waterman algorithm, to compare the globin fragment with the matching region of each full-length target. A global sequence alignment would penalize long target flanks absent from the fragment. For a linear gap penalty the dynamic programming recurrence includes a restart option:Initialize boundary entries to zero, start traceback at the largest entry anywhere, and stop at zero. If the entire query is known to be valid and must be retained, fitting sequence alignment is a suitable refinement: align the whole fragment while leaving target flanks unpenalized.

First check the sequence alphabet: the prompt calls the targets “genes”. If the FASTA format file contains DNA rather than protein, translate the relevant coding sequences into amino acids before a protein-to-protein sequence alignment. If coding frames are unknown, search translated frames. A file in FASTA format specifies storage syntax, not whether the contents are nucleotides or residues.

The gain comes from avoiding a full dynamic programming table for every target. BLAST searches for short, sufficiently high-scoring word matches between query and database, then extends promising seeds. Modern versions reserve expensive gapped refinement for candidates passing earlier filters. FASTA sequence search similarly starts from short word matches, collects promising runs along alignment diagonals, and refines or joins selected regions. Word lookup and database indexing make candidate discovery inexpensive compared with considering every pair of positions.^[ref][ref]

These are heuristic algorithms: a good local sequence alignment can be missed if it has no acceptable seed, contains many interruptions, or is filtered before extension. Longer words and higher seed-score thresholds usually speed the search but reduce sensitivity; shorter words and looser thresholds recover weaker similarities at greater cost. Repetitive or compositionally biased regions can create many uninformative hits and erode both speed and specificity. No guarantee remains of finding the optimal sequence alignment for the chosen scoring scheme.

Thus the practical work depends on word length, thresholds, database composition and the number of hits. The speed-up is substantial on typical data, but does not establish a universal improved worst-case bound over $O(mn)$: an input producing many candidate hits can still require extensive refinement.