Applied - Multiple Sequence Alignment Analysis

In one line: running an aligner, converting the alignment to a NumPy array, and computing the conservation statistics that make it useful.

Uses: AlignIO · PairwiseAligner · Substitution Matrices · ndarray · Axes and Reductions · Bio.Phylo

Running an aligner

Bio.Align.Applications was removed in Biopython 1.86. Use subprocess:

import subprocess, shutil
from pathlib import Path
from Bio import AlignIO

def run_mafft(in_fasta, out_fasta, threads=4, mode="--auto"):
    if shutil.which("mafft") is None:
        raise FileNotFoundError("mafft not on PATH")
    with open(out_fasta, "w") as fh:
        subprocess.run(["mafft", mode, "--thread", str(threads), str(in_fasta)],
                       stdout=fh, stderr=subprocess.PIPE, text=True, check=True)
    return AlignIO.read(out_fasta, "fasta")

aln = run_mafft("orthologs.faa", "orthologs.aln")

check=True is essential — a silently failed aligner leaves an empty file that parses "successfully" into a zero-sequence alignment.

Aligner	Use
MAFFT	general purpose; `--auto` picks a strategy by input size
MUSCLE	good accuracy on small sets
Clustal Omega	scales to very large sets
T-Coffee	slow, high accuracy
HMMER `hmmalign`	align to an existing profile — the right tool for domain families
MACSE / pal2nal	codon-aware for coding sequences

Alignment → NumPy array

The move that makes everything fast:

import numpy as np

arr = np.array([list(str(rec.seq)) for rec in aln])    # (n_seqs, n_cols) of chars
ids = [rec.id for rec in aln]
n_seq, n_col = arr.shape

Now every column statistic is a vectorised reduction:

gap_frac = (arr == "-").mean(axis=0)                    # per column
identity_to_ref = (arr == arr[0]).mean(axis=0)
n_variants = np.array([len(set(c) - {"-"}) for c in arr.T])

axis=0 collapses sequences, giving one value per column. Getting this backwards gives per-sequence statistics — also useful, entirely different. See Axes and Reductions.

Conservation measures

Shannon entropy

from collections import Counter

def column_entropy(col, ignore_gaps=True):
    c = [x for x in col if x != "-"] if ignore_gaps else list(col)
    if not c:
        return np.nan
    p = np.array(list(Counter(c).values()), dtype=float)
    p /= p.sum()
    return float(-(p * np.log2(p)).sum())

entropy = np.array([column_entropy(arr[:, i]) for i in range(n_col)])

Low entropy = conserved. Maximum is log₂(20) ≈ 4.32 bits for protein, 2 bits for DNA.

ignore_gaps is a real decision. A column that is 90% gaps has low entropy over its three residues and is not "conserved" in any meaningful sense — it is barely observed. Report the gap fraction alongside entropy, always.

Sum-of-pairs with a substitution matrix

from Bio.Align import substitution_matrices
from itertools import combinations

M = substitution_matrices.load("BLOSUM62")

def column_sop(col):
    res = [x for x in col if x != "-" and x in M.alphabet]
    if len(res) < 2:
        return np.nan
    return float(np.mean([M[a, b] for a, b in combinations(res, 2)]))

sop = np.array([column_sop(arr[:, i]) for i in range(n_col)])

Better than entropy for proteins, because it knows that leucine→isoleucine is conservative and leucine→proline is not. Entropy treats every substitution as equally different. See Substitution Matrices.

Jensen-Shannon divergence from background

AA = "ACDEFGHIKLMNPQRSTVWY"
BG = np.array([0.074,0.025,0.054,0.054,0.047,0.074,0.026,0.068,0.058,0.099,
               0.025,0.045,0.039,0.034,0.052,0.057,0.051,0.073,0.013,0.032])

def jsd(col, bg=BG, pseudocount=1e-6):
    res = [x for x in col if x in AA]
    if len(res) < 2:
        return np.nan
    p = np.array([res.count(a) for a in AA], float) + pseudocount
    p /= p.sum()
    m = (p + bg) / 2
    kl = lambda a, b: np.sum(a * np.log2(a / b))
    return float(0.5 * kl(p, m) + 0.5 * kl(bg, m))

conservation = np.array([jsd(arr[:, i]) for i in range(n_col)])

JSD against a background amino-acid distribution is generally the best-performing simple conservation score. It rewards columns that are unusual relative to background, not merely uniform — a column of all-alanine (a common residue) is less informative than a column of all-tryptophan (rare).

Plotting the profile

import matplotlib.pyplot as plt, seaborn as sns
sns.set_theme(style="ticks", context="paper")

fig, axes = plt.subplots(2, 1, figsize=(11, 4.5), sharex=True,
                         gridspec_kw=dict(height_ratios=[2, 1]),
                         layout="constrained")

axes[0].plot(conservation, lw=0.7, c="navy")
axes[0].set_ylabel("JSD conservation")

