Beginner·4 graded tests

Substitution Matrices

Core Principle: Substitution matrices parameterize sequence alignment algorithms with empirically derived log-odds substitution scores, where matrix selection (e.g. BLOSUM62 versus PAM250) reflects evolutionary divergence.

Loading

[1]python

out1ms

#  Matrix made by matblas from blosum62.iij
#  * column uses minimum score
#  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
#  Blocks Database = /data/blocks_5.0/blocks.dat
#  Cluster Percentage: >= 62
#  Entropy =   0.6979, Expected =  -0.5209
     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.0 -1.0 -2.0 -2.0  0.0 -1.0 -1.0  0.0 -2.0 -1.0 -1.0 -1.0 -1.0 -2.0 -1.0  1.0  0.0 -3.0 -2.0  0.0 -2.0 -1.0  0.0 -4.0
R -1.0  5.0  0.0 -2.0 -3.0  1.0  0.0 -2.0  0.0 -3.0 -2.0  2.0 -1.0 -3.0 -2.0 -1.0 -1.0 -3.0 -2.0 -3.0 -1.0  0.0 -1.0 -4.0
N -2.0  0.0  6.0  1.0 -3.0  0.0  0.0  0.0  1.0 -3.0 -3.0  0.0 -2.0 -3.0 -2.0  1.0  0.0 -4.0 -2.0 -3.0  3.0  0.0 -1.0 -4.0
D -2.0 -2.0  1.0  6.0 -3.0  0.0  2.0 -1.0 -1.0 -3.0 -4.0 -1.0 -3.0 -3.0 -1.0  0.0 -1.0 -4.0 -3.0 -3.0  4.0  1.0 -1.0 -4.0
C  0.0 -3.0 -3.0 -3.0  9.0 -3.0 -4.0 -3.0 -3.0 -1.0 -1.0 -3.0 -1.0 -2.0 -3.0 -1.0 -1.0 -2.0 -2.0 -1.0 -3.0 -3.0 -2.0 -4.0
Q -1.0  1.0  0.0  0.0 -3.0  5.0  2.0 -2.0  0.0 -3.0 -2.0  1.0  0.0 -3.0 -1.0  0.0 -1.0 -2.0 -1.0 -2.0  0.0  3.0 -1.0 -4.0
E -1.0  0.0  0.0  2.0 -4.0  2.0  5.0 -2.0  0.0 -3.0 -3.0  1.0 -2.0 -3.0 -1.0  0.0 -1.0 -3.0 -2.0 -2.0  1.0  4.0 -1.0 -4.0
G  0.0 -2.0  0.0 -1.0 -3.0 -2.0 -2.0  6.0 -2.0 -4.0 -4.0 -2.0 -3.0 -3.0 -2.0  0.0 -2.0 -2.0 -3.0 -3.0 -1.0 -2.0 -1.0 -4.0
H -2.0  0.0  1.0 -1.0 -3.0  0.0  0.0 -2.0  8.0 -3.0 -3.0 -1.0 -2.0 -1.0 -2.0 -1.0 -2.0 -2.0  2.0 -3.0  0.0  0.0 -1.0 -4.0
I -1.0 -3.0 -3.0 -3.0 -1.0 -3.0 -3.0 -4.0 -3.0  4.0  2.0 -3.0  1.0  0.0 -3.0 -2.0 -1.0 -3.0 -1.0  3.0 -3.0 -3.0 -1.0 -4.0
L -1.0 -2.0 -3.0 -4.0 -1.0 -2.0 -3.0 -4.0 -3.0  2.0  4.0 -2.0  2.0  0.0 -3.0 -2.0 -1.0 -2.0 -1.0  1.0 -4.0 -3.0 -1.0 -4.0
K -1.0  2.0  0.0 -1.0 -3.0  1.0  1.0 -2.0 -1.0 -3.0 -2.0  5.0 -1.0 -3.0 -1.0  0.0 -1.0 -3.0 -2.0 -2.0  0.0  1.0 -1.0 -4.0
M -1.0 -1.0 -2.0 -3.0 -1.0  0.0 -2.0 -3.0 -2.0  1.0  2.0 -1.0  5.0  0.0 -2.0 -1.0 -1.0 -1.0 -1.0  1.0 -3.0 -1.0 -1.0 -4.0
F -2.0 -3.0 -3.0 -3.0 -2.0 -3.0 -3.0 -3.0 -1.0  0.0  0.0 -3.0  0.0  6.0 -4.0 -2.0 -2.0  1.0  3.0 -1.0 -3.0 -3.0 -1.0 -4.0
P -1.0 -2.0 -2.0 -1.0 -3.0 -1.0 -1.0 -2.0 -2.0 -3.0 -3.0 -1.0 -2.0 -4.0  7.0 -1.0 -1.0 -4.0 -3.0 -2.0 -2.0 -1.0 -2.0 -4.0
S  1.0 -1.0  1.0  0.0 -1.0  0.0  0.0  0.0 -1.0 -2.0 -2.0  0.0 -1.0 -2.0 -1.0  4.0  1.0 -3.0 -2.0 -2.0  0.0  0.0  0.0 -4.0
T  0.0 -1.0  0.0 -1.0 -1.0 -1.0 -1.0 -2.0 -2.0 -1.0 -1.0 -1.0 -1.0 -2.0 -1.0  1.0  5.0 -2.0 -2.0  0.0 -1.0 -1.0  0.0 -4.0
W -3.0 -3.0 -4.0 -4.0 -2.0 -2.0 -3.0 -2.0 -2.0 -3.0 -2.0 -3.0 -1.0  1.0 -4.0 -3.0 -2.0 11.0  2.0 -3.0 -4.0 -3.0 -2.0 -4.0
Y -2.0 -2.0 -2.0 -3.0 -2.0 -1.0 -2.0 -3.0  2.0 -1.0 -1.0 -2.0 -1.0  3.0 -3.0 -2.0 -2.0  2.0  7.0 -1.0 -3.0 -2.0 -1.0 -4.0
V  0.0 -3.0 -3.0 -3.0 -1.0 -2.0 -2.0 -3.0 -3.0  3.0  1.0 -2.0  1.0 -1.0 -2.0 -2.0  0.0 -3.0 -1.0  4.0 -3.0 -2.0 -1.0 -4.0
B -2.0 -1.0  3.0  4.0 -3.0  0.0  1.0 -1.0  0.0 -3.0 -4.0  0.0 -3.0 -3.0 -2.0  0.0 -1.0 -4.0 -3.0 -3.0  4.0  1.0 -1.0 -4.0
Z -1.0  0.0  0.0  1.0 -3.0  3.0  4.0 -2.0  0.0 -3.0 -3.0  1.0 -1.0 -3.0 -1.0  0.0 -1.0 -3.0 -2.0 -2.0  1.0  4.0 -1.0 -4.0
X  0.0 -1.0 -1.0 -1.0 -2.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -2.0  0.0  0.0 -2.0 -1.0 -1.0 -1.0 -1.0 -1.0 -4.0
* -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0 -4.0  1.0

