| name | bio-core-comparative-genomics |
| description | 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. |
| tool_type | python |
| primary_tool | MUMmer |
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.
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(((seq2), (seq1)))
rev = np.zeros(((seq2), (seq1)))
seq2_rc = rc(seq2)
i ((seq2) - window + ):
j ((seq1) - window + ):
fm = (seq1[j + k] == seq2[i + k] k (window))
fm >= threshold:
fwd[i + window // , j + window // ] = fm / window
rm = (seq1[j + k] == seq2_rc[i + k] k (window))
rm >= threshold:
rev[(seq2) - - i - window // , j + window // ] = rm / window
fwd, rev
random_dna = n: .join(np.random.choice((), n))
seg_A, seg_B, seg_C, seg_D = [random_dna() _ ()]
comp = {: , : , : , : }
rc = s: .join(comp.get(c, c) c (s))
reference = seg_A + seg_B + seg_C + seg_D
rearranged = seg_A + rc(seg_C) + seg_B + seg_D
fwd, rev = dna_dotplot(reference, rearranged, window=, stringency=)
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.
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.
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
