Skip to content

All-vs-All Dotplot Between Two Sequence Indices

This tutorial demonstrates how to compare sequences from two separate FASTA files (e.g. two genome assemblies) using CrossIndex. Each FASTA file is assigned to its own group, and pairwise alignments are computed between all sequences in group A and all sequences in group B.

Overview

  1. Create a CrossIndex and load two FASTA files into groups A and B.
  2. Compute cross-group k-mer matches with compute_matches().
  3. Retrieve PAF alignments with get_paf() and optionally export them.
  4. Build a PafAlignment from those records for contig reordering.
  5. Plot the sorted all-vs-all dotplot with DotPlotter.
from dot_explorer import SequenceIndex
from dot_explorer.dotplot import DotPlotter
from dot_explorer.paf_io import CrossIndex, PafAlignment, PafRecord

1. Create example sequences

For demonstration we build two small sets of sequences in memory. In real usage you would call cross.load_fasta('genome_a.fasta', group='a') and cross.load_fasta('genome_b.fasta', group='b').

# Sequences for "genome A" — three contigs of different lengths
genome_a = {
    'contigA1': 'ACGTACGTACGTACGTACGT' * 10,  # 200 bp
    'contigA2': 'TACGTACGTACGTACGTACG' * 5,  # 100 bp
    'contigA3': 'GCGCGCGCGCGCGCGCGCGC' * 3,  # 60 bp
}

# Sequences for "genome B" — three contigs
genome_b = {
    'contigB1': 'ACGTACGTACGTACGTACGT' * 8,  # 160 bp  (similar to contigA1)
    'contigB2': 'GCGCGCGCGCGCGCGCGCGC' * 4,  # 80 bp   (similar to contigA3)
    'contigB3': 'TTTTAAAAAGGGGCCCCTTTT' * 2,  # 42 bp   (no matches)
}

2. Build a CrossIndex

CrossIndex holds one internal SequenceIndex that contains sequences from both groups, using internal prefixes (a: / b:) to prevent name collisions.

cross = CrossIndex(k=10)

for name, seq in genome_a.items():
    cross.add_sequence(name, seq, group='a')

for name, seq in genome_b.items():
    cross.add_sequence(name, seq, group='b')

print(cross)

CrossIndex (k=10)
  Total sequences : 6
  Group 'a'         :      3 sequences
  Group 'b'         :      3 sequences
  Computed pairs  : none (call compute_matches() first)
  PAF records     : 0

3. Compute cross-group k-mer matches

Call compute_matches() before retrieving PAF lines or reordering contigs. This is the explicit computation step — matches are not calculated automatically when loading sequences.

cross.compute_matches()  # compute k-mer matches between group 'a' and group 'b'

print('Computed group pairs:', cross.computed_group_pairs)

paf_lines = cross.get_paf(merge=True)

print(f'Total PAF lines: {len(paf_lines)}')
for line in paf_lines[:5]:
    print(line)

Computed group pairs: [('a', 'b')]
Total PAF lines: 11108
contigA1    200 188 200 +   contigB1    160 0   12  12  12  255
contigA1    200 184 200 +   contigB1    160 0   16  16  16  255
contigA1    200 180 200 +   contigB1    160 0   20  20  20  255
contigA1    200 176 200 +   contigB1    160 0   24  24  24  255
contigA1    200 172 200 +   contigB1    160 0   28  28  28  255

4. Build a PafAlignment for contig reordering

Parse the raw PAF strings into PafRecord objects so we can use the reorder_contigs method to maximise collinearity.

records = [PafRecord.from_line(line) for line in paf_lines]
aln = PafAlignment.from_records(records)

q_sorted, t_sorted = aln.reorder_contigs(
    query_names=cross.query_names,
    target_names=cross.target_names,
)

print('Sorted query (genome A) contigs:', q_sorted)
print('Sorted target (genome B) contigs:', t_sorted)

