pavel-kravchenko Multiple Sequence Alignment
When to Use
- Aligning 3+ homologous protein or nucleotide sequences (gene families, orthologs) to find conserved motifs or domain boundaries.
- Choosing between ClustalW/MUSCLE/MAFFT/T-Coffee based on the number and length of input sequences.
- Building a guide tree for progressive alignment, or a quick alignment-free distance matrix.
- Computing alignment quality metrics (sum-of-pairs, column conservation) or a consensus/PSSM/sequence logo.
- Preparing an MSA as input for phylogenetics (
bio-core-phylogenetics) or protein family analysis.
Version Compatibility
- Biopython ≥ 1.80 (uses
Bio.Align.MultipleSeqAlignment,Bio.Align.AlignInfo,Bio.Align.substitution_matrices) - MAFFT ≥ 7.5, MUSCLE ≥ 5.1 (v5 CLI syntax:
-align/-output, not v3's-in/-out), Clustal Omega ≥ 1.2 - Python ≥ 3.9, NumPy ≥ 1.24
Prerequisites
pip install biopython numpy matplotlib- External aligners on PATH:
conda install -c bioconda mafft muscle clustalo(skill degrades to a pure-Python fallback if none are installed) - Familiarity with FASTA I/O (
bio-core-biopython-essentials) and pairwise alignment concepts (bio-core-pairwise-sequence-alignment)
Running an External Aligner and Loading the Result
Goal: Align a FASTA file of sequences with whichever MSA tool is available, then load it as a Biopython alignment object.
Approach: Try tools in order of speed/quality trade-off (MUSCLE/MAFFT are fast and scale well; Clustal Omega is a solid fallback); shell out with subprocess, then parse with AlignIO.
import subprocess
import shutil
from Bio import AlignIO
def run_msa_tool(tool_name, input_fasta, output_path, timeout=60):
"""Run an external MSA tool on a FASTA file. Returns True on success.
tool_name: one of 'mafft', 'muscle', 'clustalo'
input_fasta: path to unaligned FASTA
output_path: where to write the aligned FASTA
"""
commands = {
# MUSCLE v5 syntax (v3 used -in/-out -- check `muscle -version` first)
'muscle': ['muscle', '-align', input_fasta, '-output', output_path],
'mafft': ['mafft', '--auto', input_fasta],
'clustalo': ['clustalo', '-i', input_fasta, '-o', output_path,
'--outfmt=fasta', '--force'],
}
if tool_name not in commands:
raise ValueError(f"Unknown tool: {tool_name}")
if not shutil.which(commands[tool_name][0]):
print(f"{tool_name} not found. Install with: conda install -c bioconda {tool_name}")
return False
try:
result = subprocess.run(commands[tool_name], capture_output=True,
text=True, timeout=timeout, check=False)
if tool_name == 'mafft':
# mafft writes the alignment to stdout, not a file
with open(output_path, 'w') as f:
f.write(result.stdout)
return result.returncode == 0
except (subprocess.TimeoutExpired, FileNotFoundError) as e:
print(f"{tool_name} failed: {e}")
return False
def align_fasta(input_fasta, output_path, preferred=('muscle', 'mafft', 'clustalo')):
"""Try each tool in `preferred` order; return a Bio.Align.MultipleSeqAlignment or None."""
for tool in preferred:
if run_msa_tool(tool, input_fasta, output_path):
return AlignIO.read(output_path, 'fasta')
return None
Tool choice by dataset size: ClustalW is O(N^2), impractical beyond ~200 sequences. MUSCLE handles ~1K. MAFFT --auto scales to tens of thousands (--parttree beyond ~100K). T-Coffee gives the highest quality via library-based consistency but is too slow beyond ~500 sequences.
Alignment-Free Guide Tree (k-mer Distance + UPGMA)
Goal: Build a guide tree for progressive alignment without needing a full pairwise-alignment step, and produce a Newick string usable by phylogenetics tools.
Approach: Represent each sequence as a k-mer set, use Jaccard distance, then cluster with UPGMA (assumes a molecular clock — use neighbor-joining instead if rates vary across lineages).
import numpy as np
def kmer_distance(seq1, seq2, k=3):
"""Alignment-free Jaccard distance: fraction of k-mers unique to either sequence."""
kmers1 = set(seq1[i:i + k] for i in range(len(seq1) - k + 1))
kmers2 = set(seq2[i:i + k] for i in range(len(seq2) - k + 1))
if not kmers1 and not kmers2:
return 0.0
return 1.0 - len(kmers1 & kmers2) / len(kmers1 | kmers2)
def build_distance_matrix(sequences, k=3):
"""Pairwise k-mer distance matrix (symmetric, zero diagonal)."""
n = len(sequences)
matrix = np.zeros((n, n))
for i in range(n):
for j in range(i + 1, n):
d = kmer_distance(sequences[i], sequences[j], k)
matrix[i, j] = matrix[j, i] = d
return matrix
def upgma(dist_matrix, names):
"""UPGMA clustering -> Newick string. Assumes a constant molecular clock."""
n = len(names)
D = dist_matrix.copy()
node_names = list(names)
cluster_sizes = [1] * n
while len(node_names) > 1:
m = len(node_names)
min_dist, mi, mj = np.inf, 0, 1
for i in range(m):
for j in range(i + 1, m):
if D[i, j] < min_dist:
min_dist, mi, mj = D[i, j], i, j
bl = min_dist / 2
new_name = f"({node_names[mi]}:{bl:.4f},{node_names[mj]}:{bl:.4f})"
new_size = cluster_sizes[mi] + cluster_sizes[mj]
idx_map = [x for x in range(m) if x != mi and x != mj]
new_D = np.zeros((len(idx_map) + 1, len(idx_map) + 1))
for a, ka in enumerate(idx_map):
for b, kb in enumerate(idx_map):
new_D[a, b] = D[ka, kb]
new_D[a, -1] = new_D[-1, a] = (
cluster_sizes[mi] * D[ka, mi] + cluster_sizes[mj] * D[ka, mj]
) / new_size
node_names = [node_names[x] for x in idx_map] + [new_name]
cluster_sizes = [cluster_sizes[x] for x in idx_map] + [new_size]
D = new_D
return node_names[0] + ";"
Scoring an Alignment: Consensus, Conservation, Sum-of-Pairs
Goal: Quantify how good an MSA is and extract biological signal (consensus sequence, conserved blocks) from it.
Approach: Sum-of-pairs rewards every correct pair per column (used to compare candidate alignments); per-column conservation and information content identify functionally important, invariant residues.
from collections import Counter
import math
def sum_of_pairs_score(alignment, match=1, mismatch=-1, gap_penalty=-2):
"""Sum-of-Pairs score: for each column, sum the score of every sequence pair."""
n_seqs = len(alignment)
aln_length = len(alignment[0])
total_score, column_scores = 0, []
for pos in range(aln_length):
col_score = 0
for i in range(n_seqs):
for j in range(i + 1, n_seqs):
a, b = alignment[i][pos], alignment[j][pos]
if a == '-' or b == '-':
col_score += gap_penalty
elif a == b:
col_score += match
else:
col_score += mismatch
column_scores.append(col_score)
total_score += col_score
return total_score, column_scores
def compute_consensus(alignment, threshold=0.5):
"""Majority-rule consensus sequence plus per-position conservation fraction."""
n_seqs = len(alignment)
aln_length = len(alignment[0])
consensus, conservation = [], []
for pos in range(aln_length):
col = [alignment[s][pos] for s in range(n_seqs)]
counts = Counter(c for c in col if c != '-')
if not counts:
consensus.append('-')
conservation.append(0.0)
continue
char, count = counts.most_common(1)[0]
freq = count / n_seqs
conservation.append(freq)
consensus.append(char if freq >= threshold else 'x')
return ''.join(consensus), conservation
def information_content(alignment, alphabet_size=20):
"""Per-position Shannon information content in bits (max = log2(alphabet_size)).
Use alphabet_size=4 for DNA/RNA, 20 for protein.
"""
n_seqs = len(alignment)
max_entropy = math.log2(alphabet_size)
total_ic = []
for pos in range(len(alignment[0])):
col = [alignment[s][pos] for s in range(n_seqs)]
counts = Counter(c for c in col if c != '-')
total = sum(counts.values())
if total == 0:
total_ic.append(0.0)
continue
entropy = -sum((c / total) * math.log2(c / total) for c in counts.values())
total_ic.append(max_entropy - entropy)
return total_ic
Biopython shortcuts for the same job: Bio.Align.AlignInfo.SummaryInfo(alignment).dumb_consensus(threshold=0.5) and .pos_specific_score_matrix() for a PSSM.
Protein-Coding Genes: Align by Codon
## Back-translate a protein alignment to codons in R with seqinr,
## then export for downstream dN/dS analysis (e.g. PAML, HyPhy).
library(seqinr)
prot_aln <- read.alignment("protein_msa.fasta", format = "fasta")
nt_seqs <- read.fasta("nucleotide_cds.fasta")
## Map each non-gap protein column back to its 3-nt codon per sequence,
## inserting "---" for alignment gaps, to produce a codon-aware alignment.
Pitfalls
- Progressive alignment errors propagate: early misalignments are locked in. MAFFT iterative refinement (
--maxiterate) and T-Coffee's consistency library mitigate this. - Tool selection by dataset size: see the guidance above — using ClustalW on thousands of sequences will hang.
- Gap treatment: columns with >50% gaps are unreliable — mask with trimAl or Gblocks before phylogenetic analysis.
- SP score vs. Column Score: SP rewards each correct pair; Column Score requires all pairs correct in a column, making it far stricter for benchmarking.
- Protein-coding genes: align protein sequences, then back-translate to codons. Direct nucleotide alignment is misleading — synonymous substitutions saturate the 3rd codon position.
- UPGMA assumes a molecular clock: if evolutionary rates differ across lineages, use neighbor-joining instead for the guide tree.
- Exact MSA does not scale: the 3D dynamic-programming solution is O(n^3) for 3 sequences and O(n^k) for k sequences — only usable as a teaching demo (~10 residues), never for real data.
See Also
bio-core-pairwise-sequence-alignment— pairwise DP algorithms and substitution matrices underlying progressive MSAbio-core-phylogenetics— build and interpret trees from an MSA (guide trees vs. final phylogenies)bio-core-domains— identify conserved domains once sequences are alignedbio-core-biopython-essentials— FASTA I/O and Seq/SeqRecord basics for preparing MSA input