AlignIO

Core Principle: Bio.AlignIO provides standardized input/output interfaces for parsing, converting, column slicing, and analyzing Multiple Sequence Alignments (MSAs) across standard bioinformatic formats.

The interface

from Bio import AlignIO

aln = AlignIO.read("family.aln", "clustal")       # ONE alignment
for a in AlignIO.parse("multi.sto", "stockholm"): # many alignments
    ...
AlignIO.write(aln, "out.fasta", "fasta")
AlignIO.convert("in.aln", "clustal", "out.sto", "stockholm")

Same read/parse/write/convert shape as SeqIO: read for exactly one, parse for many.

Formats

Format string	Notes
"fasta"	aligned FASTA (gaps as -)
"clustal"	ClustalW/Omega .aln
"stockholm"	Pfam/Rfam; carries per-column annotation
"phylip", "phylip-relaxed"	strict form truncates names to 10 chars!
"nexus"	used by MrBayes, PAUP
"maf"	multiple genome alignment
"emboss", "msf", "mauve"	others

"phylip" truncates sequence names to 10 characters. Silently. Use "phylip-relaxed" unless a downstream tool demands strict PHYLIP.

Working with the alignment

len(aln)                          # number of SEQUENCES (rows)
aln.get_alignment_length()        # number of COLUMNS
aln[0]                            # first SeqRecord
aln[:, 10]                        # column 10 as a string   ← the key idiom
aln[:, 10:20]                     # a sub-alignment (columns)
aln[0:5, 10:20]                   # rows and columns
aln[0].seq                        # gapped sequence of the first record

for rec in aln:
    print(rec.id, str(rec.seq)[:50])

aln[:, i] returning a plain string of that column is the workhorse for conservation analysis.

Building and combining

from Bio.Align import MultipleSeqAlignment
aln = MultipleSeqAlignment([rec1, rec2, rec3])
aln.append(rec4)
aln.extend([rec5, rec6])

aln[:, :100] + aln[:, 200:]       # concatenate column ranges
aln.sort(key=lambda r: r.id)

Every record in a MultipleSeqAlignment must be the same length. Adding a record of a different length raises: which is the point.

Per-column analysis

import numpy as np
from collections import Counter

L = aln.get_alignment_length()

# gap fraction per column
gaps = np.array([aln[:, i].count("-") / len(aln) for i in range(L)])

# Shannon entropy per column (lower = more conserved)
def entropy(col):
    col = col.replace("-", "")
    if not col:
        return np.nan
    counts = np.array(list(Counter(col).values()), dtype=float)
    p = counts / counts.sum()
    return -(p * np.log2(p)).sum()

ent = np.array([entropy(aln[:, i]) for i in range(L)])

# consensus
consensus = "".join(Counter(aln[:, i]).most_common(1)[0][0] for i in range(L))

# columns present in ≥80% of sequences
keep = np.flatnonzero(gaps <= 0.2)
trimmed = aln[:, int(keep[0]):int(keep[-1]) + 1]     # contiguous range only

Note that slicing columns with a list of indices is not supported directly: you concatenate ranges, or convert to an array (below).

Alignment as a NumPy array

The fast route for column statistics:

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

gap_frac = (arr == "-").mean(axis=0)                     # per column
identity = (arr == arr[0]).mean(axis=0)                  # identity to the first seq
np.unique(arr[:, 42], return_counts=True)                # composition of column 42

# select non-gappy columns: arbitrary column sets work here
arr_trimmed = arr[:, gap_frac <= 0.2]

For a 1000-sequence × 5000-column alignment, the NumPy route is orders of magnitude faster than a Python loop over aln[:, i]. See Vectorization and Axes and Reductions.

Running an aligner

Bio.Align.Applications was removed in Biopython 1.86: no more MuscleCommandline, MafftCommandline, ClustalOmegaCommandline. Use subprocess:

import subprocess
from pathlib import Path