The object is a subclass of numpy.ndarray with string-index lookup, so NumPy operations work on it directly.

The bundled families

Family	Range	Best for
BLOSUM	BLOSUM45, 50, 62, 80, 90	protein; from observed blocks
PAM	PAM30, 70, 250	protein; from an evolutionary model
NUC	NUC.4.4	nucleotide
Others	GONNET, RAO, RISLER, MDM78, JOHNSON, SCHNEIDER	specialised

The BLOSUM number is the clustering identity threshold, and the relationship is counter-intuitive:

Higher number (BLOSUM80, 90) → for closely related sequences. Derived from blocks clustered at ≥80% identity.
Lower number (BLOSUM45, 50) → for distant homologues. More tolerant of substitution.

PAM runs the opposite way: higher PAM = more divergent (PAM250 for distant, PAM30 for close). Mixing up the two directions is the classic error.

BLOSUM62 is the default for BLASTP and a good general choice. Reach for BLOSUM45 when you are hunting remote homology and BLOSUM62 finds nothing.

Using with PairwiseAligner

[ ]python

out1ms

This code completed without producing output.

Gap penalties must match the matrix. They are on the same log-odds scale, so a gap penalty tuned for BLOSUM62 is wrong for PAM250. The BLAST defaults are a safe pairing:

