Bio Alignment Io

Version Compatibility

Reference examples tested with: BioPython 1.83+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Alignment File I/O

Read, write, and convert multiple sequence alignment files in various formats.

Required Import

Goal: Load modules for reading, writing, and manipulating multiple sequence alignments.

Approach: Import AlignIO for file I/O and supporting classes for programmatic alignment construction.

from Bio import AlignIO
from Bio.Align import MultipleSeqAlignment
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq

Format Coverage Map

Three Python libraries cover the alignment-format space, with overlapping but non-identical support. Pick by what is actually required.

Formats NOT in BioPython (use dedicated tools):

Recommend Bio.Align (modern API) over Bio.AlignIO (legacy) for new code; it returns Alignment objects with built-in .counts() and .substitutions properties. For multi-gigabyte Stockholm databases such as Pfam-A.full, pyhmmer.easel.MSAFile streams record-by-record where Bio.AlignIO.parse works but at higher per-record cost.

Reading Alignments

"Read an alignment file" -> Parse an alignment file into an alignment object with sequences and metadata accessible.

Goal: Load alignment data from files in various formats (Clustal, PHYLIP, Stockholm, FASTA).

Approach: Use AlignIO.read() for single-alignment files or AlignIO.parse() for files containing multiple alignments.

Single Alignment File

from Bio import AlignIO

alignment = AlignIO.read('alignment.aln', 'clustal')
print(f'Alignment length: {alignment.get_alignment_length()}')
print(f'Number of sequences: {len(alignment)}')

Multiple Alignments in One File

for alignment in AlignIO.parse('multi_alignment.sto', 'stockholm'):
    print(f'Alignment with {len(alignment)} sequences, length {alignment.get_alignment_length()}')

Read as List

alignments = list(AlignIO.parse('alignments.phy', 'phylip'))
print(f'Read {len(alignments)} alignments')

Writing Alignments

Goal: Save alignment data to files in standard formats for downstream tools or archival.

Approach: Use AlignIO.write() with the target format specifier, supporting single or multiple alignments and file handles.

Write Single Alignment

AlignIO.write(alignment, 'output.fasta', 'fasta')

Write Multiple Alignments

alignments = [alignment1, alignment2, alignment3]
count = AlignIO.write(alignments, 'output.sto', 'stockholm')
print(f'Wrote {count} alignments')

Write to Handle

with open('output.aln', 'w') as handle:
    AlignIO.write(alignment, handle, 'clustal')

Format Conversion

"Convert alignment format" -> Transform an alignment file from one format to another (e.g., Clustal to PHYLIP).

Goal: Convert alignment files between formats for compatibility with different analysis tools.

Approach: Use AlignIO.convert() for direct one-step conversion, or read-modify-write for cases requiring intermediate manipulation.

Direct Conversion (Most Efficient)

AlignIO.convert('input.aln', 'clustal', 'output.phy', 'phylip')

With Alphabet Specification

AlignIO.convert('input.sto', 'stockholm', 'output.nex', 'nexus', molecule_type='DNA')

Manual Conversion (When Modification Needed)

alignment = AlignIO.read('input.aln', 'clustal')
# ... modify alignment ...
AlignIO.write(alignment, 'output.fasta', 'fasta')

Accessing Alignment Data

Goal: Navigate and extract data from alignment objects including sequences, columns, and slices.

Approach: Use iteration, indexing, and column slicing on the alignment object.

alignment = AlignIO.read('alignment.aln', 'clustal')

# Iterate over sequences
for record in alignment:
    print(f'{record.id}: {record.seq}')

# Access by index
first_seq = alignment[0]
last_seq = alignment[-1]

# Slice columns
column_slice = alignment[:, 10:20]  # Columns 10-19

# Get specific column
column = alignment[:, 5]  # Column 5 as string

Working with Alignment Objects

Get Alignment Properties

alignment = AlignIO.read('alignment.aln', 'clustal')

length = alignment.get_alignment_length()
num_seqs = len(alignment)
seq_ids = [record.id for record in alignment]

Slice Alignments

# Get subset of sequences
subset = alignment[0:5]  # First 5 sequences

# Get subset of columns
trimmed = alignment[:, 50:150]  # Columns 50-149

# Combine slicing
region = alignment[0:5, 50:150]  # 5 sequences, columns 50-149

Creating Alignments Programmatically

Goal: Build an alignment object from sequences defined in code rather than read from a file.

Approach: Construct SeqRecord objects with gap characters and wrap them in a MultipleSeqAlignment.

from Bio.Align import MultipleSeqAlignment
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq

records = [
    SeqRecord(Seq('ACTGACTGACTG'), id='seq1'),
    SeqRecord(Seq('ACTGACT-ACTG'), id='seq2'),
    SeqRecord(Seq('ACTG-CTGACTG'), id='seq3'),
]
alignment = MultipleSeqAlignment(records)
AlignIO.write(alignment, 'new_alignment.fasta', 'fasta')

Format Selection for Downstream Tools

Choosing the output format depends on which downstream tool consumes the alignment:

Annotation Preservation

Not all formats support annotations. Converting between formats can silently discard metadata:

Converting Stockholm to FASTA or PHYLIP discards all annotations, secondary structure markup, and per-residue quality scores. If annotations matter, keep a Stockholm master copy.

Format-Specific Notes

PHYLIP Format Pitfalls

PHYLIP has two incompatible variants (interleaved vs sequential) and two name-length modes (strict vs relaxed). Confusing these causes silent data corruption.

Strict PHYLIP truncates sequence names to exactly 10 characters. This can silently merge distinct sequences whose names share a 10-character prefix (e.g., Homo_sapiens_chr1 and Homo_sapiens_chr2 both become Homo_sapie).

# Strict PHYLIP (10-char names, interleaved) -- only for tools requiring it
alignment = AlignIO.read('file.phy', 'phylip')

# Sequential PHYLIP (10-char names, one sequence at a time) -- PAML/codeml
alignment = AlignIO.read('file.phy', 'phylip-sequential')

# Relaxed PHYLIP (no name limit) -- RAxML-NG, IQ-TREE (recommended default)
alignment = AlignIO.read('file.phy', 'phylip-relaxed')

# Always prefer phylip-relaxed for writing unless the downstream tool
# specifically requires strict format
AlignIO.write(alignment, 'output.phy', 'phylip-relaxed')

PHYLIP-Relaxed Dialect Mismatches Between Tree Tools

Biopython's 'phylip-relaxed' writes a single space between name and sequence. RAxML-NG and IQ-TREE accept this; PhyML rejects sequence names containing colons or parentheses; PAML's codeml expects sequential format with name-truncation behaviour distinct from interleaved. Common silent failures:

Always verify by running the downstream tool's "validate input only" mode (e.g. iqtree2 -s file.phy --check) before committing to a long compute.

MAF Block Coordinate Conventions

UCSC MAF (read via AlignIO.parse(file, 'maf')) returns blocks with per-row annotations:

• start (0-based; converts directly to BED but is off-by-one vs GFF)
• size (length on src strand)
• strand (+ or -)
• srcSize (length of source chromosome)

For minus-strand rows, start is measured from the END of the source contig: the corresponding plus-strand start is srcSize - start - size. Without this conversion, lifting MAF to genome coordinates places minus-strand blocks at the wrong locus. Reference: UCSC MAF spec at genome.ucsc.edu/FAQ/FAQformat.html#format5.

def maf_to_plus_strand_coords(row_anno):
    if row_anno['strand'] == '-':
        return row_anno['srcSize'] - row_anno['start'] - row_anno['size']
    return row_anno['start']

Stockholm Format Annotations

Stockholm format (used by Pfam, Rfam, HMMER) supports four annotation line types:

Common GC annotations: SS_cons (consensus secondary structure), RF (reference coordinates), seq_cons (consensus sequence).

WUSS notation in RNA #=GC SS_cons lines uses nested bracket pairs (<>, (), [], {}) for paired bases and characters like _, -, ,, :, ., ~ for unpaired regions; pseudoknots use upper/lower-case letter pairs (Aa, Bb). Consult the Infernal user guide for the full character table before writing or parsing custom SS_cons strings.

alignment = AlignIO.read('pfam.sto', 'stockholm')

for record in alignment:
    print(record.id, record.annotations)
    if 'secondary_structure' in record.letter_annotations:
        print(f'  SS: {record.letter_annotations["secondary_structure"]}')

ss_cons = alignment.column_annotations.get('secondary_structure')

Round-trip caveat: AlignIO.write(alignment, 'out.fasta', 'fasta') discards every Stockholm annotation silently. Re-reading and re-writing as Stockholm preserves GC/GR but dropped/added sequences invalidate the per-residue annotations -- regenerate annotations after edits.

Pfam-style name/start-end identifier convention: Pfam, Rfam, and Dfam Stockholm IDs (e.g. Q9Y6Y0/45-198) encode a 1-based inclusive region. Biopython does not split this; before passing to RAxML or IQ-TREE, parse the suffix into record.annotations['start'] / ['end'] and strip from record.id, then restore it after.