def run_mafft(in_fasta, out_fasta, threads=4):
    with open(out_fasta, "w") as fh:
        subprocess.run(
            ["mafft", "--auto", "--thread", str(threads), str(in_fasta)],
            stdout=fh, check=True,
        )
    return AlignIO.read(out_fasta, "fasta")

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

check=True raises on a non-zero exit code, which is what you want: a silently failed aligner produces an empty or truncated file that then parses "successfully".

Capture stderr for diagnostics:

p = subprocess.run([...], capture_output=True, text=True, check=True)
if p.stderr:
    print(p.stderr[-2000:])

Codon-aware alignment

Aligning coding sequences at the nucleotide level introduces frameshifting gaps. The correct procedure is to align the proteins and then back-translate the alignment onto the codons:

# 1. translate → 2. align proteins → 3. map gaps back onto the DNA
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)

Tools like pal2nal and MACSE do this properly, including handling frameshifts and stops. Codon-aware alignment matters for any dN/dS analysis: a nucleotide-level alignment will put gaps mid-codon and destroy the reading frame.

Bioinformatics examples

# conservation profile for a protein family
aln = AlignIO.read("kinase_domain.aln", "clustal")
arr = np.array([list(str(r.seq)) for r in aln])
gap = (arr == "-").mean(axis=0)
ent = np.array([entropy(aln[:, i]) for i in range(aln.get_alignment_length())])

fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(ent, lw=0.8, color="steelblue")
ax.fill_between(range(len(gap)), 0, gap * ent.max(), alpha=0.2,
                color="grey", label="gap fraction")
ax.set(xlabel="Alignment column", ylabel="Shannon entropy (bits)")
ax.legend()

# pairwise identity matrix from an MSA
n = len(aln)
pid = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        both = (arr[i] != "-") & (arr[j] != "-")
        pid[i, j] = (arr[i][both] == arr[j][both]).mean()

sns.clustermap(pd.DataFrame(pid, index=[r.id for r in aln],
                            columns=[r.id for r in aln]),
               cmap="rocket", vmin=0, vmax=1)

# is a variant position conserved?
col = aln[:, aln_position]
print(f"{col.count(ref_aa)}/{len(col)} sequences carry the reference residue")

# alignment → tree
AlignIO.convert("family.aln", "clustal", "family.phy", "phylip-relaxed")
subprocess.run(["iqtree2", "-s", "family.phy", "-m", "MFP"], check=True)

Common mistakes

• Bio.Align.Applications. Removed in 1.86; use subprocess.
• "phylip" truncating names to 10 characters. Use "phylip-relaxed".
• AlignIO.read on a multi-alignment file → ValueError.
• Nucleotide alignment of coding sequences. Align proteins, back-translate.
• Python loops over columns on a large alignment. Convert to a NumPy array.
• Forgetting gap characters in composition counts: - is not an amino acid.
• Comparing alignment columns to ungapped sequence positions. They are different coordinate systems; you must map through the gaps.
• subprocess.run without check=True, silently accepting a failed alignment.

See also

SeqIO · PairwiseAligner · Substitution Matrices · Bio.Phylo · Applied - Multiple Sequence Alignment Analysis · Biopython Deprecations

Exercise

Biopython is not installed in this sandbox, so rather than calling AlignIO.read, you will parse a Stockholm alignment yourself. The point is the same one the lesson makes about aln[:, i]: an alignment is a rectangle of characters, and every conservation statistic is a loop over its columns.

Write summarize(text), taking the whole Stockholm file as one string and returning a list of output lines. Ignore blank lines, any line beginning with #, and the terminating //; every other line is name followed by whitespace and a block of residues. A name may appear more than once, because Stockholm is interleaved: append each new block to that sequence in the order the names were first seen. Print sequences=<n> columns=<L>, then consensus=<string>, then conserved=<count>. The consensus takes the most frequent non-gap character in each column, breaking ties in favour of the character that sorts first, and uses - for a column that is entirely gaps. A column counts as conserved only when it contains no gap and every character is identical. If the file contains no sequences at all, print only the sequences=0 columns=0 line.