PAF Workflow Tutorial¶

This notebook demonstrates the full PAF (Pairwise mApping Format) workflow in rusty-dot:

1. Generating PAF output from a SequenceIndex
2. Writing PAF output to a file
3. Reading PAF files back with parse_paf_file and PafAlignment
4. Parsing CIGAR strings for alignment statistics
5. Filtering and reordering contigs with PafAlignment.reorder_contigs
6. Using the low-level py_coords_to_paf function

PAF is a tab-separated, 12-column format used by many alignment tools. It is the primary output format of rusty-dot.

In [ ]:

import os
import tempfile

from rusty_dot import SequenceIndex, py_coords_to_paf
from rusty_dot.paf_io import (
    PafAlignment,
    PafRecord,
    parse_paf_file,
)

import os import tempfile from rusty_dot import SequenceIndex, py_coords_to_paf from rusty_dot.paf_io import ( PafAlignment, PafRecord, parse_paf_file, )

1. Generate PAF output from a SequenceIndex¶

The easiest way to produce PAF lines is SequenceIndex.get_paf.

In [ ]:

# Build a small index
idx = SequenceIndex(k=10)
unit = 'ACGTACGTACGT'
idx.add_sequence('query', unit * 8)  # 96 bp
idx.add_sequence('target', ('T' + unit * 7 + 'T'))  # 98 bp

# Get PAF lines
paf_lines = idx.get_paf('query', 'target', merge=True)
print(f'Generated {len(paf_lines)} PAF lines')
print('\nPAF columns (first record):')
col_names = [
    'query_name',
    'query_len',
    'query_start',
    'query_end',
    'strand',
    'target_name',
    'target_len',
    'target_start',
    'target_end',
    'residue_matches',
    'block_len',
    'mapq',
]
if paf_lines:
    for col, val in zip(col_names, paf_lines[0].split('\t')):
        print(f'  {col:20s}: {val}')

# Build a small index idx = SequenceIndex(k=10) unit = 'ACGTACGTACGT' idx.add_sequence('query', unit * 8) # 96 bp idx.add_sequence('target', ('T' + unit * 7 + 'T')) # 98 bp # Get PAF lines paf_lines = idx.get_paf('query', 'target', merge=True) print(f'Generated {len(paf_lines)} PAF lines') print('\nPAF columns (first record):') col_names = [ 'query_name', 'query_len', 'query_start', 'query_end', 'strand', 'target_name', 'target_len', 'target_start', 'target_end', 'residue_matches', 'block_len', 'mapq', ] if paf_lines: for col, val in zip(col_names, paf_lines[0].split('\t')): print(f' {col:20s}: {val}')

2. Write PAF to a file¶

PAF is a plain text format — just write each line with a newline.

In [ ]:

with tempfile.NamedTemporaryFile(mode='w', suffix='.paf', delete=False) as fh:
    paf_path = fh.name
    for line in paf_lines:
        fh.write(line + '\n')

print(f'PAF file written: {paf_path}')
print(f'File size: {os.path.getsize(paf_path)} bytes')

with tempfile.NamedTemporaryFile(mode='w', suffix='.paf', delete=False) as fh: paf_path = fh.name for line in paf_lines: fh.write(line + '\n') print(f'PAF file written: {paf_path}') print(f'File size: {os.path.getsize(paf_path)} bytes')

3. Read a PAF file with parse_paf_file¶