Matrix	Gap open	Gap extend
BLOSUM62	−11	−1
BLOSUM80	−10	−1
BLOSUM45	−15	−2
PAM250	−14	−2

Using match/mismatch scores instead of a matrix for protein alignment throws away most of the biological signal: it treats a leucine→isoleucine substitution (conservative, score +2) the same as leucine→proline (disruptive, score −3).

Reading the values

The entries are log-odds:

score(a,b) = (1/λ) · log( P(a,b observed in alignments) / (P(a)·P(b)) )

Positive: the pair occurs more often in real alignments than chance predicts; a conservative substitution.
Zero: as expected by chance.
Negative: disfavoured.

[3]python

out1ms

4.0

That last contrast is worth sitting with: the diagonal is not constant, because matching a rare residue is stronger evidence of homology than matching a common one. This is exactly what a log-odds score encodes, and it is why substitution matrices beat percent identity as a similarity measure.

Building your own

from Bio.Align import substitution_matrices

# empirical substitution counts from an alignment
aln = aligner.align(seq1, seq2)[0]
counts = aln.substitutions
freq = counts / counts.sum()

You can also read a matrix from a file in the standard NCBI format, or construct one from an Array of counts. Useful for organism-specific or domain-specific scoring, though in practice the standard matrices are hard to beat without a lot of data.

Nucleotide scoring

Nucleotide alignment usually uses simple match/mismatch rather than a matrix, because there are only four letters and no chemistry to encode:

[ ]python

out1ms

This code completed without producing output.

Those are the BLASTN defaults. A transition/transversion-aware matrix (NUC.4.4, or a custom one) is worth it when you are working at the edge of detectability, since transitions (A↔G, C↔T) are ~2× more frequent than transversions.

Bio.SubsMat is gone

Removed in Biopython 1.80.

from Bio.SubsMat import MatrixInfo            # ImportError
MatrixInfo.blosum62                            # gone

# modern
from Bio.Align import substitution_matrices
substitution_matrices.load("BLOSUM62")

Old code and tutorials use MatrixInfo.blosum62, a plain dict keyed by tuples with only one triangle populated (so ("A","W") existed but ("W","A") did not: a genuinely annoying API). The replacement is symmetric and array-backed. See Biopython Deprecations.

Bioinformatics examples

# conservation score for each column of an MSA
from Bio import AlignIO
aln = AlignIO.read("family.aln", "clustal")
m = substitution_matrices.load("BLOSUM62")

def column_score(col):
    residues = [c for c in col if c != "-"]
    pairs = [(a, b) for i, a in enumerate(residues) for b in residues[i+1:]]
    return np.mean([m[a, b] for a, b in pairs]) if pairs else np.nan

scores = [column_score(aln[:, i]) for i in range(aln.get_alignment_length())]
plt.plot(scores)

# is a missense variant conservative or radical?
def severity(ref_aa, alt_aa):
    s = m[ref_aa, alt_aa]
    return "conservative" if s >= 0 else "radical"

variants["aa_score"] = [m[r, a] for r, a in zip(variants["ref_aa"], variants["alt_aa"])]
sns.boxplot(data=variants, x="clinical_significance", y="aa_score")

# a heatmap of the matrix itself
sns.heatmap(pd.DataFrame(np.array(m), index=list(m.alphabet),
                         columns=list(m.alphabet)),
            cmap="vlag", center=0, square=True, annot=False)

The BLOSUM score of a missense change is a cheap, interpretable feature for variant prioritisation: not as good as a dedicated predictor, but free and easy to explain.

Common mistakes

Bio.SubsMat. Removed in 1.80.
BLOSUM direction reversed. Higher = closer relatives.
Confusing BLOSUM and PAM directions. They run opposite ways.
Gap penalties not matched to the matrix.
Match/mismatch scores for protein alignment.
Comparing raw scores across different matrices. They are on different scales; use bit scores or E-values.
Forgetting X, B, Z and * exist in the alphabet: a lookup on an unexpected character raises KeyError.

