Writing Reordered and Reoriented FASTA¶

When you optimise a dotplot for collinearity, you often want to persist that layout as an actual FASTA file: contigs written in the collinearity-sorted order, with reverse-oriented contigs reverse-complemented so they read forward.

CrossIndex supports this via:

reorder_for_colinearity(query_group, target_group, reorder_target=...) — sort (and detect reverse orientation for) the query group; set reorder_target=False to leave the target group's order fixed.
reorder_by_length(group) — sort a group by descending length.
reversed_contigs(group) — the contigs detected as reverse-oriented.
write_fasta(path, group, order=..., reverse=...) — write a group to FASTA in the current (or a given) order, reverse-complementing the flagged contigs.

Orientation is relative. A contig is only "reversed" with respect to a reference. Throughout this tutorial Assembly A is the forward reference: its bases are never flipped; Assembly B is reoriented to match it.

This tutorial shows three scenarios.

Setup: two toy assemblies¶

Assembly A has three contigs (chrA, chrB, chrC). Assembly B contains the same three sequences but shuffled, with the copy of chrC reverse complemented — exactly the kind of mess a fresh assembly produces.

import random
import tempfile
from pathlib import Path

from dot_explorer.paf_io import CrossIndex, reverse_complement


def random_dna(length, seed):
    rng = random.Random(seed)
    return ''.join(rng.choice('ACGT') for _ in range(length))


# Assembly A: three distinct, non-repetitive contigs.
chrA = random_dna(400, seed=1)
chrB = random_dna(320, seed=2)
chrC = random_dna(260, seed=3)

assembly_a = {'chrA': chrA, 'chrB': chrB, 'chrC': chrC}

# Assembly B: same sequences, shuffled, with the chrC copy reverse-complemented.
assembly_b = {
    'utg1': chrB,                      # matches chrB (forward)
    'utg2': reverse_complement(chrC),  # matches chrC (reverse-oriented!)
    'utg3': chrA,                      # matches chrA (forward)
}

print('Assembly A order:', list(assembly_a))
print('Assembly B order:', list(assembly_b))

def build_cross():
    """Fresh CrossIndex with A in group 'a' and B in group 'b'."""
    cross = CrossIndex(k=13)
    for name, seq in assembly_a.items():
        cross.add_sequence(name, seq, group='a')
    for name, seq in assembly_b.items():
        cross.add_sequence(name, seq, group='b')
    # Matches with B as query and A as target (A is the reference axis).
    cross.compute_matches(query_group='b', target_group='a')
    return cross


def show_fasta(path):
    """Print FASTA headers and sequence lengths from a file."""
    for line in Path(path).read_text().splitlines():
        if line.startswith('>'):
            print(line)


tmp = Path(tempfile.mkdtemp())

Case 1 — Reorder + reorient B to match A, leaving A untouched¶

reorder_for_colinearity('b', 'a', reorder_target=False) sorts and orients B against A while keeping A exactly as loaded. We then write A verbatim and B in its new order with utg2 reverse-complemented.

cross = build_cross()
cross.reorder_for_colinearity('b', 'a', reorder_target=False)

print('A order (unchanged):', cross.contig_order['a'])
print('B order (optimised):', cross.contig_order['b'])
print('B reverse-oriented :', cross.reversed_contigs('b'))

cross.write_fasta(tmp / 'case1_A.fasta', 'a')  # forward reference, original order
cross.write_fasta(tmp / 'case1_B.fasta', 'b')  # reordered + reoriented

print('\n--- case1_A.fasta ---')
show_fasta(tmp / 'case1_A.fasta')
print('--- case1_B.fasta ---')
show_fasta(tmp / 'case1_B.fasta')

# The written utg2 is reverse-complemented, so it now equals the forward chrC.
written = {}
for line in (tmp / 'case1_B.fasta').read_text().split('>'):
    if not line.strip():
        continue
    head, *seq = line.splitlines()
    written[head.split()[0]] = ''.join(seq)
assert written['utg2'] == chrC
print('\nutg2 was flipped back to the forward orientation of chrC ✅')

Case 2 — A sorted by length, B reordered + reoriented against it¶

Here we first sort A by descending length (a deliberate, fixed layout), then align B to that fixed order. A's order changes but its strands do not.

cross = build_cross()
cross.reorder_by_length('a')                                  # A: fixed, by length
cross.reorder_for_colinearity('b', 'a', reorder_target=False)  # B: aligned to A

print('A order (by length):', cross.contig_order['a'])
print('B order (optimised):', cross.contig_order['b'])
print('B reverse-oriented :', cross.reversed_contigs('b'))

cross.write_fasta(tmp / 'case2_A.fasta', 'a')
cross.write_fasta(tmp / 'case2_B.fasta', 'b')
print('\n--- case2_A.fasta ---')
show_fasta(tmp / 'case2_A.fasta')
print('--- case2_B.fasta ---')
show_fasta(tmp / 'case2_B.fasta')

Case 3 — Both assemblies reordered (A is the reference frame)¶

With reorder_target=True (the default) both contig orders are optimised for collinearity. Orientation still needs a fixed frame, so A stays forward and B is reverse-complemented where needed. Both FASTA files are rewritten.

cross = build_cross()
cross.reorder_for_colinearity('b', 'a', reorder_target=True)

print('A order (optimised):', cross.contig_order['a'])
print('B order (optimised):', cross.contig_order['b'])
print('B reverse-oriented :', cross.reversed_contigs('b'))

cross.write_fasta(tmp / 'case3_A.fasta', 'a')
cross.write_fasta(tmp / 'case3_B.fasta', 'b')
print('\n--- case3_A.fasta ---')
show_fasta(tmp / 'case3_A.fasta')
print('--- case3_B.fasta ---')
show_fasta(tmp / 'case3_B.fasta')

Visual check¶

A dotplot of the case-3 layout, with B's reverse-oriented contigs rendered flipped, should show every block on the main diagonal. DotPlotter pulls the reversed set automatically when reverse_contigs is omitted.

import matplotlib.pyplot as plt

from dot_explorer.dotplot import DotPlotter

plotter = DotPlotter(cross)
fig = plotter.plot(
    query_group='b',
    target_group='a',
    scale_sequences=True,
    title='Case 3 — reordered + reoriented (B vs A)',
    dpi=100,
)
plt.close(fig)
print('Plotted collinear layout (reverse_contigs auto-detected).')