Dotplot Visualization Tutorial¶

import os
import tempfile

import matplotlib.pyplot as plt

from dot_explorer import SequenceIndex
from dot_explorer.dotplot import DotPlotter

# Helper to create a reverse complement
def revcomp(seq):
    table = str.maketrans('ACGTacgt', 'TGCAtgca')
    return seq.translate(table)[::-1]


unit = 'ACGTACGTACGT'  # 12 bp repeat unit
seq_a = unit * 10  # 120 bp  — the reference
seq_b = 'T' + unit * 9 + 'T'  # 120 bp  — shifted by 1
seq_c = revcomp(unit * 5) + unit * 5  # 120 bp  — half inverted

idx = SequenceIndex(k=8)
idx.add_sequence('reference', seq_a)
idx.add_sequence('shifted', seq_b)
idx.add_sequence('partial_inv', seq_c)

print(f'Index: {idx}')

plotter = DotPlotter(idx)

# No output_path: the figure is returned and displayed inline in Jupyter
fig = plotter.plot(title='All vs All — inline display')
plt.close(fig)  # free memory when no longer needed

# Inline display for a single pair
fig = plotter.plot_single(
    query_name='reference',
    target_name='partial_inv',
    title='reference vs partial_inv — inline',
)
plt.close(fig)

with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    all_vs_all_path = fh.name

fig = plotter.plot(
    output_path=all_vs_all_path,
    title='All vs All',
)
plt.close(fig)
print(f'Saved: {all_vs_all_path}  ({os.path.getsize(all_vs_all_path)} bytes)')

with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    subset_path = fh.name

plotter.plot(
    query_names=['reference', 'shifted'],
    target_names=['partial_inv'],
    output_path=subset_path,
    title='Reference & Shifted vs Partial Inversion',
)
print(f'Subset plot saved: {subset_path}')

with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    single_path = fh.name

plotter.plot_single(
    query_name='reference',
    target_name='shifted',
    output_path=single_path,
    figsize=(5, 5),
    title='reference vs shifted',
)
print(f'Single-pair plot saved: {single_path}')

with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    styled_path = fh.name

plotter.plot(
    output_path=styled_path,
    dot_size=1.5,
    cap_style='round',
    dot_color='crimson',
    dpi=200,
    title='Custom style: crimson, round caps, dpi=200',
)
print(f'Styled plot saved: {styled_path}')

with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    unmerged_path = fh.name

plotter.plot_single(
    query_name='reference',
    target_name='shifted',
    output_path=unmerged_path,
    merge=False,
    title='reference vs shifted (unmerged k-mer hits)',
)
print(f'Unmerged plot saved: {unmerged_path}')

for dpi in [72, 150, 300]:
    with tempfile.NamedTemporaryFile(suffix=f'_dpi{dpi}.png', delete=False) as fh:
        path = fh.name
    plotter.plot_single(
        'reference',
        'shifted',
        output_path=path,
        dpi=dpi,
        title=f'DPI = {dpi}',
    )
    size_kb = os.path.getsize(path) / 1024
    print(f'DPI={dpi:4d}  file size={size_kb:.1f} kB  path={path}')

with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    large_path = fh.name

plotter.plot(
    output_path=large_path,
    figsize_per_panel=6.0,  # each panel is 6×6 inches
    title='Large panels (6 inches each)',
)
print(f'Large-panel plot saved: {large_path}')

# PNG (default raster format)
with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    png_path = fh.name

fig = plotter.plot(output_path=png_path, title='PNG output')
plt.close(fig)
print(f'PNG:  {png_path}  ({os.path.getsize(png_path)} bytes)')

# SVG via file extension
with tempfile.NamedTemporaryFile(suffix='.svg', delete=False) as fh:
    svg_ext_path = fh.name

fig = plotter.plot(output_path=svg_ext_path, title='SVG via extension')
plt.close(fig)
print(f'SVG (ext):     {svg_ext_path}  ({os.path.getsize(svg_ext_path)} bytes)')