Exercise

Biopython is not installed here, so read the matrix out of its NCBI-format text yourself and score an alignment with it. The format is comment lines beginning with #, then a header row of the alphabet, then one row per residue whose first field is that residue's letter. The matrix is symmetric, so store both orientations and the half-populated dict problem that Bio.SubsMat had never arises.

Write score_alignment(matrix, alphabet, top, bottom, gap_open, gap_extend) returning (score, identities, mismatches, gaps). Score each column in which neither sequence has a gap by its matrix entry, replacing any residue outside the alphabet with X before the lookup. Score a run of k consecutive gaps in one sequence as -(gap_open + gap_extend * (k - 1)), so an affine penalty is charged once for opening and once per extension; charging gap_open per position instead is the mistake this is testing for. A column gapped in both sequences contributes nothing and does not interrupt a run. Print score=<x> rounded to one decimal, then identities=<n>, mismatches=<n> and gaps=<n> counting gapped columns. The input is the matrix block terminated by a line containing only //, then a line holding gap_open gap_extend as positive numbers, then the two aligned sequences.

exerciseEvery submission runs independently in a fresh namespace.

Test cases · 4

#	via	input	expected stdout
1	stdin	# a reduced BLOSUM62 A C D G I L P S T V W X A 4 0 -2 0 -1 -1 -1 1 0 0 -3 0 C 0 9 -3 -3 -1 -1 -3 -1 -1 -1 -2 -2 D -2 -3 6 -1 -3 -4 -1 0 -1 -3 -4 -1 G 0 -3 -1 6 -4 -4 -2 0 -2 -3 -2 -1 I -1 -1 -3 -4 4 2 -3 -2 -1 3 -3 -1 L -1 -1 -4 -4 2 4 -3 -2 -1 1 -2 -1 P -1 -3 -1 -2 -3 -3 7 -1 -1 -2 -4 -2 S 1 -1 0 0 -2 -2 -1 4 1 -2 -3 0 T 0 -1 -1 -2 -1 -1 -1 1 5 0 -2 0 V 0 -1 -3 -3 3 1 -2 -2 0 4 -3 -1 W -3 -2 -4 -2 -3 -2 -4 -3 -2 -3 11 -2 X 0 -2 -1 -1 -1 -1 -2 0 0 -1 -2 -1 // 11 1 WTAPILSGCD WTAPVLSGCD	score=59.0 identities=9 mismatches=1 gaps=0
2	stdin	# a reduced BLOSUM62 A C D G I L P S T V W X A 4 0 -2 0 -1 -1 -1 1 0 0 -3 0 C 0 9 -3 -3 -1 -1 -3 -1 -1 -1 -2 -2 D -2 -3 6 -1 -3 -4 -1 0 -1 -3 -4 -1 G 0 -3 -1 6 -4 -4 -2 0 -2 -3 -2 -1 I -1 -1 -3 -4 4 2 -3 -2 -1 3 -3 -1 L -1 -1 -4 -4 2 4 -3 -2 -1 1 -2 -1 P -1 -3 -1 -2 -3 -3 7 -1 -1 -2 -4 -2 S 1 -1 0 0 -2 -2 -1 4 1 -2 -3 0 T 0 -1 -1 -2 -1 -1 -1 1 5 0 -2 0 V 0 -1 -3 -3 3 1 -2 -2 0 4 -3 -1 W -3 -2 -4 -2 -3 -2 -4 -3 -2 -3 11 -2 X 0 -2 -1 -1 -1 -1 -2 0 0 -1 -2 -1 // 11 1 WTAPILSGCD WT---LSGCD	score=32.0 identities=7 mismatches=0 gaps=3
3	stdin	# a reduced BLOSUM62 A C D G I L P S T V W X A 4 0 -2 0 -1 -1 -1 1 0 0 -3 0 C 0 9 -3 -3 -1 -1 -3 -1 -1 -1 -2 -2 D -2 -3 6 -1 -3 -4 -1 0 -1 -3 -4 -1 G 0 -3 -1 6 -4 -4 -2 0 -2 -3 -2 -1 I -1 -1 -3 -4 4 2 -3 -2 -1 3 -3 -1 L -1 -1 -4 -4 2 4 -3 -2 -1 1 -2 -1 P -1 -3 -1 -2 -3 -3 7 -1 -1 -2 -4 -2 S 1 -1 0 0 -2 -2 -1 4 1 -2 -3 0 T 0 -1 -1 -2 -1 -1 -1 1 5 0 -2 0 V 0 -1 -3 -3 3 1 -2 -2 0 4 -3 -1 W -3 -2 -4 -2 -3 -2 -4 -3 -2 -3 11 -2 X 0 -2 -1 -1 -1 -1 -2 0 0 -1 -2 -1 // 11 1 WTAP-L-SGCUD WT-PILAS--CD	score=-10.0 identities=6 mismatches=1 gaps=5
4	stdin	# a reduced BLOSUM62 A C D G I L P S T V W X A 4 0 -2 0 -1 -1 -1 1 0 0 -3 0 C 0 9 -3 -3 -1 -1 -3 -1 -1 -1 -2 -2 D -2 -3 6 -1 -3 -4 -1 0 -1 -3 -4 -1 G 0 -3 -1 6 -4 -4 -2 0 -2 -3 -2 -1 I -1 -1 -3 -4 4 2 -3 -2 -1 3 -3 -1 L -1 -1 -4 -4 2 4 -3 -2 -1 1 -2 -1 P -1 -3 -1 -2 -3 -3 7 -1 -1 -2 -4 -2 S 1 -1 0 0 -2 -2 -1 4 1 -2 -3 0 T 0 -1 -1 -2 -1 -1 -1 1 5 0 -2 0 V 0 -1 -3 -3 3 1 -2 -2 0 4 -3 -1 W -3 -2 -4 -2 -3 -2 -4 -3 -2 -3 11 -2 X 0 -2 -1 -1 -1 -1 -2 0 0 -1 -2 -1 // 11 1 WTAPILSGCD WTAPILSGCD	score=60.0 identities=10 mismatches=0 gaps=0