axes[1].fill_between(range(n_col), gap_frac, step="mid", alpha=0.5, color="grey")
axes[1].set(ylabel="Gap fraction", xlabel="Alignment column", ylim=(0, 1))

for s, e in domain_bounds:                 # highlight annotated domains
    for ax in axes:
        ax.axvspan(s, e, alpha=0.12, color="steelblue", zorder=0)

Plotting gap fraction beneath conservation is what keeps the figure honest — high conservation in a mostly-gapped column is an artefact, and the reader can see it.

Trimming

keep_cols = gap_frac <= 0.2
arr_trim = arr[:, keep_cols]
print(f"kept {keep_cols.sum()} of {n_col} columns")

# poorly-aligned sequences
seq_gap = (arr == "-").mean(axis=1)
keep_seqs = seq_gap < 0.5
arr_trim = arr_trim[keep_seqs]

For real trimming use trimAl or Gblocks — they use column-neighbourhood criteria rather than a per-column threshold, which matters for phylogenetics.

Trimming changes your tree. Over-trimming removes signal; under-trimming adds noise. Report what you did.

Pairwise identity matrix

n = len(arr)
pid = np.zeros((n, n))
for i in range(n):
    for j in range(i, n):
        both = (arr[i] != "-") & (arr[j] != "-")
        v = (arr[i][both] == arr[j][both]).mean() if both.any() else np.nan
        pid[i, j] = pid[j, i] = v

import pandas as pd
sns.clustermap(pd.DataFrame(pid, index=ids, columns=ids),
               cmap="rocket", vmin=0, vmax=1, figsize=(8, 8),
               cbar_kws=dict(label="Pairwise identity"))

Note the denominator: identity over positions where both sequences have a residue. The alternative — identity over all alignment columns, counting gaps as mismatches — gives a substantially lower number. Both are used; say which.

Mapping to ungapped coordinates

Alignment columns are not sequence positions. To ask "is residue 175 of my protein conserved?":

def col_to_pos(seq_row):
    """alignment column index → 1-based ungapped position (or None for gaps)."""
    pos, out = 0, []
    for c in seq_row:
        if c == "-":
            out.append(None)
        else:
            pos += 1
            out.append(pos)
    return out

def pos_to_col(seq_row):
    return {p: i for i, p in enumerate(col_to_pos(seq_row)) if p is not None}

ref_map = pos_to_col(arr[0])
col = ref_map[175]
print(f"residue 175 → column {col}, conservation {conservation[col]:.2f}, "
      f"gap fraction {gap_frac[col]:.2f}")
print("observed residues:", Counter(arr[:, col]))

This mapping is where variant-interpretation code goes wrong. Reporting "column 412 is highly conserved" when the user asked about residue 175 is a silent, confident error.

Codon-aware alignment

Aligning coding DNA directly introduces frameshifting gaps. Align proteins, then back-translate:

def back_translate(prot_aln, dna_records):
    dna = {r.id: str(r.seq) for r in dna_records}
    out = []
    for rec in prot_aln:
        nt, i, parts = dna[rec.id], 0, []
        for aa in str(rec.seq):
            if aa == "-":
                parts.append("---")
            else:
                parts.append(nt[i:i+3]); i += 3
        out.append(SeqRecord(Seq("".join(parts)), id=rec.id, description=""))
    return MultipleSeqAlignment(out)

Essential for any dN/dS analysis — a nucleotide-level alignment puts gaps mid-codon and destroys the reading frame. pal2nal and MACSE do this properly, including frameshifts and internal stops.

Feeding a tree

AlignIO.convert("family.aln", "clustal", "family.phy", "phylip-relaxed")
subprocess.run(["iqtree2", "-s", "family.phy", "-m", "MFP",
                "-B", "1000", "-T", "AUTO"], check=True)

from Bio import Phylo
tree = Phylo.read("family.phy.treefile", "newick")
tree.ladderize()

"phylip-relaxed", not "phylip" — strict PHYLIP truncates names to 10 characters, silently. See AlignIO and Bio.Phylo.

Common mistakes

Bio.Align.Applications. Removed in 1.86.
Nucleotide alignment of coding sequences.
Column index confused with residue position.
Conservation without reporting gap fraction.
Ambiguous percent identity denominator.
Python loops over columns on a large alignment. Use the array.
"phylip" truncating names.
Not stating the trimming procedure.
subprocess.run without check=True.
Treating a conservation score as functional importance — conservation reflects constraint, which is necessary but not sufficient.

Applied - Multiple Sequence Alignment Analysis

Applied - Multiple Sequence Alignment Analysis

Running an aligner

Alignment → NumPy array

Conservation measures

Shannon entropy

Sum-of-pairs with a substitution matrix

Jensen-Shannon divergence from background

Plotting the profile

Trimming

Pairwise identity matrix

Mapping to ungapped coordinates

Codon-aware alignment

Feeding a tree

Common mistakes

See also