parse_paf_file is a generator that yields PafRecord objects one at a time. Comment lines (starting with #) and blank lines are skipped.

In [ ]:

records = list(parse_paf_file(paf_path))
print(f'Loaded {len(records)} records')

if records:
    rec = records[0]
    print('\nFirst record:')
    print(
        f'  query:   {rec.query_name}  [{rec.query_start}-{rec.query_end}]  len={rec.query_len}'
    )
    print(
        f'  target:  {rec.target_name}  [{rec.target_start}-{rec.target_end}]  len={rec.target_len}'
    )
    print(f'  strand:  {rec.strand}')
    print(f'  matches: {rec.residue_matches}')
    print(f'  q_aligned_len: {rec.query_aligned_len}')
    print(f'  t_aligned_len: {rec.target_aligned_len}')

records = list(parse_paf_file(paf_path)) print(f'Loaded {len(records)} records') if records: rec = records[0] print('\nFirst record:') print( f' query: {rec.query_name} [{rec.query_start}-{rec.query_end}] len={rec.query_len}' ) print( f' target: {rec.target_name} [{rec.target_start}-{rec.target_end}] len={rec.target_len}' ) print(f' strand: {rec.strand}') print(f' matches: {rec.residue_matches}') print(f' q_aligned_len: {rec.query_aligned_len}') print(f' t_aligned_len: {rec.target_aligned_len}')

4. Parsing a PAF record from a string¶

PafRecord.from_line parses a single PAF line string.

In [ ]:

line = 'query\t96\t5\t90\t+\ttarget\t98\t5\t90\t80\t85\t255\ttp:A:P\tcg:Z:10=5X20=3I30='
rec = PafRecord.from_line(line)

print('Parsed PafRecord:')
print(f'  cigar:            {rec.cigar}')
print(f'  alignment_length: {rec.alignment_length}')
print(f'  n_matches:        {rec.n_matches}')
print(f'  n_mismatches:     {rec.n_mismatches}')
print(f'  n_gaps:           {rec.n_gaps}')
print(f'  n_gap_bases:      {rec.n_gap_bases}')
print(f'  tags:             {rec.tags}')

line = 'query\t96\t5\t90\t+\ttarget\t98\t5\t90\t80\t85\t255\ttp:A:P\tcg:Z:10=5X20=3I30=' rec = PafRecord.from_line(line) print('Parsed PafRecord:') print(f' cigar: {rec.cigar}') print(f' alignment_length: {rec.alignment_length}') print(f' n_matches: {rec.n_matches}') print(f' n_mismatches: {rec.n_mismatches}') print(f' n_gaps: {rec.n_gaps}') print(f' n_gap_bases: {rec.n_gap_bases}') print(f' tags: {rec.tags}')

5. Serialise a PafRecord back to a string¶

to_line() recreates the 12-column PAF string (tags are not re-emitted).

In [ ]:

print(rec.to_line())

print(rec.to_line())

6. PafAlignment — collection with contig reordering¶

PafAlignment wraps a list of PafRecord objects and provides filtering and gravity-based contig reordering.

In [ ]:

# Load from the PAF file we created earlier
aln = PafAlignment.from_file(paf_path)
print(aln)
print(f'Query names:  {aln.query_names}')
print(f'Target names: {aln.target_names}')

# Load from the PAF file we created earlier aln = PafAlignment.from_file(paf_path) print(aln) print(f'Query names: {aln.query_names}') print(f'Target names: {aln.target_names}')

In [ ]:

# Construct from a list of records
aln2 = PafAlignment.from_records(records)
print(f'from_records: {aln2}')

# Construct from a list of records aln2 = PafAlignment.from_records(records) print(f'from_records: {aln2}')

Filtering¶

In [ ]:

# Filter to a specific query
filtered_q = aln.filter_by_query(['query'])
print(f"filter_by_query(['query']): {filtered_q}")

# Filter to a specific target
filtered_t = aln.filter_by_target(['target'])
print(f"filter_by_target(['target']): {filtered_t}")

# Filter by minimum alignment length (query span)
# Only keeps records where query_end - query_start >= min_length
filtered_len = aln.filter_by_min_length(20)
print(f'filter_by_min_length(20): {filtered_len}')

# Filter to a specific query filtered_q = aln.filter_by_query(['query']) print(f"filter_by_query(['query']): {filtered_q}") # Filter to a specific target filtered_t = aln.filter_by_target(['target']) print(f"filter_by_target(['target']): {filtered_t}") # Filter by minimum alignment length (query span) # Only keeps records where query_end - query_start >= min_length filtered_len = aln.filter_by_min_length(20) print(f'filter_by_min_length(20): {filtered_len}')

Gravity-based contig reordering¶

reorder_contigs sorts query and target names so that a subsequent dotplot shows maximum collinearity. Each contig is assigned a gravity equal to the weighted mean position of its alignment blocks on the opposing axis.

In [ ]:

# Build a multi-contig index to demonstrate reordering
idx_multi = SequenceIndex(k=8)
unit = 'ACGTACGTACGT'
idx_multi.add_sequence('qA', unit * 5)
idx_multi.add_sequence('qB', 'GCGCGCGCGCGCGCGCGCGC' * 3)
idx_multi.add_sequence('tA', unit * 5)
idx_multi.add_sequence('tB', ('T' + unit) * 4)

# Write all-vs-all PAF
with tempfile.NamedTemporaryFile(mode='w', suffix='.paf', delete=False) as fh:
    multi_paf_path = fh.name
    for q in ['qA', 'qB']:
        for t in ['tA', 'tB']:
            for line in idx_multi.get_paf(q, t):
                fh.write(line + '\n')

aln_multi = PafAlignment.from_file(multi_paf_path)
print(f'Multi-contig alignment: {aln_multi}')

q_sorted, t_sorted = aln_multi.reorder_contigs(
    query_names=['qA', 'qB'],
    target_names=['tA', 'tB'],
)
print(f'Reordered queries:  {q_sorted}')
print(f'Reordered targets:  {t_sorted}')

# Build a multi-contig index to demonstrate reordering idx_multi = SequenceIndex(k=8) unit = 'ACGTACGTACGT' idx_multi.add_sequence('qA', unit * 5) idx_multi.add_sequence('qB', 'GCGCGCGCGCGCGCGCGCGC' * 3) idx_multi.add_sequence('tA', unit * 5) idx_multi.add_sequence('tB', ('T' + unit) * 4) # Write all-vs-all PAF with tempfile.NamedTemporaryFile(mode='w', suffix='.paf', delete=False) as fh: multi_paf_path = fh.name for q in ['qA', 'qB']: for t in ['tA', 'tB']: for line in idx_multi.get_paf(q, t): fh.write(line + '\n') aln_multi = PafAlignment.from_file(multi_paf_path) print(f'Multi-contig alignment: {aln_multi}') q_sorted, t_sorted = aln_multi.reorder_contigs( query_names=['qA', 'qB'], target_names=['tA', 'tB'], ) print(f'Reordered queries: {q_sorted}') print(f'Reordered targets: {t_sorted}')

7. Low-level: py_coords_to_paf¶

py_coords_to_paf converts raw coordinate tuples to PAF strings without going through a SequenceIndex. It is useful when you have coordinate data from another source.

In [ ]:

matches = [
    (0, 50, 10, 60),
    (55, 90, 65, 100),
]

paf_strings = py_coords_to_paf(
    matches,
    query_name='my_query',
    query_len=100,
    target_name='my_target',
    target_len=120,
)

for line in paf_strings:
    print(line)

matches = [ (0, 50, 10, 60), (55, 90, 65, 100), ] paf_strings = py_coords_to_paf( matches, query_name='my_query', query_len=100, target_name='my_target', target_len=120, ) for line in paf_strings: print(line)

8. Using SequenceIndex.optimal_contig_order¶

The SequenceIndex also exposes a built-in gravity-based ordering via optimal_contig_order, which uses the cached comparison results.

In [ ]:

# Pre-compute to populate the cache
idx_multi.precompute_all_pairs(merge=True)

q_order, t_order = idx_multi.optimal_contig_order(
    query_names=['qA', 'qB'],
    target_names=['tA', 'tB'],
)
print(f'Optimal query order:  {q_order}')
print(f'Optimal target order: {t_order}')

# Pre-compute to populate the cache idx_multi.precompute_all_pairs(merge=True) q_order, t_order = idx_multi.optimal_contig_order( query_names=['qA', 'qB'], target_names=['tA', 'tB'], ) print(f'Optimal query order: {q_order}') print(f'Optimal target order: {t_order}')

Summary¶

Task	API
Generate PAF lines	SequenceIndex.get_paf(query, target)
Parse a PAF file	parse_paf_file(path) or PafAlignment.from_file(path)
Parse a single line	PafRecord.from_line(line)
Serialise a record	PafRecord.to_line()
Filter by query/target	PafAlignment.filter_by_query(names) / filter_by_target(names)
Filter by alignment length	PafAlignment.filter_by_min_length(min_length)
Reorder contigs	PafAlignment.reorder_contigs() or SequenceIndex.optimal_contig_order()
Convert coords to PAF	py_coords_to_paf(matches, ...)