# Bio Core Comparative Genomics

> Build dot plots, detect synteny/rearrangements, classify orthologs vs paralogs, compute pan-genomes; pick MUMmer/LASTZ/minimap2 for alignment. Use when comparing genomes, reading a dot plot, or finding synteny/orthologs.

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

---


# Comparative Genomics: Dot Plots, Synteny, and Orthology

## When to Use
- Visualizing sequence similarity and structural rearrangements between two sequences or genomes (dot plots)
- Detecting conserved gene order (synteny blocks) between species, or where synteny breaks
- Classifying a pair of genes as orthologs (speciation) vs. paralogs (duplication)
- Computing a bacterial pan-genome (core / accessory / unique genes) across strains
- Choosing a whole-genome alignment tool (MUMmer, LASTZ, minimap2, Cactus) by divergence level

## Version Compatibility
Python ≥3.10, NumPy ≥1.24, matplotlib ≥3.7. External tool references: MUMmer4 (nucmer/promer), LASTZ 1.04, minimap2 ≥2.26, OrthoFinder ≥2.5, fastANI ≥1.34, Progressive Mauve. The dot-plot/synteny/pan-genome functions below are pure Python and have no external dependency beyond NumPy.

## Prerequisites
- `pip install numpy matplotlib`
- Optional CLI tools for real genome-scale work: `mummer`, `lastz`, `minimap2` (install via bioconda)
- Helpful prior skills: bio-sequence-io-read-sequences (loading FASTA), bio-alignment-pairwise-alignment (scoring models)

## Key Concepts

**Ortholog vs. paralog**: Orthologs diverged via speciation (same function across species). Paralogs diverged via duplication (may gain new functions). Reliable ortholog inference requires sequence similarity + syntenic context, not identity alone.

**Dot plot diagonals encode structure**:
- Main diagonal (bottom-left → top-right) = conserved linear segment
- Anti-diagonal = inverted segment (reverse-complement match)
- Off-main-diagonal block = transposition
- Repeated parallel diagonals = tandem duplication/repeat

**Core genome vs. pan-genome**: Core = genes in all strains. Pan = core + accessory (some strains) + unique (single strain). Bacterial pan-genomes are "open" (grow with more sequencing); eukaryotic pan-genomes are "closed".

**Identity thresholds are heuristic**: >90% ANI = same species; >70% identity usable with LASTZ; <30% identity = "twilight zone" where synteny is required to call orthology. ANI < 95% across strains generally implies different species.

## Tool Selection

| Scenario | Tool |
|----------|------|
| Closely related genomes (>90% ANI), fast | MUMmer/nucmer |
| Moderate divergence (>70% identity) | LASTZ |
| Any pairwise comparison, long reads | minimap2 |
| Multiple sequence alignment | MAFFT, MUSCLE |
| Synteny visualization | MCScanX, SynMap, JCVI, Genomicus |
| Multi-genome alignment with rearrangements | Cactus, Progressive Mauve |

**Goal:** Visualize local sequence similarity and detect structural rearrangements (inversions, translocations, duplications) between two sequences.
**Approach:** Slide a window across both sequences, score the fraction of matching positions per window, and separately test the reverse complement so inversions surface as anti-diagonals instead of gaps.

```python
import numpy as np

def filtered_dotplot(seq1, seq2, window=5, stringency=0.6):
    """Sliding-window dot plot. Returns a 2D match-score matrix (len(seq2) x len(seq1))."""
    threshold = int(window * stringency)
    matrix = np.zeros((len(seq2), len(seq1)))
    for i in range(len(seq2) - window + 1):
        for j in range(len(seq1) - window + 1):
            matches = sum(seq1[j + k] == seq2[i + k] for k in range(window))
            if matches >= threshold:
                matrix[i + window // 2, j + window // 2] = matches / window
    return matrix

def dna_dotplot(seq1, seq2, window=11, stringency=0.7):
    """Dot plot with reverse-complement detection. Returns (forward, revcomp) matrices."""
    comp = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'}
    rc = lambda s: ''.join(comp.get(c, c) for c in reversed(s))
    threshold = int(window * stringency)
    fwd = np.zeros((len(seq2), len(seq1)))
    rev = np.zeros((len(seq2), len(seq1)))
    seq2_rc = rc(seq2)
    for i in range(len(seq2) - window + 1):
        for j in range(len(seq1) - window + 1):
            fm = sum(seq1[j + k] == seq2[i + k] for k in range(window))
            if fm >= threshold:
                fwd[i + window // 2, j + window // 2] = fm / window
            rm = sum(seq1[j + k] == seq2_rc[i + k] for k in range(window))
            if rm >= threshold:
                rev[len(seq2) - 1 - i - window // 2, j + window // 2] = rm / window
    return fwd, rev

random_dna = lambda n: ''.join(np.random.choice(list('ACGT'), n))  # example: a reference A-B-C-D, C inverted+shifted
seg_A, seg_B, seg_C, seg_D = [random_dna(150) for _ in range(4)]
comp = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'}
rc = lambda s: ''.join(comp.get(c, c) for c in reversed(s))
reference = seg_A + seg_B + seg_C + seg_D          # A-B-C-D
rearranged = seg_A + rc(seg_C) + seg_B + seg_D     # A-C'-B-D (inversion + transposition)
fwd, rev = dna_dotplot(reference, rearranged, window=11, stringency=0.7)
```

