# Bio Core Multiple Sequence Alignment

> Align FASTA sequences via MAFFT/MUSCLE/Clustal Omega, parse with Biopython AlignIO, build UPGMA guide trees, score sum-of-pairs/consensus/conservation. Use for MSA, MSA tool choice, consensus/logo, or guide-tree building.

- Skill: `pavel-kravchenko/bio-core-multiple-sequence-alignment` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-multiple-sequence-alignment`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-multiple-sequence-alignment/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Design & Media
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-core-multiple-sequence-alignment

---


# 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`.

```python
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).

```python
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.

```python
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

```r
## 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 MSA
- `bio-core-phylogenetics` — build and interpret trees from an MSA (guide trees vs. final phylogenies)
- `bio-core-domains` — identify conserved domains once sequences are aligned
- `bio-core-biopython-essentials` — FASTA I/O and Seq/SeqRecord basics for preparing MSA input