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

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((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.

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