Substitution Matrices

In one line: log-odds scores for aligning one residue against another — BLOSUM62 by default, and the choice matters more than people assume.

Loading

from Bio.Align import substitution_matrices

substitution_matrices.load()                    # list all bundled matrices
m = substitution_matrices.load("BLOSUM62")

m["A", "A"]        # 4.0
m["W", "C"]        # -2.0
m.alphabet         # 'ARNDCQEGHILKMFPSTWYVBZX*'
m.shape
print(m)

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

from Bio import Align

aligner = Align.PairwiseAligner()
aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
aligner.open_gap_score = -11
aligner.extend_gap_score = -1

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.

m = substitution_matrices.load("BLOSUM62")
m["L", "I"]     # +2  — both aliphatic hydrophobic, freely interchangeable
m["L", "P"]     # -3  — proline breaks helices; rarely tolerated
m["W", "W"]     # +11 — tryptophan is rare, so a match is highly informative
m["A", "A"]     # +4  — alanine is common, so a match says less

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:

aligner = Align.PairwiseAligner(match_score=2, mismatch_score=-3,
                                open_gap_score=-5, extend_gap_score=-2)

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.

See also

PairwiseAligner · AlignIO · Bio.Blast · Biopython Deprecations · Applied - Multiple Sequence Alignment Analysis