Sorted query (genome A) contigs: ['contigA1', 'contigA2', 'contigA3']
Sorted target (genome B) contigs: ['contigB1', 'contigB2', 'contigB3']

Note: Contigs with no cross-group matches (e.g. contigB3) are placed at the end, sorted by descending length.

5. Build a combined index for plotting (optional)

Note: You can plot directly from a CrossIndex using query_group and target_group parameters — see step 6 — and DotPlotter also accepts a PafAlignment directly (DotPlotter(aln)), which renders straight from the alignment records. The combined index below is only one more alternative, useful when you want plain (unprefixed) sequence names without groups.

combined_idx = SequenceIndex(k=10)

for name, seq in genome_a.items():
    combined_idx.add_sequence(name, seq)

for name, seq in genome_b.items():
    combined_idx.add_sequence(name, seq)

print(combined_idx)

SequenceIndex(k=10, sequences=6)

6. Plot the all-vs-all dotplot with CrossIndex directly

DotPlotter accepts query_group and target_group parameters so you can plot directly from a CrossIndex without building a separate combined index. Pre-computed merged alignments from compute_matches() are used automatically.

By default rows and columns render in load order. To display the collinearity-sorted layout, pass contig_order='colinearity' — the plot-time ordering runs the same gravity algorithm as reorder_contigs above.

Pass scale_sequences=True so that each subplot's width and height are proportional to the lengths of the compared sequences.

import matplotlib.pyplot as plt

# Plot directly from CrossIndex — no separate combined_idx needed.
plotter = DotPlotter(cross)

fig = plotter.plot(
    query_group='a',     # sequence names looked up from group 'a'
    target_group='b',    # sequence names looked up from group 'b'
    contig_order='colinearity',  # sort rows/columns by the gravity algorithm
    figsize_per_panel=4.0,
    scale_sequences=True,
    title='Genome A vs Genome B — collinearity-sorted contigs',
    dpi=100,
)
plt.close(fig)  # free memory when no longer needed

Using CrossIndex.reorder_contigs directly

CrossIndex.reorder_contigs() runs the same gravity algorithm and returns the sorted (query, target) name lists; it also records which query contigs are reverse-oriented (exposed via reversed_contigs). It does not change the CrossIndex's stored contig order — pass the returned lists to plot(query_names=…), use contig_order='colinearity' at plot time, or call reorder_for_colinearity() to apply the ordering to the index in place (as the reorder/reorient tutorial does). compute_matches() must have been called first — it is the explicit step that computes and caches the cross-group matches.

# compute_matches was already called above.  reorder_contigs returns the
# sorted orders (and records orientation flags) without mutating the index.
q_opt, t_opt = cross.reorder_contigs()
print('Optimal query order:', q_opt)
print('Optimal target order:', t_opt)

Optimal query order: ['contigA1', 'contigA2', 'contigA3']
Optimal target order: ['contigB1', 'contigB2', 'contigB3']

7. Flip reverse-oriented contigs

reorder_contigs() also detects query contigs that align in reverse orientation against their best-matching target chromosome (using the d-genies orientation check). These are exposed via cross.reversed_contigs(group).

Passing that set to DotPlotter.plot(reverse_contigs=...) renders those contigs reverse-complemented so their alignments read along the main diagonal instead of the anti-diagonal — the underlying coordinates are never modified, only the rendering. Omitting the argument auto-pulls the detected set for query_group.

# Query contigs detected as reverse-oriented against their best chromosome
reversed_a = cross.reversed_contigs('a')
print('Reverse-oriented query contigs:', reversed_a)

# Render them flipped so they read along the main diagonal.  (Omit the
# reverse_contigs argument to auto-pull this same set for query_group='a'.)
fig = plotter.plot(
    query_group='a',
    target_group='b',
    figsize_per_panel=4.0,
    scale_sequences=True,
    title='Reverse-oriented contigs flipped onto the main diagonal',
    dpi=100,
    reverse_contigs=reversed_a,
)
plt.close(fig)

Reverse-oriented query contigs: {'contigA1', 'contigA3'}