import sys

def load_matrix(block):
    """Parse an NCBI-format matrix block into (matrix_dict, alphabet_set)."""
    header = None
    matrix = {}
    for line in block.splitlines():
        if not line.strip() or line.lstrip().startswith("#"):
            continue
        fields = line.split()
        if header is None:
            header = fields
            continue
        for column, value in zip(header, fields[1:]):
            matrix[(fields[0], column)] = float(value)
            matrix[(column, fields[0])] = float(value)
    return matrix, set(header or [])

def score_alignment(matrix, alphabet, top, bottom, gap_open, gap_extend):
    """Return (score, identities, mismatches, gaps) for a gapped pair.

Gap runs are charged gap_open once plus gap_extend per extra position.
    """
    # TODO: walk the columns, scoring pairs and collapsing gap runs
    pass

if __name__ == "__main__":
    block, _, rest = sys.stdin.read().partition("\n//\n")
    rows = [r for r in rest.split("\n") if r.strip()]
    matrix, alphabet = load_matrix(block)
    opening, extending = (float(v) for v in rows[0].split())
    score, ident, mism, gaps = score_alignment(
        matrix, alphabet, rows[1].strip(), rows[2].strip(), opening, extending)
    print("score=%.1f" % score)
    print("identities=%d" % ident)
    print("mismatches=%d" % mism)
    print("gaps=%d" % gaps)

Hints · 3

01Hint

The first non-comment line is the alphabet; every later row starts with its own letter.

02Hint

Store matrix[(a, b)] and matrix[(b, a)] as you read each value.

03Hint

Walk the columns with a while loop so you can advance past a whole gap run at once.

Substitution Matrices

Loading

The bundled families

Using with PairwiseAligner

Reading the values

Building your own

Nucleotide scoring

Bio.SubsMat is gone

Bioinformatics examples

Common mistakes

See also

Exercise

Test cases · 4

Hints · 3