# SVG via explicit format parameter (output path need not end in .svg)
with tempfile.NamedTemporaryFile(suffix='.out', delete=False) as fh:
    svg_fmt_path = fh.name

fig = plotter.plot(output_path=svg_fmt_path, format='svg', title='SVG via format param')
plt.close(fig)
print(f'SVG (fmt):     {svg_fmt_path}  ({os.path.getsize(svg_fmt_path)} bytes)')

# PDF — vector format suitable for LaTeX / Word
with tempfile.NamedTemporaryFile(suffix='.pdf', delete=False) as fh:
    pdf_path = fh.name

fig = plotter.plot(output_path=pdf_path, title='PDF output')
plt.close(fig)
print(f'PDF:  {pdf_path}  ({os.path.getsize(pdf_path)} bytes)')

# Verify SVG header
with open(svg_ext_path) as f:
    snippet = f.read(80)
print(f'\nSVG file header: {snippet!r}')

# Without filtering: all merged hits are drawn
with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    unfiltered_path = fh.name
plotter.plot_single(
    'reference', 'shifted', output_path=unfiltered_path, title='No min_length filter'
)

# With filtering: only hits of at least 24 bp are drawn
with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    filtered_path = fh.name
plotter.plot_single(
    'reference',
    'shifted',
    output_path=filtered_path,
    min_length=24,
    title='min_length=24',
)

print(f'Unfiltered: {unfiltered_path}  ({os.path.getsize(unfiltered_path)} bytes)')
print(f'Filtered:   {filtered_path}  ({os.path.getsize(filtered_path)} bytes)')

import random

from dot_explorer.paf_io import PafAlignment, PafRecord

# Build synthetic PAF records with varying identity values to illustrate
# the colour-by-identity feature without needing a real aligner.
random.seed(0)
paf_records = []
for i in range(8):
    block_len = random.randint(10, 30)
    identity = 0.6 + 0.04 * i  # 60 % … 88 %
    residue_matches = round(block_len * identity)
    q_start = i * 12
    paf_records.append(
        PafRecord(
            query_name='reference',
            query_len=144,
            query_start=q_start,
            query_end=q_start + block_len,
            strand='+',
            target_name='shifted',
            target_len=144,
            target_start=q_start,
            target_end=q_start + block_len,
            residue_matches=residue_matches,
            alignment_block_len=block_len,
            mapping_quality=255,
        )
    )

paf_aln = PafAlignment(paf_records)
print(f'{len(paf_aln)} records loaded')

# Create a DotPlotter with the PAF alignment attached
identity_plotter = DotPlotter(idx, paf_alignment=paf_aln)

# ── Plot coloured by identity ─────────────────────────────────────────────
with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    identity_path = fh.name

fig = identity_plotter.plot_single(
    'reference',
    'shifted',
    output_path=identity_path,
    color_by_identity=True,
    identity_palette='viridis',
    title='Coloured by identity (viridis)',
)
plt.close(fig)
print(f'Identity plot: {identity_path}  ({os.path.getsize(identity_path)} bytes)')

# ── Colorbar (standalone scale figure) ───────────────────────────────────
with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    colorbar_path = fh.name

fig_cb = identity_plotter.plot_identity_colorbar(
    palette='viridis',
    output_path=colorbar_path,
)
plt.close(fig_cb)
print(f'Colorbar:      {colorbar_path}  ({os.path.getsize(colorbar_path)} bytes)')

# ── Try a different palette ───────────────────────────────────────────────
with tempfile.NamedTemporaryFile(suffix='.png', delete=False) as fh:
    plasma_path = fh.name

fig_p = identity_plotter.plot_single(
    'reference',
    'shifted',
    output_path=plasma_path,
    color_by_identity=True,
    identity_palette='plasma',
    title='Coloured by identity (plasma)',
)
plt.close(fig_p)
print(f'Plasma plot:   {plasma_path}  ({os.path.getsize(plasma_path)} bytes)')

import random