`fwd` shows a diagonal for A (top-left) and D (bottom-right) with an off-diagonal block for B; `rev` shows an anti-diagonal for the inverted C segment.

**Goal:** Detect conserved gene order (synteny blocks) between two genomes from a list of matched genes.
**Approach:** Pair genes by ortholog identity, then walk consecutive pairs and start a new block whenever the relative order or strand relationship flips — a proxy for detecting inversions/transpositions.

```python
def find_ortholog_pairs(genes_a, genes_b):
    """Match genes by name (proxy for ortholog identification). Returns list of (gene_a, gene_b) tuples."""
    name_to_b = {g['name']: g for g in genes_b}
    return [(ga, name_to_b[ga['name']]) for ga in genes_a if ga['name'] in name_to_b]

def detect_synteny_blocks(pairs):
    """Split ortholog pairs into collinear/inverted runs based on order + strand relationship."""
    blocks, current = [], [pairs[0]]
    for i in range(1, len(pairs)):
        prev_a, prev_b = pairs[i - 1]
        curr_a, curr_b = pairs[i]
        same_order = curr_b['start'] > prev_b['start']
        same_strand_rel = ((curr_a['strand'] == curr_b['strand']) ==
                           (prev_a['strand'] == prev_b['strand']))
        if same_order and same_strand_rel:
            current.append(pairs[i])
        else:
            blocks.append(current)
            current = [pairs[i]]
    blocks.append(current)
    return blocks
```

Synteny block interpretation: a **collinear block** (genes in same order/relative strand) is a conserved region; an **inverted block** (same genes, reversed order, flipped strands) marks an inversion event; **orphan pairs** outside any block indicate a transposition or gene loss/gain.

**Goal:** Classify genes across bacterial strains into core (all strains), accessory (some strains), and unique (one strain) sets.
**Approach:** Union all strains' gene sets, intersect across every strain for the core, and subtract the rest to find genes unique to one strain.

```python
def pan_genome_analysis(strains):
    """
    strains: dict of strain_name -> {'genes': set(gene_ids), ...}
    Returns dict with 'core', 'accessory' (sets), 'unique' (dict strain->set), and size counters.
    """
    all_genes = set()
    for s in strains.values():
        all_genes |= s['genes']

    core = all_genes.copy()
    for s in strains.values():
        core &= s['genes']

    unique = {}
    for name, s in strains.items():
        others = set()
        for other_name, other_s in strains.items():
            if other_name != name:
                others |= other_s['genes']
        unique[name] = s['genes'] - others

    unique_all = set().union(*unique.values()) if unique else set()
    accessory = all_genes - core - unique_all
    return {
        'core': core, 'accessory': accessory, 'unique': unique,
        'core_size': len(core), 'accessory_size': len(accessory),
    }
```

## Pitfalls

- **Dot plot window size trade-off**: small windows (1–3) show all matches including noise; large windows (>15) miss short conserved regions. Start with window=11, stringency=0.7 for genomic DNA.
- **Self-comparison diagonal is uninformative**: mask the main diagonal ± window/2 to reveal internal repeats instead of the trivial identity line.
- **Ortholog thresholds are heuristic**: >30% identity for bacterial orthologs, >70% for strain-level, but synteny + phylogeny (or an orthogroup tool like OrthoFinder) are required for confident calls.
- **ANI vs. percent identity**: ANI averages over all shared genomic regions; BLAST percent identity is local to an alignment. Use ANI (fastANI/PyANI) for species-level comparisons, not raw BLAST identity.
- **Inversions only appear as anti-diagonals when the reverse complement is checked**: a dot plot using only forward-strand matches shows a gap, not a diagonal, for inverted regions.
- **Pan-genome openness**: never extrapolate a bacterial core genome from <10 strains — the core shrinks as more diverse strains are added, and open pan-genomes never plateau.
- **Coordinate systems**: MUMmer output is 1-based; BED/BEDTools are 0-based half-open. Always check a tool's coordinate convention before merging intervals downstream.

## See Also
- bio-comparative-genomics-synteny-analysis
- bio-comparative-genomics-ortholog-inference
- bio-comparative-genomics-ancestral-reconstruction
- bio-phylogenetics-modern-tree-inference