Applied - Multiple Sequence Alignment Analysis

Core Principle: Multiple Sequence Alignment (MSA) analysis evaluates evolutionary sequence constraint, positional conservation scores, and pairwise similarity metrics across homologous sequences.

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.

See also

AlignIO · PairwiseAligner · Substitution Matrices · Bio.Phylo · Axes and Reductions · Applied - Sequence Composition Analysis · Applied - Annotating a Variant Table

Exercise

Biopython is not installed here, so parse the aligned FASTA by hand rather than through AlignIO, and compute the conservation profile directly. Everything the lesson does with np.frombuffer on the alignment array you can do with indexing into the strings.

Write analyse(text), returning a list of output lines. The first input line holds the maximum tolerated gap fraction for a column. What follows is an aligned FASTA: a > header line begins each record, and its sequence may be wrapped over several lines that must be concatenated. Uppercase every residue and treat - as the gap character. If the records are not all the same length the file is not an alignment: print error: ragged alignment and nothing else.

Compute for each column its gap fraction and its Shannon entropy in bits over the non-gap characters only, treating an all-gap column as entropy zero. Keep the columns whose gap fraction is at or below the threshold. Print sequences=<n> columns=<L> kept=<n>, then mean_entropy=<v>, the mean entropy over the kept columns rounded to three decimals and 0.0 when no column is kept, then invariant=<n>, the number of kept columns that carry no gap and a single distinct residue. Then print one line per record in input order holding the record name -- the first whitespace-separated token of the header, without the > -- and the number of non-gap residues it retains in the trimmed alignment, tab separated. An input with no records must print sequences=0 columns=0 kept=0 followed by no alignment.