# Reference and a 3%-divergent copy.  SNPs break the diagonal into many short
# exact-match k-mer blocks, producing a segment-rich panel.
random.seed(0)
bases = 'ACGT'
ref = ''.join(random.choice(bases) for _ in range(8000))
mut = ''.join(random.choice(bases) if random.random() < 0.03 else b for b in ref)

div_idx = SequenceIndex(k=15)
div_idx.add_sequence('ref', ref)
div_idx.add_sequence('mut', mut)
div_plotter = DotPlotter(div_idx)


def count_segments(fig):
    """Count drawn match segments across all panels of a figure."""
    from matplotlib.collections import LineCollection

    return sum(
        len(c.get_segments())
        for ax in fig.axes
        for c in ax.collections
        if isinstance(c, LineCollection)
    )


def svg_is_vector(path):
    """True when the SVG has no embedded raster <image> (match layer is vector)."""
    with open(path) as fh:
        return '<image' not in fh.read()


# Same data, three rasterization policies (plus an 'auto' with a low threshold
# to show the automatic switch to a rasterised layer).
cases = [
    ('auto (default)', {}),
    ('force vector', {'rasterized': False}),
    ('force raster', {'rasterized': True}),
    ('auto, thresh=100', {'rasterization_threshold': 100}),
]
for label, kwargs in cases:
    with tempfile.NamedTemporaryFile(suffix='.svg', delete=False) as fh:
        path = fh.name
    fig = div_plotter.plot_single('ref', 'mut', output_path=path, **kwargs)
    n = count_segments(fig)
    plt.close(fig)
    kb = os.path.getsize(path) / 1024
    kind = 'vector' if svg_is_vector(path) else 'raster'
    print(f'{label:18s}: {n:6d} segments  {kb:8.1f} kB  {kind}')

# Sweep chain_gap; keep rasterized=False so file size tracks the segment count.
for gap in [0, 25, 100, 1000]:
    with tempfile.NamedTemporaryFile(suffix='.svg', delete=False) as fh:
        path = fh.name
    fig = div_plotter.plot_single(
        'ref',
        'mut',
        output_path=path,
        chain_gap=gap,
        rasterized=False,
        title=f'chain_gap = {gap}',
    )
    n = count_segments(fig)
    plt.close(fig)
    kb = os.path.getsize(path) / 1024
    print(f'chain_gap={gap:5d}: {n:6d} segments  {kb:8.1f} kB (vector)')

# Build a small demo index and a GFF3 annotation.
import random

from dot_explorer import DotPlotter, SequenceIndex
from dot_explorer.annotation import GffAnnotation

random.seed(11)
seq = "".join(random.choice("ACGT") for _ in range(4000))
ann_idx = SequenceIndex(k=13)
ann_idx.add_sequence("chrA", seq)
ann_idx.add_sequence("chrB", seq[2000:] + seq[:2000])

gff_text = """\
chrA\tdemo\tgene\t201\t1400\t.\t+\t.\tID=gene1;Name=GeneA
chrA\tdemo\tCDS\t201\t600\t.\t+\t0\tID=cds1;Parent=gene1
chrA\tdemo\tCDS\t901\t1400\t.\t+\t0\tID=cds1;Parent=gene1
chrA\tdemo\tgene\t1201\t2400\t.\t-\t.\tID=gene2;Name=GeneB
chrA\tdemo\trepeat_region\t2801\t3600\t.\t.\t.\tID=rep1
chrB\tdemo\tgene\t501\t1800\t.\t-\t.\tID=gene3;Name=GeneC
"""
ann = GffAnnotation.from_text(gff_text)
ann_plotter = DotPlotter(ann_idx)

# Diagonal squares: features shade self-vs-self panels *behind* the
# alignments (one colour per feature type).
fig = ann_plotter.plot(annotation=ann, title="Annotated all-vs-all grid")

# Focused single-pair view with side annotation tracks: lane-packed
# features left of the y axis and below the x axis, direction arrows for
# stranded types and connectors joining multi-part CDS groups.
fig = ann_plotter.plot(
    query_names=["chrB"],
    target_names=["chrA"],
    annotation_query=ann,
    annotation_target=ann,
    annotation_tracks=True,
)