API reference

Library entry point

Library-facing API.

A stable, import-friendly surface for embedding arda in other Python tools. The heavy lifting lands in Phase 2 (arda.annotate); this module keeps the public signature stable.

arda.adapter.annotate_sequences(sequences, seqtype='nt', organism='human', map_d=True)

Annotate FR/CDR regions for a batch of sequences.

Parameters:

• sequences (Iterable[str] | Iterable[tuple[str, str]]) – Either raw sequence strings or (id, sequence) pairs.

• seqtype (Literal['nt', 'aa']) – "nt" for nucleotide input, "aa" for amino acid.

• organism (str) – One of the supported organisms (human, mouse, rat, rabbit, rhesus_monkey).

• map_d (bool) – Map D segments (d_call/d2_call/np*) for VDJ-locus hits; False skips D mapping. Applies to nucleotide input only.

Returns:

A list of AIRR-style annotation record dicts (one per input sequence).

Runtime annotation

Runtime annotation: map input sequences to the reference and transfer markup.

Pipeline: read input (FASTA/FASTQ) -> MMseqs2 search against the curated scaffold DB -> best hit per query -> project reference region markup onto the query (C++ hot path) -> AIRR TSV.

arda.annotate.mapper.annotate_file(input, output, organism='human', seqtype='nt', *, threads=0, sensitivity=None, strand='both', chunk_size=50000, map_d=True, d_max_evalue=None, shm='framework', complete_junction_nt=0)

Annotate a FASTA/FASTQ file and stream an AIRR TSV.

The input is processed in bounded chunks with a background reader thread that prefetches the next chunk while the current one is annotated (mmseqs releases the GIL during its subprocess), so memory stays flat for arbitrarily large FASTQ and read parsing overlaps compute. The reference + target DB are loaded once and reused across all chunks.

Parameters:

• input (str | Path)

• output (str | Path)

• organism (str)

• seqtype (str)

• threads (int)

• sensitivity (float | None)

• strand (str)

• chunk_size (int)

• map_d (bool)

• d_max_evalue (float | None)

• shm (str)

• complete_junction_nt (int)

Return type:

Path

arda.annotate.mapper.annotate_records(records, organism='human', seqtype='nt', *, threads=0, sensitivity=None, strand='both', map_d=True, d_max_evalue=None, shm='framework', complete_junction_nt=0)

Annotate in-memory (id, sequence) records; return AIRR record dicts.

Parameters:

• strand (str) – "both" (default, nt only) searches both strands and re-orients reverse-complement hits; "forward" searches the plus strand only. Ignored for protein input.

• map_d (bool) – True (default) maps D segments into the junction of VDJ-locus hits (d_call/d2_call/np*); False skips D mapping (nt input only — D mapping never runs for protein input).

• d_max_evalue (float | None) – E-value gate on the D call(s); None keeps the shipped 0.2. Lower is stricter (see arda.annotate.transfer._map_d()).

• shm (str) – SHM scoping — "framework" (default), "both" or "off"; see arda.shm.

• records (list[tuple[str, str]])

• organism (str)

• seqtype (str)

• threads (int)

• sensitivity (float | None)

• complete_junction_nt (int)

Return type:

list[dict]

arda.annotate.mapper.build_index(organism='all', *, force=False)

(Re)build the precompiled mmseqs DBs shipped under database/.

Writes database/vdj/<org>/mmseqs/<seqtype>/db* + a VERSION marker so the runtime can use them out of the box (and detect a mmseqs-version mismatch). Skips up-to-date DBs unless force.

Parameters:

• organism (str)

• force (bool)

Return type:

None

Project reference region markup onto a query via the C++ hot path.

Takes a parsed mmseqs hit plus the reference entry for the matched scaffold and returns an AIRR-style record dict for the query.

Junction handling follows AIRR strictly: junction spans the conserved Cys104 through the [FW]118 that opens FR4; junction_aa starts with C and ends with F/W for a canonical rearrangement. A junction is reported even when not canonical (out-of-frame, missing the conserved residues). For an out-of-frame junction (V and J in different frames) the amino-acid translation inserts 1-2 N bases after the V germline end to restore the J frame; the codon that then contains an inserted N is rendered as _. The V/J split inside the junction is located from the transferred v_sequence_end / j_sequence_start.

arda.annotate.transfer.transfer_hit(query_id, query_seq, hit, ref, seqtype='nt', rev_comp=False, d_germlines=None, submitted_seq=None, anchors=None, d_max_evalue=None, shm='framework', complete_junction_nt=0)

Build an AIRR record by projecting ref region coords onto the query.

query_seq is the coding-strand sequence all markup/coords/CIGARs are computed on. submitted_seq is the read AS SUBMITTED, stored verbatim in the AIRR sequence field; for a reverse-strand hit it is the reverse complement of query_seq and rev_comp is set, per AIRR (“if rev_comp is True, all output data are based on the reverse complement of sequence”). Defaults to query_seq (forward reads, where the two are identical).

d_max_evalue overrides the D-call E-value gate; see _map_d().

shm scopes the SHM fields — framework (default), both or off; see arda.shm.

complete_junction_nt > 0 lets a read that reached Cys104 but stopped before [FW]118 have its junction finished from the called J’s germline, up to that many nt; the count lands in junction_completed_nt. 0 (the default) emits observed junctions only. See _germline_completed_junction().

Parameters:

• query_id (str)

• query_seq (str)

• hit (dict)

• ref (RefEntry)

• seqtype (str)

• rev_comp (bool)

• d_germlines (list[tuple[str, str]] | None)

• submitted_seq (str | None)

• anchors (dict | None)

• d_max_evalue (float | None)

• shm (str)

• complete_junction_nt (int)

Return type:

dict

Load the curated reference markup for runtime projection.

For nucleotide annotation we use markup.tsv + alleles.fasta; for amino acid annotation markup.aa.tsv + alleles.aa.fasta. Both expose region *_start/*_end columns in the same coordinate space (nt or aa), so the projection code is identical.

class arda.annotate.reference.RefEntry(locus, v_call, j_call, starts, ends, v_sequence_end=0, j_sequence_start=0, c_call='', vj_end=0)

Bases: object

Per-scaffold reference markup: region coords (in target space) + calls.

Parameters:

• locus (str)

• v_call (str)

• j_call (str)

• starts (list[int])

• ends (list[int])

• v_sequence_end (int)

• j_sequence_start (int)

• c_call (str)

• vj_end (int)

locus: str

v_call: str

j_call: str

starts: list[int]

ends: list[int]

v_sequence_end: int

j_sequence_start: int

c_call: str

vj_end: int

property is_jc: bool

a J followed by the CH1 exon, with no V.

Type:

A constant-region scaffold

class arda.annotate.reference.Reference(organism, seqtype, target_fasta, entries, d_germlines, anchors=<factory>, _jc_combos=None)

Bases: object

In-memory reference for one (organism, seqtype).

Parameters:

• organism (str)

• seqtype (str)

• target_fasta (Path)

• entries (dict[str, RefEntry])

• d_germlines (dict[str, list[tuple[str, str]]])

• anchors (dict)

• _jc_combos (dict[tuple[str, str], str] | None)

organism: str

seqtype: str

target_fasta: Path

entries: dict[str, RefEntry]

d_germlines: dict[str, list[tuple[str, str]]]

anchors: dict

jc_combinations()

(j_allele, c_allele) -> J+C scaffold id, the C-side twin of combinations.tsv.

The segment reference carries the J and the constant region as SEPARATE targets, because a constant sequence shared across a locus’ J+C scaffolds is a cross-product and copying it through cost 76.4 % of the segment search’s alignments. So a J→C read now names its home the same way a V→J read does — by the pair it hit, resolved through this table.

Derived from the loaded markup rather than from a file: the J+C scaffolds are already there, and a second generated artifact is one more thing that can go stale against it.

⛔ Keyed per ALLELE, never by the comma-joined group string. A j_call is a group of alleles arda could not tell apart, and the two sides of this lookup group them by different rules: a J+C scaffold’s j_call collapses alleles with an identical J sequence (refbuild.constant), while a J| segment target inherits the V×J scaffold’s j_call, which collapses alleles with an identical assembled V + pad + J and folds in reading frame (refbuild.combinations). They genuinely disagree: human IGL yields J|IGLJ2*01,IGLJ3*01 and J|IGLJ2A*01 on one side against IGLJ2*01,IGLJ2A*01,IGLJ3*01 on the other, so a group-string key leaves 24 J+C scaffolds unreachable from any J hit — and an unreachable scaffold means the J→C contest never fires for those reads, which is exactly the invented-junction/destroyed-c_call bug the contest exists to prevent. Splitting both sides removes the need for the two collapses to agree at all.

⛔ Segment targets are excluded by their `|` prefix. load_reference loads segments.markup.tsv into this same key space after markup.tsv, so a pre-2.8.0 file – whose 345 rows are J+C scaffolds copied verbatim as JC|<sid> with the same j_call and c_call – would collide with each base scaffold and, being later, win. Every value would become an id the full target DB does not contain, _align_implied would drop them all as unknown targets, and the contest would be silently off for the whole run: no error, plausible output. Reachable on the first run after any upgrade, because the self-heal regenerates the files without reloading entries.

Return type:

dict[tuple[str, str], str]

segment_j_call(name)

J allele for a JC| segment target, which is named by SCAFFOLD id, not by allele.

Feeding the raw target name into a (V, J) combination lookup silently fails for every J->C read; measured, that collapsed the two-pass fast path from 85.3 % to 0.1 %.

Retained for references built before the constant region became its own C| target — a build and a mapper of different vintages must not silently mis-resolve a J call.

Parameters:

name (str)

Return type:

str

get(scaffold_id)

Parameters:

scaffold_id (str)

Return type:

RefEntry | None

arda.annotate.reference.load_reference(organism, seqtype='nt')

Load reference markup + target FASTA path for an organism.

Parameters:

• organism (str)

• seqtype (str)

Return type:

Reference

Sequence I/O: streaming FASTA/FASTQ readers and chunking.

Native parsing (no BioPython). Transparently handles gzip by .gz extension.

arda.annotate.io.open_text(path)

Open a (possibly gzipped) text file for reading.

Parameters:

path (str | Path)

arda.annotate.io.read_sequences(path, *, with_qual=False)

Yield (id, sequence) from a FASTA or FASTQ file (auto-detected).

with_qual=True yields (id, sequence, qual) instead – the FASTQ Phred string, or None for FASTA (which has no quality). The default is unchanged and pays nothing: the quality line is consumed either way, only kept when asked (needed solely by the paired overlap-merge, arda.rnaseq.map.merge_pair()).

Parameters:

• path (str | Path)

• with_qual (bool)

Return type:

Iterator[tuple]

arda.annotate.io.detect_format(path)

Return "fasta" or "fastq" by peeking at the first non-empty char.

Parameters:

path (str | Path)

Return type:

str

arda.annotate.io.write_fasta(records, path)

Write (id, sequence) records to a FASTA file.

Parameters:

• records (Iterator[tuple[str, str]])

• path (str | Path)

Return type:

Path

arda.annotate.io.chunked(it, size)

Yield lists of up to size items from an iterator.

Parameters:

• it (Iterator)

• size (int)

Return type:

Iterator[list]

Per-segment AIRR CIGAR strings from the mmseqs scaffold alignment.

arda aligns a query to a V + N*pad + J [+ C] scaffold, not to each germline segment separately. AIRR wants a CIGAR per segment (v_cigar/j_cigar/c_cigar) whose reference is that segment’s germline; since the scaffold’s V part IS the V germline (target position == germline position), its J part IS the J allele, and its C part the CH1 exon, each segment’s CIGAR is the sub-walk of the one query->scaffold alignment whose target falls in that segment’s range.

CIGAR operators follow the AIRR spec (SAM subset):

• S – query positions before the alignment starts (query 5’ offset). Required, precedes N.

• N – reference positions before the alignment starts (germline 5’ offset). Required.

• M / I (gap in reference) / D (gap in query) – the aligned body.

• trailing S (query 3’ remainder) is emitted; trailing N (germline 3’ remainder) is optional per the spec and omitted (arda does not always know the full germline length, e.g. the C-region CH1 exon is longer than the shipped stub).

segment_cigars builds all three in a SINGLE pass over the aligned strings.

SOMATIC HYPERMUTATION (v_mutations / j_mutations). The same walk answers “which germline positions does this read differ at”, which is the input an SHM or lineage-tree tool needs, so it is emitted from here rather than re-derived downstream: G45A,C112T — germline base, 1-based position in that segment’s OWN germline, read base.

The information was already recoverable — sequence_alignment and germline_alignment carry every column, and on a real bulk IG library the germline they report matches the shipped allele on 28,365 of 28,365 mapped reads (66,526 V mismatches, zero disagreements). What it was not, was usable: recovering it needs the scaffold geometry (mmseqs2_tstart, mmseqs2_t_vend, mmseqs2_t_jstart, mmseqs2_t_vjend) and the knowledge that arda aligns to a V + pad + J scaffold rather than to a germline. A consumer that does the obvious thing instead — diff the two AIRR alignment strings — gets 100,091 mismatches on that library of which 20,140 (20.1 %) are N-pad or constant-region columns, i.e. it attributes junction positions to a germline.

⛔ Which is the reason the scoping here is STRUCTURAL and not a filter. A mutation inside the V..J interior is not attributable to any germline: V(D)J recombination chews the segment ends back and adds non-templated N/P bases, so the V-end / NDN / J-start partition of a junction is often not identifiable from the sequence at all. The mutation lists are built only for seg_key V and J — the pad is not a segment, so a pad column has no germline coordinate to be recorded under and cannot enter the list by any code path. Nothing downstream needs to remember to exclude it.

Substitutions only. Indels are in the CIGAR (I/D) and stay there: an SHM indel is one event of unbounded length, not a per-position call, and its representation would have to be re-interpreted by every consumer. The germline coordinates on the far side of an indel are still correct, because the walk tracks the target position across the gap columns.

Correcting cigars for CONTIGS (Stage 3). A contig is just a long query, so BOTH ways to get its cigars end in segment_cigars and produce the same record (see arda.annotate.contig):

• RE-ANNOTATE the assembled contig through mapper.annotate_records – one mmseqs alignment, then segment_cigars. No cigar arithmetic; check_cigar validates it.

• MERGE the reads’ existing alignments column-by-column into the contig’s (C++ _markup.merge_alignment), skipping the alignment pass.

Both are built and proven byte-for-byte equal (tests/unit/test_contig_merge.py on 29 real GenBank contigs). MEASURED (arda-benchmark scripts/bench_contig_cigars.py): at scRNA-seq scale (~10^5 contigs/sample) merge is ~9x faster – the whole gap is mmseqs; the C++ stitch is ~3 % of merge’s wall and barely grows with read depth. Prefer merge when the assembly layout is available (the reads carry their scaffold + offset); re-annotate is the fallback when it is not.

arda.annotate.cigar.parse_cigar(cigar)

"57S291M1054S" -> [(57,"S"), (291,"M"), (1054,"S")]. Inverse of build_cigar().

Parameters:

cigar (str)

Return type:

list[tuple[int, str]]

arda.annotate.cigar.cigar_query_length(cigar)

Query (read/contig) bases the CIGAR spans – M/I/S/=/X; D and N are reference-side.

Parameters:

cigar (str)

Return type:

int

arda.annotate.cigar.cigar_reference_length(cigar)

Reference (germline) bases the CIGAR spans – M/D/N/=/X; I and S are query-side.

Parameters:

cigar (str)

Return type:

int

arda.annotate.cigar.check_cigar(cigar, query_len)

A CIGAR is consistent with a query of query_len iff its query-side ops sum to it.

This is the invariant a corrected/re-annotated sequence (a read OR an assembled contig – a contig is just a long query) must satisfy: v_cigar/j_cigar/c_cigar each lay over the WHOLE sequence, soft-clipping the parts outside their own segment. Use it to validate a cigar after correcting or re-deriving it.

Parameters:

• cigar (str)

• query_len (int)

Return type:

bool

arda.annotate.cigar.build_cigar(q_lead, g_lead, ops, q_trail)

Assemble one AIRR CIGAR: {q_lead}S {g_lead}N <body> {q_trail}S (parts of length 0 are dropped). ops is the per-column M/I/D sequence of the aligned body; consecutive equal operators are run-length encoded. Trailing germline N is intentionally omitted (optional).

Parameters:

• q_lead (int)

• g_lead (int)

• ops (list[str])

• q_trail (int)

Return type:

str

arda.annotate.cigar.segment_cigars(qaln, taln, qstart, tstart, qlen, t_vend, t_jstart, t_vjend)

Per-segment AIRR fields from ONE walk of the alignment: v_cigar/j_cigar/c_cigar and v_mutations/j_mutations (only the keys that have content).

qaln/taln are the mmseqs aligned strings (- for gaps), qstart/tstart their 1-based query/target start, qlen the full query length. Boundaries are 1-based scaffold positions; pass 0 for an absent segment.

The mutation lists are G45A,C112T — germline base, 1-based position in that segment’s own germline, read base — and they are what makes the SHM in an arda record usable without re-deriving it. See arda.annotate.cigar’s module docstring for why the segment scoping is structural rather than a filter.

Parameters:

• qaln (str)

• taln (str)

• qstart (int)

• tstart (int)

• qlen (int)

• t_vend (int)

• t_jstart (int)

• t_vjend (int)

Return type:

dict[str, str]

D-segment mapping on a bare nucleotide junction — no read, no mmseqs search.

transfer._map_d already maps D into the V..J interior of a query, but it is fed v_sequence_end / j_sequence_start projected from an mmseqs scaffold hit. A VDJdb-style record has no read to align: it has a junction and a V/J call. The per-allele germlines shipped in database/vdj/<org>/cdr3_anchors.tsv close that gap, so the interior can be located directly and the existing mapper reused.

Junction space, as everywhere in arda.cdr3fix: the input runs Cys104 -> Phe/Trp118 inclusive.

Finding the interior. The V germline is exact at the junction’s 5’ end (V/D/J are not somatically mutated in TCR, and IGH mutation is rare this close to the anchor), so the V contribution is the longest common prefix of the junction and the V’s CDR3-region germline; the J contribution is the longest common suffix. Validated against OLGA ground truth on 1300 junctions across human IGH/TRB/TRD and mouse TRB: the prefix length is exact for 80-85 % of records and never underestimates (it can overshoot by 1-2 nt when the first N-region base happens to match), and the derived interior contains the whole true D segment in 94-99 % of records.

Only IGH, TRB and TRD have D germlines; VJ loci return an empty call.

class arda.annotate.dmap.DCall(locus='', d_call='', d_sequence_start=-1, d_sequence_end=-1, d_support='', d2_call='', d2_sequence_start=-1, d2_sequence_end=-1, d2_support='', np1='', np2='', np3='', v_sequence_end=-1, j_sequence_start=-1, extra=<factory>)

Bases: object

D mapping of one junction. Coordinates are 1-based closed, junction space.

Parameters:

• locus (str)

• d_call (str)

• d_sequence_start (int)

• d_sequence_end (int)

• d_support (str)

• d2_call (str)

• d2_sequence_start (int)

• d2_sequence_end (int)

• d2_support (str)

• np1 (str)

• np2 (str)

• np3 (str)

• v_sequence_end (int)

• j_sequence_start (int)

• extra (dict)

locus: str = ''

d_call: str = ''

d_sequence_start: int = -1

d_sequence_end: int = -1

d_support: str = ''

d2_call: str = ''

d2_sequence_start: int = -1

d2_sequence_end: int = -1

d2_support: str = ''

np1: str = ''

np2: str = ''

np3: str = ''

v_sequence_end: int = -1

j_sequence_start: int = -1

extra: dict

property called: bool

property is_dd: bool

markup(junction_nt)

The junction cut into labelled parts, 5’->3’.

[("V", ...), ("np1", ...), (d_call, ...), ("np2", ...), (d2_call, ...), ("np3", ...), ("J", ...)] for a tandem D-D, without the last two entries for a single D, and [("V", ...), ("N", ...), ("J", ...)] when no D was called.

The parts concatenate back to junction_nt EXACTLY – that is the contract this method exists to make checkable, and it is what a D-D markup consumer needs: the AIRR columns alone give it as a set of coordinates and three np strings that it must re-derive the D-observed sequence from.

⛔ The V-end / np / D-start boundaries INSIDE the junction are not identifiable from sequence – exonuclease chew-back and N/P addition make the partition probabilistic. This is one consistent reading of the junction, not ground truth. Empty when the V/J split could not be located at all.

Parameters:

junction_nt (str)

Return type:

list[tuple[str, str]]

arda.annotate.dmap.map_d_junction(junction_nt, v_call, j_call, species='human', d_max_evalue=None)

Map D (and a tandem second D) into a bare nucleotide junction.

d_max_evalue overrides the shipped E-value gate on the D call(s); see arda.annotate.transfer._map_d(). None keeps the shipped 0.2.

Parameters:

• junction_nt (str)

• v_call (str)

• j_call (str)

• species (str)

• d_max_evalue (float | None)

Return type:

DCall

Two ways to give an assembled contig its AIRR cigars — and they agree.

A Stage-3 contig is a consensus of reads that Stage 1 already aligned to a scaffold. Its v_cigar/j_cigar/c_cigar + alignment strings can be produced two ways:

• reannotate_contigs() — treat the contig as one long query and run it back through annotate_records() (one mmseqs alignment, then segment_cigars). Simple, exact, no new code; the cost is a second alignment pass.

• merge_contigs() — stitch the reads’ existing alignments into the contig’s (the C++ _markup.merge_alignment per-column consensus over N reads), skipping the alignment pass. Wins when a sample has ~10^5 contigs (scRNA-seq).

Both converge on the same synthetic hit and reuse transfer_hit(), so their output is field-for-field comparable. Which is optimal is a measured question — see tests/unit/test_contig_merge.py and the arda-benchmark Phase-D benchmark.

class arda.annotate.contig.ReadPlacement(qaln, taln, qstart, tstart, offset)

Bases: object

One read’s placement in a contig: its scaffold alignment + contig offset.

qaln/taln are the read’s coding-strand aligned strings vs the scaffold (- for gaps, as Stage 1 emits). qstart/tstart are 1-based starts in the read / scaffold. offset is the 0-based position of the read’s first base within the contig (the assembly layout), in contig orientation.

Parameters:

• qaln (str)

• taln (str)

• qstart (int)

• tstart (int)

• offset (int)

qaln: str

taln: str

qstart: int

tstart: int

offset: int

class arda.annotate.contig.Contig(sequence_id, sequence, target, reads=<factory>)

Bases: object

An assembled contig + the reads it was built from, all hitting one scaffold.

Parameters:

• sequence_id (str)

• sequence (str)

• target (str)

• reads (list[ReadPlacement])

sequence_id: str

sequence: str

target: str

reads: list[ReadPlacement]

arda.annotate.contig.reannotate_contigs(records, organism='human', seqtype='nt', *, threads=0, sensitivity=None, strand='both', map_d=True, d_max_evalue=None)

Annotate assembled contigs by re-aligning them (baseline path).

records are (contig_id, contig_seq). A thin wrapper over annotate_records(): a contig is just a long query.

Parameters:

• records (list[tuple[str, str]])

• organism (str)

• seqtype (str)

• threads (int)

• sensitivity (float | None)

• strand (str)

• map_d (bool)

• d_max_evalue (float | None)

Return type:

list[dict]

arda.annotate.contig.merge_contig(contig, reference, *, map_d=True)

Annotate one contig by stitching its reads’ alignments (merge path).

reference is a preloaded Reference (load it once for a whole sample; see merge_contigs()). Raises KeyError if the contig’s target scaffold is absent from the reference.

Parameters:

• contig (Contig)

• reference (Reference)

• map_d (bool)

Return type:

dict

arda.annotate.contig.merge_contigs(contigs, organism='human', seqtype='nt', *, reference=None, map_d=True)

Annotate contigs by the merge path; loads the reference once for all of them.

Parameters:

• contigs (list[Contig])

• organism (str)

• seqtype (str)

• reference (Reference | None)

• map_d (bool)

Return type:

list[dict]

Segment shortlist → implied V×J scaffold, with a rescue set so no read is ever lost.

The fast path: align against the 1,244-target segment reference (V, J and J+C separately, see arda.refbuild.segments), take each read’s best V and best J, and look the pair up in combinations.tsv. That names exactly one V×J scaffold in the full reference — so the second alignment is one target per read instead of ~277, measured at 5.36 s → 0.044 s (122×) on a 20 k-read TRA amplicon, and it lands back in the scaffold coordinate system arda’s markup transfer already speaks.

The invariant this module exists to enforce. A fast path that silently drops reads is not an optimisation, it is a different tool: arda’s whole claim is near-zero Stage-1 false negatives. So every read is accounted for. shortlist() partitions the input into

• implied — a V and a J both hit, and the pair names a real scaffold. The fast path.

• rescue — anything else: only a V hit, only a J hit, an unknown V×J pair, or a read whose second-pass alignment failed. These go back to the full reference, exactly as today.

implied ∪ rescue == every read that hit anything, by construction and by assertion. The rescue set is small and cheap to realign — measured 1.9 % of amplicon reads and ~10 % of bulk, i.e. 11 % and 20 % of the new total cost — so exactness costs almost nothing:

amplicon fast 0.84 s + rescue 0.10 s = 0.94 s vs 5.36 s → 5.7× bulk fast 0.52 s + rescue 0.13 s = 0.65 s vs 1.26 s → 2.0×

Why reads land in rescue, measured on a TRA amplicon:

• V only (12.6 %) — the read never reaches a J, so no pair exists. The baseline picks an arbitrary J for these; the rescue pass reproduces that behaviour exactly rather than inventing a different arbitrary answer.

• J only (2.1 %) — J→C reads with no V.

• failed second alignment (1.9 %) — the synthesized diagonal is derived from the V hit, so it is wrong for a read whose scaffold alignment does not begin in V.

α/δ is not ambiguity, it is the answer. TRD is TRAV/DV + TRDJ: the J (and C) decides the locus. An earlier draft rescued TRAV/DV reads whose best J crossed the locus, which discarded real rearrangements — arda’s reference already carries 45 TRAV/DV + TRDJ scaffolds under locus TRD alongside 1,005 TRAV/DV + TRAJ ones under TRA. So the “114 of 212 residual V disagreements are TRA→TRD” finding is most likely the two-pass being more correct than a baseline that picks by whole-scaffold bit score with an arbitrary J half.

class arda.annotate.shortlist.Shortlist(implied=<factory>, rescue=<factory>, reasons=<factory>, reason_of=<factory>)

Bases: object

Which reads take the fast path, which must be realigned, and why.

Parameters:

• implied (dict[str, str])

• rescue (list[str])

• reasons (dict[str, int])

• reason_of (dict[str, str])

implied: dict[str, str]

read id -> scaffold id implied by (best V, best J)

rescue: list[str]

read ids that must go back to the full reference

reasons: dict[str, int]

reason -> count, for the run report

reason_of: dict[str, str]

a caller can only re-route one rescue class (e.g. v_only onto its own V segment) if it knows which reads are in it.

Type:

read id -> why it was rescued. The same information as reasons, per read

property n_total: int

property fast_fraction: float

as_dict()

Return type:

dict

arda.annotate.shortlist.load_combinations(path)

combinations.tsv → {(v_allele, j_allele): scaffold_id}.

v_calls/j_calls may be comma-separated ambiguity lists; every listed allele maps to the scaffold, so a lookup succeeds whichever member the segment pass reported.

Parameters:

path (str | Path)

Return type:

dict[tuple[str, str], str]

arda.annotate.shortlist.shortlist(best_v, best_j, combos, *, failed=None)

Partition reads into the fast path and the rescue set.

Parameters:

• best_v (dict[str, str]) – read id -> best V| allele (absent if the read hit no V target).

• best_j (dict[str, str]) – read id -> best J|/JC| allele.

• combos (dict[tuple[str, str], str]) – from load_combinations().

• failed (set[str] | None) – read ids whose second-pass alignment produced nothing; always rescued.

Returns:

Shortlist. Every read appearing in best_v or best_j lands in exactly one of implied / rescue — asserted, not assumed.

Return type:

Shortlist

⚠ best_j must hold a J allele. JC| targets are named by scaffold id, not by allele (see arda.refbuild.segments), so a caller feeding the raw target name straight through gets no_such_combination for every J→C read. Measured cost of getting this wrong: the fast path collapsed from 85.3 % to 0.1 %, with 9,388 reads needlessly rescued — correct output, silently no faster. Resolve via segments.markup.tsv’s j_call column.

Exact k-mer prefilter — drop reads that cannot align before MMseqs2 sees them.

On bulk RNA-seq, mmseqs search spends nearly all of its time proving that reads are not receptor reads: 4 M reads of SRR10611239 take 48.9 s to find 947 hits (0.024 %). The fitted cost model says why — wall ~ reads/46,353 + hits/350, so the dominant term is set by the read count and not by the answer. A read can only align to a V(D)J scaffold if it shares an exact k-mer with one; testing that is a lookup, proving it is Smith-Waterman.

The index is built from target_fasta — the same FASTA MMseqs2 searches. That is deliberate: the design’s largest single finding is that a prefilter built over V+pad+J alone loses 16.29 % of real reads, 69.27 % of them J->C reads, and that indexing the constant region takes the loss to 0.53 % (OPTIMIZATION.md §3.3). Deriving the index from a hand-listed set of segments is exactly how that hole would come back; deriving it from the search target makes it structurally impossible for the two to disagree about what is indexable.

The filter is off by default. It trades a measured ~0.5 % of real reads for ~6x on bulk, and arda’s near-zero Stage-1 false-negative rate is the thing it is not allowed to trade silently.

arda.prefilter.K = 16

Seed length. k=12 passes 62-65 % of reads (no filtering); k>=18 adds nothing over 16.

arda.prefilter.MIN_HITS = 1

Windows a read must share with the reference to survive. A real read sharing one exact 16-mer usually shares many, so >=2 barely moves the pass rate (4.64 % -> 3.91 %) while FN climbs.

arda.prefilter.MAX_USEFUL_PASS_RATE = 0.3

Above this pass rate the filter costs more than it saves — amplicon libraries run 46-49 % receptor, where MMseqs2 has to look at nearly every read anyway.

arda.prefilter.build(target_fasta, k=16)

Index every k-mer of target_fasta and its reverse complement.

Cached per FASTA: the index is ~1-3 MB and takes well under a second to build, but a per-chunk rebuild would charge every chunk of a 4 M-read run for it.

Parameters:

• target_fasta (Path)

• k (int)

arda.prefilter.available()

Is the native extension importable? A source tree without a built ext is a normal state.

Return type:

bool

arda.prefilter.keep_records(records, target_fasta, *, threads=1, min_hits=1, k=16)

The (id, sequence) records worth handing to MMseqs2.

Filtering happens entirely inside the extension. Returning a per-read mask instead would put two full Python passes around the scan – building a list of sequences to pass in, unpacking a list of ints coming back – and those cost more than the scan does: the same C++ change measured 2.66x in isolation and 1.16x through a mask-shaped API. Here the only Python objects created are the survivors, which on bulk is 0.5-2 % of the input.

Parameters:

• records (list[tuple[str, str]])

• target_fasta (Path)

• threads (int)

• min_hits (int)

• k (int)

Return type:

list[tuple[str, str]]

Junction markup and repair

Working from a bare (junction_aa, v_call, j_call, species) record — a VDJdb row, with no read to align — rather than from a sequenced fragment.

Markup and repair of bare (junction_aa, V, J) records — the VDJdb case.

Coordinate convention. Everything here is junction space: Cys104 through the Phe/Trp118 that opens FR4, both anchors included. That is what VDJdb’s cdr3 column actually holds (CASSARSGELFF with vEnd=4, jStart=7), and it is NOT arda’s cdr3 (which excludes both anchors). Conflating the two silently corrupts every coordinate emitted here.

The V and J germlines each template a known run of residues into the junction, and database/vdj/<organism>/cdr3_anchors.tsv ships them per allele. So marking up a record needs no germline search: align the junction’s 5’ end against the V’s templated residues (anchored at Cys104) and its 3’ end against the J’s (anchored at [FW]118), and read off the edit operations.

Both alignments are one semi-global Needleman-Wunsch anchored at the conserved residue with free end gaps on the junction-interior side. The free end gap is what makes the result honest: the germline templated run is an upper bound (V and J are exonuclease-trimmed), so the alignment stops wherever germline agreement stops paying for itself, and the untemplated N/D region is never scored. Concretely:

germline CASS    vs CCSS...   -> sub at index 1, d=1 -> repaired to CASS...
germline CASS    vs CGGS...   -> v_end = 1, no error (that is the V/N boundary)
germline TNEKLFF vs ...NEKLF  -> deletion, d=0 -> repaired to ...NEKLFF
germline TNEKLFF vs ...NNKLFF -> sub at index 8, d=4 -> REPORTED, not repaired

Detection and repair are deliberately separate: every germline disagreement is reported (with its position, extent and distance from the anchor), but only edits adjacent to the conserved anchor are applied. See _MAX_REPLACE for why.

Repair always targets a canonical junction. cdr3_repaired is only accepted when it opens with Cys104 and closes with Phe/Trp118 (_canonicalise); a repair that would hand back a junction missing either anchor is refused and the submission returned untouched. A repair exists to restore the anchors, and downstream every consumer trusts cdr3_repaired. So good implies canonical, by construction rather than by luck.

Fix-type names mirror VDJdb’s Cdr3Fixer so its cdr3fix JSON is directly comparable — on the committed 250-row fixture arda now reproduces VDJdb’s repair on all 100 records it flags, and agrees with its good/vCanonical/jCanonical verdicts on every row. The per-position errors list is arda’s addition.

class arda.cdr3fix.Cdr3Error(side, kind, pos, length, frm, to, dist=0, applied=False)

Bases: object

One edit between the observed junction and the germline-templated run.

pos indexes the observed junction and length is how far the error extends. frm is what the record has, to what the germline says. dist is the distance from the conserved anchor.

applied is true only when this edit was actually written into Cdr3Markup.cdr3_repaired. Being within _MAX_REPLACE of the anchor makes an edit eligible; the whole side’s repair is still discarded if its fix type comes back Failed* — no alignment, more than _MAX_FIX invented residues, more than _MAX_TRIM trimmed ones, or a result that would not be canonical. An error can therefore be reported with applied=False and the junction left alone — which is the point: detection and repair are separate decisions.

Parameters:

• side (str)

• kind (str)

• pos (int)

• length (int)

• frm (str)

• to (str)

• dist (int)

• applied (bool)

side: str

kind: str

pos: int

length: int

frm: str

to: str

dist: int = 0

applied: bool = False

class arda.cdr3fix.Cdr3Markup(cdr3, cdr3_repaired, v_call='', j_call='', locus='', species='', v_end=-1, j_start=-1, v_fix='FailedBadSegment', j_fix='FailedBadSegment', errors=<factory>, sequence_id='')

Bases: object

Result of marking up one (junction_aa, V, J) record.

Parameters:

• cdr3 (str)

• cdr3_repaired (str)

• v_call (str)

• j_call (str)

• locus (str)

• species (str)

• v_end (int)

• j_start (int)

• v_fix (str)

• j_fix (str)

• errors (list[Cdr3Error])

• sequence_id (str)

cdr3: str

cdr3_repaired: str

v_call: str = ''

j_call: str = ''

locus: str = ''

species: str = ''

v_end: int = -1

j_start: int = -1

v_fix: str = 'FailedBadSegment'

j_fix: str = 'FailedBadSegment'

errors: list[Cdr3Error]

sequence_id: str = ''

property v_canonical: bool

Does the junction as repaired open with the conserved Cys104?

Read off cdr3_repaired, not the submission – restoring the anchor is the whole point of the repair, and VDJdb’s vCanonical/jCanonical mean the same thing. Reading the submission instead disagreed with VDJdb on 76 of 250 fixture rows.

property j_canonical: bool

Does the junction as repaired close with the conserved Phe/Trp118?

property good: bool

Both sides repaired, and the result carries both conserved anchors.

A repair exists to produce a canonical junction. One that ends up without its Cys104 or its Phe/Trp118 has not repaired the record, it has invented a junction nobody submitted – so it can never be good (see _canonicalise).

property fix_needed: bool

to_cdr3fix()

The VDJdb cdr3fix JSON object, key-for-key.

Return type:

dict

explain()

One human-readable line: what happened, why, and where.

Return type:

str

class arda.cdr3fix.Anchor(locus, segment, templated_aa, functionality, status, source, anchor_nt=-1, partial_nt=0, germline_nt='')

Bases: object

A germline segment’s contribution to the junction.

Parameters:

• locus (str)

• segment (str)

• templated_aa (str)

• functionality (str)

• status (str)

• source (str)

• anchor_nt (int)

• partial_nt (int)

• germline_nt (str)

locus: str

segment: str

templated_aa: str

functionality: str

status: str

source: str

anchor_nt: int = -1

partial_nt: int = 0

germline_nt: str = ''

arda.cdr3fix.load_anchors(organism)

{(segment, allele): Anchor} for one organism; {} if not built.

Parameters:

organism (str)

Return type:

dict[tuple[str, str], Anchor]

arda.cdr3fix.markup_cdr3(cdr3, v_call, j_call, species='human', *, anchors=None, sequence_id='', max_replace=1)

Mark up and repair one junction. cdr3 is junction space (C..[FW]).

max_replace is how far from the conserved anchor an edit may sit and still be repaired; edits beyond it are reported with applied=False. Raising it repairs more, at the cost of rewriting N-region residues that merely look like germline (see _MAX_REPLACE).

Parameters:

• cdr3 (str)

• v_call (str)

• j_call (str)

• species (str)

• anchors (dict | None)

• sequence_id (str)

• max_replace (int)

Return type:

Cdr3Markup

arda.cdr3fix.markup_records(df, *, cdr3='cdr3', v='v', j='j', species='species', sequence_id=None, organism=None, max_replace=1)

Mark up a whole table. Anchors are loaded (and cached) once per organism.

Parameters:

• df (DataFrame)

• cdr3 (str)

• v (str)

• j (str)

• species (str)

• sequence_id (str | None)

• organism (str | None)

• max_replace (int)

Return type:

list[Cdr3Markup]

arda.cdr3fix.markup_batch(df, **kw)

markup_records + to_frame.

Parameters:

df (DataFrame)

Return type:

DataFrame

arda.cdr3fix.to_frame(records)

Records -> a TSV-ready frame with the vdjdb-compatible cdr3fix column.

Parameters:

records (Iterable[Cdr3Markup])

Return type:

DataFrame

arda.cdr3fix.format_report(records, *, show_ok=False)

Human-readable log: a summary table, then a line per fixed/failed record.

show_ok=True also lists the records that needed no repair.

Parameters:

• records (Iterable[Cdr3Markup])

• show_ok (bool)

Return type:

str

Which D, and where — from an amino-acid junction alone.

A VDJdb-style record has no nucleotides, and D segments are short and trimmed at both ends, so the D is often invisible in the translated junction. Two independent sources of information remain, and they are complementary:

Where. The junction’s nucleotide length is known (3x its amino-acid length), so insVD + |D surviving| + insDJ is pinned. Marginalising the generative model’s insertion-length and D-trimming distributions therefore places the D even when the sequence says nothing at all about it. Measured against OLGA ground truth, the MAP d_start is a median 1 nt off for mouse TRB, 2 nt for human TRB, and 3 nt for TRD and IGH.

Which. The length constraint is nearly useless for identity — the D length distributions overlap, so the posterior barely moves off the prior. What identity the prior does carry is P(D | J), and for TRB that is mostly genomic order: TRBD2 lies 3’ of the whole TRBJ1 cluster, so a TRBJ1 junction can only have used TRBD1 (see _mask_forbidden). What otherwise identifies a D is the amino-acid match, and only where enough D survives: median surviving D is 17 nt for IGH (~5.7 aa) but 5 nt for human TRB (~1.7 aa).

So neither source alone is enough, and which one dominates flips by locus:

locus prior only aa only combined n (held-out seed, generated) human IGH 15 % 81 % 82 % 345 human TRB 76 % 70 % 82 % 595 human TRD 86 % 88 % 87 % 699 mouse TRB 76 % 83 % 85 % 699

“prior only” is beta = 0; “aa only” is argmax of the match score under a uniform prior, ties broken by marginal usage. Combining wins at IGH and both TRB; TRD is a wash, because one D gene (TRDD3) accounts for 85 % of rearrangements and the aa match already finds it.

The combination is log P(D | M, J) + beta * s_D: the length-and-J prior, tempered by the best gapless local alignment score s_D of the D’s three-frame translations against the non-templated middle of the junction. beta is fitted per locus and shipped in database/vdj/<org>/d_prior.tsv with the distributions themselves, so nothing here needs OLGA at runtime. It is flat above ~1.25 for TRB, so the shipped values are not delicate.

Honesty about the numbers. The table is measured on junctions drawn from the same generative model that supplies the prior, so the prior’s contribution is flattered. The amino-acid contribution is not — it is germline matching. Rearrangements that genomic order forbids are excluded from the truth: OLGA’s human TRB model emits TRBD2 x TRBJ1 in 8.7 % of draws, and scoring against those measures agreement with a model artifact.

Out of model, against nucleotide D calls (E <= 0.05) on the real GenBank fixtures: human TRB 94 %, IGH 85 %, TRD 91 %, mouse TRB 85 %. On TRB, note that both this posterior and the nucleotide caller enforce the same D2-x-J1 constraint, so their agreement on TRBJ1 records is guaranteed rather than earned; the TRBJ2 rows, where both D genes stay possible, score 91 % (human) and 81 % (mouse).

Priors exist only for the (organism, locus) pairs with a published model: human IGH, TRB and TRD, and mouse TRB. Everything else returns None rather than guessing.

class arda.dpost.DPosterior(locus, d_call, posterior, entropy, by_gene=<factory>, support_aa=0, d_start=-1, d_start_ci90=(-1, -1), n_middle_nt=0)

Bases: object

Posterior over the D gene, and over where it sits in the junction.

Parameters:

• locus (str)

• d_call (str)

• posterior (float)

• entropy (float)

• by_gene (dict[str, float])

• support_aa (int)

• d_start (int)

• d_start_ci90 (tuple[int, int])

• n_middle_nt (int)

locus: str

d_call: str

posterior: float

entropy: float

by_gene: dict[str, float]

support_aa: int = 0

d_start: int = -1

d_start_ci90: tuple[int, int] = (-1, -1)

n_middle_nt: int = 0

property confident: bool

A hard call. 0.9 keeps ~the top decile of TRB records and most of IGH.

arda.dpost.posterior_d(junction_aa, v_call, j_call, species='human')

Posterior over the D gene (and its position) for an amino-acid junction.

junction_aa is junction space (Cys104 .. Phe/Trp118, both included), as in arda.cdr3fix. Returns None when the locus has no D, no shipped model, or the junction cannot be marked up.

Parameters:

• junction_aa (str)

• v_call (str)

• j_call (str)

• species (str)

Return type:

DPosterior | None

arda.dpost.load_d_prior(organism)

{locus: DPrior}; empty when the organism has no shipped model.

Parameters:

organism (str)

Return type:

dict[str, DPrior]

Bulk RNA-seq

Stage 1 — RNA-seq filter + map.

Reuses the streaming, memory-bounded annotator (annotate.mapper._prep + _annotate_chunk + the background-reader/bounded-queue loop of annotate_file): MMseqs2 is the parallel layer and its k-mer prefilter rejects non-receptor reads before alignment, so mostly-non-receptor RNA-seq is cheap. The difference from arda annotate is that we write only the reads that map (keyed by read id, so the AIRR TSV is the read-id → junction map), plus an optional candidate FASTA and a run report.

Paired FASTQ mates are streamed independently, tagged <id>/1 / <id>/2 so query ids stay unique; a pair is kept if either mate maps (recall-first — the base id recovers the pair).

arda.rnaseq.map.map_rnaseq(r1, output, *, r2=None, organism='human', seqtype='nt', threads=0, sensitivity=None, strand='both', chunk_size=400000, map_d=True, d_max_evalue=None, reconstruct=False, min_score=75.0, max_seqs=300, kmer=-1, drop_constant_only=True, limit=None, emit_reads=None, report_path=None, two_pass=False, adaptive=False, fast_segments=False, indel_rescue=False, segment_only_v=False, prefilter=False, with_junction_quality=False, with_mutation_quality=False, shm='framework', complete_junction_nt=0)

Filter + map an RNA-seq FASTQ (single or paired); write mapped reads as AIRR.

Parameters:

• r1 (str | Path) – FASTA/FASTQ (gzip by .gz). Single-end, or R1 of a pair.

• output (str | Path) – AIRR TSV of the mapped reads only (keyed by sequence_id).

• r2 (str | Path | None) – R2 FASTQ for paired input; None for single-end.

• min_score (float) – drop mapped reads below this MMseqs2 bit score. 0 disables the filter (recall-max). See _MIN_SCORE for the calibration.

• kmer (int | None) – MMseqs2 -k. The memory knob: the nucleotide prefilter allocates 4**k index entries, so the tool default k=15 costs ~8.4 GB peak RSS whatever else you set. arda defaults to 13 (~0.7 GB, and never slower). None = MMseqs2’s default.

• max_seqs (int) – MMseqs2 target hits per read. Does not change which reads are kept, only which V/J scaffold wins. See arda.annotate.mapper._MAX_SEQS.

• limit (int | None) – analyse only the first limit reads (single-end) / read pairs (paired), then stop — a native head, so a subsample no longer needs an external zcat | head | gzip round-trip. None maps the whole file.

• with_junction_quality (bool) – also emit a junction_quality column — the read’s Phred+33 string over exactly the bases of junction, same orientation (see junction_quality()). OFF by default: it appends a non-schema column, so the default output stays byte-identical. This is the only place the FASTQ quality is still in hand — Stage 1 otherwise discards it — and it is what correct --min-junction-q gates on. Refused with reconstruct (a merged fragment has no single input quality string).

• emit_reads (str | Path | None) – optional path — write the mapped reads’ sequences as FASTA (coding-strand oriented) for downstream handoff.

• report_path (str | Path | None) – optional path — write the RnaseqReport as JSON.

• fast_segments (bool) – with two_pass, answer the segment pass structurally instead of with mmseqs search – 37x on that step, agreeing with it on .9997 of V alleles and .9998 of J. It only NOMINATES candidates; the winner is still aligned against the full scaffold by MMseqs2, so the contract is that the AIRR output does not move. Ignored without two_pass, since there is no segment pass to replace.

• indel_rescue (bool) – with fast_segments, send reads carrying the two-diagonal signature of an indel to the GAPPED rescue path instead of resolving them on the fast path. One ungapped extension scores such a read only up to the indel, so its segment score is systematically low. Measured on 341,294 real IGH mates: 3.18 % of reads carry a V indel, rising to 8.00 % below 90 % V identity. Reroutes, never drops.

• adaptive (bool) – cap alignments per read and re-search only the reads whose capped score is low (arda.annotate.mapper._extend_uncertain()). Measured 2.17x on 1 M bulk reads with zero reads lost — but read preservation is not the whole guarantee. OFF by default: on the real-read fixture it also changes junction_aa on 3 of 453 reads, and two of them scored 128 and 131, far above the 90-bit trigger. So a high score does NOT certify that the best alignment was found, and the trigger cannot be calibrated on score alone. Opt in only where a junction-level difference is acceptable.

• two_pass (bool) –

  use the segment reference to shortlist a single V×J scaffold per read before aligning (arda.annotate.mapper._segment_best_hits()). Reads it cannot resolve are realigned against the full reference, so nothing is dropped — see arda.annotate.shortlist. Requires segments.fasta (written by arda build-index); silently falls back to the one-pass search when it is absent.

  ⛔ The win is set by whether reads SPAN V INTO J, not by the library type, and the predictor is fast_fraction in the report. Measured: 3.51× on a TCR amplicon (fast path 85 %), 2.96× on a 100 %-receptor human TRB set (95.6 %), 2.64× on mouse TRA (89 %) — but 1.03× slower on the human IGH leg of that same 100 %-receptor dataset (16.3 %: those reads cover V and stop short of the short IGHJ target), and 0.762× on 2.74 %-receptor bulk (5 %). Off by default because no library type predicts it; run a sample and read fast_fraction.

• organism (str)

• seqtype (str)

• threads (int)

• sensitivity (float | None)

• strand (str)

• chunk_size (int)

• map_d (bool)

• d_max_evalue (float | None)

• reconstruct (bool)

• drop_constant_only (bool)

• segment_only_v (bool)

• prefilter (bool)

• with_mutation_quality (bool)

• shm (str)

• complete_junction_nt (int)

Returns:

The run RnaseqReport (also printed by the CLI).

Return type:

RnaseqReport

arda.rnaseq.map.read_pairs(r1, r2=None, *, reconstruct=False, limit=None, with_qual=False)

Stream (id, sequence) reads for single-end (r1 only) or paired input.

with_qual yields (id, sequence, quality) instead — the Phred+33 string, "" for FASTA input. It is incompatible with reconstruct: a merged fragment’s bases come from two different reads, so no single input quality string describes it (map_rnaseq() refuses the combination rather than emitting a quality that does not belong to the sequence).

For paired input the two mates carry the same id, so they are tagged <id>/1 and <id>/2 to keep query ids unique (strip the suffix to recover the pair). With reconstruct, overlapping mates are merged into one fragment (merge_pair()) keyed by the bare id — giving a short read the mate’s V/J context; non-overlapping mates fall back to the tagged-independent form.

Parameters:

• limit (int | None) – analyse only the first limit input records — reads (single-end) or read pairs (paired) — then stop, without decompressing the rest of the file. None reads everything. The mate-order / truncation checks below still run on every record actually read; a truncation beyond limit is simply never reached — that is the intent of a head-style limit, not a hole in the check.

• r1 (str | Path)

• r2 (str | Path | None)

• reconstruct (bool)

• with_qual (bool)

Raises:

ValueError – if the two files disagree on read names or record count. This is not paranoia: a truncated R2 makes zip stop early and silently analyse a prefix, and a shuffled R2 pairs mate 1 of one fragment with mate 2 of another. Both were observed in this project’s own data and produced a published false discovery (a spurious R2-only blind spot) that had to be retracted. A pair of FASTQs is an assertion; check it.

Return type:

Iterator[tuple]

arda.rnaseq.map.merge_pair(s1, s2, *, q1=None, q2=None, min_overlap=12, max_mismatch_rate=0.1)

Overlap-merge a read pair into one fragment, or None if they don’t overlap.

Aligns s1 (R1) with reverse_complement(s2) (R2 flipped to the same strand) by finding an exact _MERGE_ANCHOR-mer from the flipped R2’s 5’ end inside R1 (C-level str.find → O(len), so non-overlapping pairs — the common RNA-seq case — cost almost nothing), then verifying the implied overlap; the mate provides the V/J context a short read lacks.

In the overlap the two mates may disagree. Given Phred qualities q1/q2 the base with the higher quality wins per position (rc(R2)’s quality is q2 reversed, not complemented); without them R2 wins the whole overlap (the historical behaviour). Outside the overlap, R1 supplies its 5’ part and R2 its 3’ tail.

Parameters:

• s1 (str)

• s2 (str)

• q1 (str | None)

• q2 (str | None)

• min_overlap (int)

• max_mismatch_rate (float)

Return type:

str | None

arda.rnaseq.map.junction_quality(rec, qual)

The read’s Phred+33 substring covering rec["junction"], or "".

⛔ The quality string belongs to the read AS SUBMITTED, while the junction and every coordinate in the record are on the CODING strand. For rev_comp == "T" the two run in opposite directions, so the quality is reversed (not complemented – a Phred char has no complement) before slicing. Getting that backwards yields a string of the right LENGTH holding the wrong bases’ qualities, which no length or format check downstream can catch. So the slice is verified against the junction it claims to describe before it is returned.

The junction is coding[cdr3_start - 3 : cdr3_end + 3] – the CDR3 flanked by the Cys104 and [FW]118 anchor codons, see arda.annotate.transfer._junction_nt() – so the coordinates place it in O(1). They can be absent (a producer that did not fill the region columns), hence the find fallback; if neither reproduces the junction, return "" rather than a misaligned string.

Parameters:

• rec (dict)

• qual (str)

Return type:

str

arda.rnaseq.map.mutation_quality(rec, qual)

Per-mutation Phred scores for rec["v_mutations"] / rec["j_mutations"].

⛔ Driven by the mutation list that was EMITTED, not by re-deriving one. Walking the alignment and scoring every mismatch reproduces what _markup.segment_cigars found – which since 2.16.0 is a SUPERSET of what the columns carry, because arda.shm then drops the junction-internal entries (measured on this repo’s own fixture: 25 of 242 V rows had more mismatches than mutations). The result lines up in length only by accident and pairs entry i with a different base’s score. So the walk builds germline position -> query position and each emitted entry looks its own position up; an entry whose position the alignment does not cover yields "" for the whole segment rather than a short, misaligned list.

⛔ Quality is oriented like the READ AS SUBMITTED; the alignment and every coordinate here are on the coding strand. Same reversal rule as junction_quality().

ponytail: a Python pass over the alignment, not a second output from the C++ walk that already visits these columns. It rides its own flag, so it never lands on a mode run; move it into segment_cigars if it ever shows up in a profile.

Parameters:

• rec (dict)

• qual (str)

Return type:

dict[str, str]

class arda.rnaseq.map.RnaseqReport(input, organism, paired=False, input_bytes=0, read_length_min=0, read_length_max=0, read_length_mean=0.0, total_reads=0, mapped_reads=0, per_locus=<factory>, constant_only_fragments=0, isotype_from_mate=0, min_score=0.0, threads=0, wall_seconds=0.0, reads_per_second=0.0, peak_rss_mb=0.0, rss_gain_mb=0.0, segment_search=<factory>, prefilter_stats=<factory>)

Bases: object

Counts + timing for one map run (written as JSON with --report).

Parameters:

• input (str)

• organism (str)

• paired (bool)

• input_bytes (int)

• read_length_min (int)

• read_length_max (int)

• read_length_mean (float)

• total_reads (int)

• mapped_reads (int)

• per_locus (dict[str, int])

• constant_only_fragments (int)

• isotype_from_mate (int)

• min_score (float)

• threads (int)

• wall_seconds (float)

• reads_per_second (float)

• peak_rss_mb (float)

• rss_gain_mb (float)

• segment_search (dict)

• prefilter_stats (dict)

input: str

organism: str

paired: bool = False

the AIRR carries only the reads that mapped, so its row count, its sequence lengths and its mate suffixes all describe the receptor subset rather than the library. arda stats reports them.

Type:

Library shape, recorded here because nothing downstream can recover it

input_bytes: int = 0

read_length_min: int = 0

read_length_max: int = 0

read_length_mean: float = 0.0

total_reads: int = 0

mapped_reads: int = 0

per_locus: dict[str, int]

constant_only_fragments: int = 0

isotype_from_mate: int = 0

min_score: float = 0.0

threads: int = 0

wall_seconds: float = 0.0

reads_per_second: float = 0.0

peak_rss_mb: float = 0.0

rss_gain_mb: float = 0.0

segment_search: dict

prefilter_stats: dict

property mapped_fraction: float

as_dict()

Return type:

dict

Stage 3 — contig assembly (anchored greedy overlap-extension).

Reconstruct the clonotypes that Stage 1 maps but cannot call: a long CDR3 (V(DD)J ultralong, ~20-40 aa) does not fit in one 100-150 bp read, so no read spans the junction and correct’s complete-junction filter drops every read of it. Assembly-based extractors recover these by assembly; this module does the same on the reads Stage 1 already mapped.

Why overlap-extension and not a de Bruijn graph: every clonotype sharing a germline V/J contributes identical k-mers, so a dBG collapses distinct clones exactly across the CDR3 (the region of interest). This is the reason the Pevzner-lab Ig assemblers use a read graph, not a k-mer graph (Safonova 2015, 10.1093/bioinformatics/btv238). We exploit arda’s own anchors instead: Stage 1 gives every V-side read a cdr3_start offset, so reads of one clone are already coordinate-aligned at the CDR3 – we seed from those and extend 3’ through the CDR3 into J, where the sequence is clone-specific, and stop before running deep into the (shared) constant region. Seeding never extends 5’ into the germline V, which is what keeps distinct clones apart and bounds the germline-k-mer blow-up.

The assembled contig physically unites a clone’s V-side reads with its J/C-side reads under one junction, so correct also gets the clone’s isotype (from the J/C mates’ c_class) for free – the long clones were previously invisible to both.

Output is a per-member-read AIRR fragment (sequence_id -> the contig’s complete junction), meant to be concatenated with the Stage-1 mapped AIRR and fed to correct in one pass: a read that was incomplete in Stage 1 is dropped there and kept here, so fragments count once.

arda.rnaseq.assemble.assemble_contigs(airr_tsv, output, *, organism='human', k=21, min_overlap=21, min_id=0.9, max_ext_past_cdr3=130, scan_cap=400, threads=0, map_d=True, d_max_evalue=None, report_path=None)

Assemble long-CDR3 contigs from Stage-1 mapped reads and attribute their junctions.

Reads the Stage-1 mapped-reads AIRR (needs sequence, rev_comp, locus, cdr3_start and the v/j/c_call columns), assembles per-locus contigs, re-annotates them (reannotate_contigs()), and writes an AIRR TSV with one row per incomplete member read carrying the contig’s complete junction (and the read’s own c_class so isotype survives). Concatenate this with the mapped AIRR and run correct once: the read’s incomplete Stage-1 row is dropped and this complete row kept, so each fragment is counted exactly once.

The contig’s D call travels with the junction (d_call, d2_call, d_support, d2_support, np1-np3). An ultralong CDR3 is the one place a tandem D-D is both most likely and least visible to a single read, so the contig is where it must be called.

Parameters:

• airr_tsv (str | Path) – Stage-1 mapped-reads AIRR TSV.

• output (str | Path) – assembled-reads AIRR TSV (header only if nothing assembles).

• map_d (bool) – map D segments on the assembled contig (default True).

• d_max_evalue (float | None) – E-value gate on the D call(s); None keeps the shipped 0.2.

• max_ext_past_cdr3 (int) – stop extending a contig once it reaches this many nt past the CDR3 start – enough to cross the junction into J without running into the shared C region.

• scan_cap (int) – per-step cap on candidate reads examined for a (germline-frequent) k-mer.

• organism (str)

• k (int)

• min_overlap (int)

• min_id (float)

• threads (int)

• report_path (str | Path | None)

Returns:

An AssembleReport.

Return type:

AssembleReport

class arda.rnaseq.assemble.AssembleReport(reads_in: 'int' = 0, seeds: 'int' = 0, contigs: 'int' = 0, contigs_complete: 'int' = 0, reads_rescued: 'int' = 0, members_without_junction: 'int' = 0, wall_seconds: 'float' = 0.0, peak_rss_mb: 'float' = 0.0, rss_gain_mb: 'float' = 0.0)

Bases: object

Parameters:

• reads_in (int)

• seeds (int)

• contigs (int)

• contigs_complete (int)

• reads_rescued (int)

• members_without_junction (int)

• wall_seconds (float)

• peak_rss_mb (float)

• rss_gain_mb (float)

reads_in: int = 0

seeds: int = 0

contigs: int = 0

contigs_complete: int = 0

reads_rescued: int = 0

members_without_junction: int = 0

Members skipped because their span did not cover the contig’s junction by _MIN_JUNCTION_COVER – recruited on germline-only overlap, so there is no evidence they belong to this clonotype. They keep their reads; only the attribution is withheld.

wall_seconds: float = 0.0

peak_rss_mb: float = 0.0

rss_gain_mb: float = 0.0

as_dict()

Return type:

dict

Stage 2 — CDR3 error correction (sequencing-error model).

Collapses sequencing-error CDR3 variants onto their parent clonotype, using seqtree neighbour search (a fast edit-bounded index) to find substitution/indel neighbours. A clonotype C is an error child of a more-abundant neighbour P (differing by n_subs substitutions and n_indel inserted/deleted bases) iff the expected number of such misread parent reads – count[P] * p_sub**n_subs * p_ind**n_indel – is at least count[C]. The rates are PER BASE and the per-mismatch probability is length-scaled (p_sub = error_rate * L): a single mismatch over a longer junction sheds proportionally more error mass, so the default error_rate = 0.001 reproduces vdjtools’ ~1/20 at a 45 nt (15 aa) junction and scales elsewhere. A multi-base (in-frame SHM) indel costs p_ind**len and is kept as a real clonotype. The count is the SPANNING read depth – reads that fully observe the junction – so the test is over the reads that actually saw the discriminating base ("2/2, not 2/200"); error_method in {binom, betabinom} instead piles up partial reads per position for extra depth at very low coverage. Children route to the parent; chains collapse to the ultimate ancestor; count[parent] * p_err >= count[child] with p_err < 1 gives strictly increasing counts along parent pointers, so there are no cycles.

arda.rnaseq.correct.correct_airr(airr_tsv, output, *, organism='human', map_d=True, d_max_evalue=None, max_subs=3, max_indel=0, error_rate=0.001, indel_rate=0.001, require_vj=True, error_method=None, ec_mode='fast', clonotype_key='full', call_level='allele', flag_chimeras=False, isotype=True, min_junction_q=None, complete_only=True, coverage=True, read_map=None, extra_airr=None, report_path=None)

Aggregate mapped reads into clonotypes and collapse CDR3 sequencing errors.

Parameters:

• airr_tsv (str | Path) – Stage-1 mapped-reads AIRR TSV (needs junction, sequence_id).

• organism (str) – reference organism, used only to map D into each clonotype’s junction.

• map_d (bool) – append the D columns (_D_STR_COLS + _D_INT_COLS), called once per clonotype on its corrected junction (see _clonotype_d()). Default True.

• d_max_evalue (float | None) – E-value gate on the D call(s); None keeps the shipped 0.2. Lower is stricter – 0.01 is the band where D agrees .9985 with IgBLAST on a TRB amplicon.

• output (str | Path) – corrected clonotype table TSV (junction, junction_aa, v_call, j_call, c_call, locus, duplicate_count, consensus_count, and with map_d the D columns), sorted by abundance. A clonotype is keyed by (locus, v_call, j_call, junction). Per the AIRR schema, duplicate_count is the number of READS supporting the clonotype (both paired mates of a molecule count) and consensus_count is the number of distinct fragment consensuses (the two mates of one molecule are one consensus). c_call is the clonotype’s dominant isotype CLASS (from c_class: IGHG, IGHA, …), preferring a resolved class over the ambiguous IGHC; empty when no read carried a constant call.

• max_subs (int) – max substitutions between an error child and its parent (seqtree neighbour search). This is a SEARCH RADIUS, not a threshold – the accept/reject decision is the length-scaled probability model above, so widening it only lets the model SEE parents it would already have accepted. The default was 2 through 2.9.0, which truncated the search below what the model would take on a deep clone: on the two monoclonal cell lines in the arda-benchmark set 2 -> 3 collapses Jurkat 74 -> 57 clonotypes and Raji 91 -> 58, while a polyclonal mouse spleen (7,942) and an oligoclonal B-LCL (13) are UNCHANGED at 2, 3 and 4 – the model refuses those collapses on abundance regardless of radius. It saturates at 3 (4 gives the same four numbers), so 3 is the default.

• max_indel (int) – max inserted/deleted bases searched for indel error children (default 0). A 1-2 bp instrument indel is a frameshift and is already dropped by complete_only, so on complete junctions the indel search only costs time (~160x slower) and collapses nothing; a multi-base in-frame SHM indel costs (indel_rate*L)**len and is kept as a real clonotype either way. Set it > 0 only with --all-junctions (frameshift indels kept).

• error_rate (float) – per-BASE substitution error rate (~Phred 30 = 0.001). The per-substitution collapse probability is length-scaled, error_rate * junction_len, so the default reproduces vdjtools’ ~1/20 at a 45 nt (15 aa) junction and scales for other lengths.

• indel_rate (float) – per-BASE indel error rate (instrument-dependent; default 0.001, length-scaled).

• ec_mode (str) – knob preset, "fast" (default, = today’s shipped behaviour) or "accurate". See EC_MODES; an explicitly passed error_method / min_junction_q overrides it.

• min_junction_q (int | None) – reassign onto its parent a read whose junction differs from that parent below this Phred score (_quality_gate()). 0 disables it. Needs the junction_quality column from map --junction-quality and RAISES without it.

• error_method (str | None) – "simple" (default) tests on spanning read counts; "binom" / "betabinom" pile up partial reads per discriminating position for extra depth at very low coverage (_error_pileup()).

• require_vj (bool) – only collapse neighbours sharing v_call and j_call (default True – a true sequencing error does not change the germline-anchored V/J call).

• complete_only (bool) – keep only reads whose junction spans both conserved anchors, is in frame, and has no stop codon (see _COMPLETE). A read that stops short of the [FW]118 anchor yields a prefix of a junction, not a clonotype. Setting this False reproduces the raw per-read behaviour and is almost never what you want. (This governs which reads DEFINE clonotypes, not how they are counted – see coverage.)

• coverage (bool) – count a clonotype’s abundance as EVERY read that encompasses its junction (aligns to it), not only the reads that span it end-to-end (default True). A long CDR3 is covered by many partial V-side / J-side reads that never reach both anchors; counting only spanning reads under-reports it non-uniformly (the deficit scales with CDR3 length). Coverage counting (_assign_coverage()) is the true expression estimate. False reverts to spanning-read counts.

• read_map (str | Path | None) – optional TSV sequence_id -> junction (the corrected clonotype a read ends up in) — the read-id → junction map after correction.

• extra_airr (str | Path | None) – optional Stage-3 assembled-reads AIRR (from assemble_contigs()), concatenated with airr_tsv before aggregation. Its rows carry a contig’s complete junction for reads whose own Stage-1 junction was incomplete, so a long-CDR3 clone no single read spans is counted once (the read’s incomplete Stage-1 row is dropped by complete_only).

• clonotype_key (str)

• call_level (str)

• flag_chimeras (bool)

• isotype (bool)

• report_path (str | Path | None)

Returns:

A CorrectReport.

Return type:

CorrectReport

class arda.rnaseq.correct.CorrectReport(clonotypes_in: 'int' = 0, clonotypes_out: 'int' = 0, reads: 'int' = 0, reads_assigned: 'int' = 0, collapsed: 'int' = 0, reads_with_junction: 'int' = 0, empty_clonotypes: 'int' = 0, reads_from_assembly: 'int' = 0, reads_incomplete: 'int' = 0, reads_low_quality: 'int' = 0, clonotypes_low_quality: 'int' = 0, rescued_clonotypes: 'int' = 0, rescued_reads: 'int' = 0, orphan_clonotypes: 'int' = 0, orphan_reads: 'int' = 0, wall_seconds: 'float' = 0.0, peak_rss_mb: 'float' = 0.0, rss_gain_mb: 'float' = 0.0)

Bases: object

Parameters:

• clonotypes_in (int)

• clonotypes_out (int)

• reads (int)

• reads_assigned (int)

• collapsed (int)

• reads_with_junction (int)

• empty_clonotypes (int)

• reads_from_assembly (int)

• reads_incomplete (int)

• reads_low_quality (int)

• clonotypes_low_quality (int)

• rescued_clonotypes (int)

• rescued_reads (int)

• orphan_clonotypes (int)

• orphan_reads (int)

• wall_seconds (float)

• peak_rss_mb (float)

• rss_gain_mb (float)

clonotypes_in: int = 0

clonotypes_out: int = 0

reads: int = 0

⛔ SPANNING reads entering Stage 2, counted BEFORE any correction runs. It is therefore invariant to everything correct does, and it is NOT the read-conservation quantity. Comparing it across --ec-mode shows 0 on every sample and reads exactly like the invariant holding – which is how a 1.39 % leak on a full-depth TRA amplicon was missed. The invariant is reads_assigned.

reads_assigned: int = 0

sum(duplicate_count) over the emitted clonotype table, i.e. every read the correction actually placed. This is what must not fall when a denoising mode is switched on – error correction MOVES reads onto a parent, it never discards them. Reported so a run is self-checking instead of relying on a docstring.

Type:

The read-conservation invariant

collapsed: int = 0

reads_with_junction: int = 0

empty_clonotypes: int = 0

Clonotypes dropped because no read was assigned to them – see the note at the emit step. Only reachable with –all-junctions, where a rescued read’s truncated Stage-1 junction is itself a clonotype key and loses its reads to the assembled row.

reads_from_assembly: int = 0

Junction-bearing rows contributed by –assemble (extra_airr). A rescued read appears in BOTH frames – its incomplete Stage-1 row and its assembled row – so these are reported separately rather than folded into reads_with_junction, which is a Stage-1 statistic.

reads_incomplete: int = 0

reads_low_quality: int = 0

clonotypes_low_quality: int = 0

rescued_clonotypes: int = 0

rescued_reads: int = 0

orphan_clonotypes: int = 0

orphan_reads: int = 0

wall_seconds: float = 0.0

peak_rss_mb: float = 0.0

rss_gain_mb: float = 0.0

as_dict()

Return type:

dict

arda.rnaseq.correct.CLONOTYPE_KEYS = ('full', 'junction')

How a clonotype is identified.

full – (locus, v_call, j_call, junction), the historical key. junction – (locus, junction): the V/J calls are canonicalised to the junction’s majority before grouping, so call splits collapse. That class is invisible to every error model – a junction byte-identical to an abundant clone’s under a different V or J call has no discriminating base to score – and on Jurkat it is the largest error class BY READS (130 of 14,531), including an allele-level TRG split. Measured cost on a POLYCLONAL TRA amplicon: 132 of 19,956 clonotypes merge (0.66 %), and in every ambiguous case inspected the minority call carried ONE read against 4-10 for the majority on a short 30-39 nt junction – a call error on a low-abundance read, not a second clone. Benefit on Jurkat: TRB 35 -> 33 clonotypes at purity .99096 -> .99696, reads unchanged.

arda.rnaseq.correct.CALL_LEVELS = ('allele', 'gene')

At what resolution a V/J call names a germline.

allele – TRGJ1*01, the historical behaviour and the default. gene – TRGJ1: the allele suffix is dropped before the clonotype key is formed, so allele-level call splits collapse. Jurkat’s largest single split is exactly that shape (TRGJ1*01 64 reads against TRGJ1*02 140 on the same junction), and no error model can see it because two identical junctions have no discriminating base.

⚠ It also collapses a TIE LIST whose members differ only by allele (TRAV1*01,TRAV1*02 -> TRAV1), which is the artifact behind the 14-point spread between v_allele_exact (median .8328) and v_allele_resolved (.9763) across 25 cluster datasets. ⚠ IGH carries 4.33 alleles/gene, ~2x every other locus, so expect the effect to be far larger there than on TR — measure per library before quoting a cost.

Stages 2-3 as one function, shared by the single-node and sharded delivery paths.

The RNA-seq pipeline is map -> assemble -> correct. Stage 1 (map) is per-read and shards perfectly; Stages 2 and 3 are global and do not shard at all:

• correct collapses sequencing-error variants onto a parent clonotype and counts distinct fragments. Run per shard, one clone split across N shards is counted N times, and error children collapse against a fraction of their parent’s depth.

• assemble grows contigs across reads. Reads that tile one long CDR3 in different shards never meet, so the contig is never built – which is precisely the class Stage 3 exists for.

So the sharded path runs map per shard, concatenates the Stage-1 AIRR in shard order, and then calls exactly the same finish() the single-node path calls. Not “the same steps” – the same function, so the two cannot drift apart in a parameter.

That plus contiguous shards (arda.cluster.split_pairs()) is what makes a sharded run byte-identical to a single-node one, rather than merely similar.

arda.rnaseq.pipeline.finish(airr, out_dir, out_prefix, *, organism='human', threads=0, assemble=True, complete_only=True, map_d=True, d_max_evalue=None, ec_mode='fast', min_junction_q=None, clonotype_key='full', call_level='allele', isotype=True, map_report=None, write_qc=True, echo=None)

Run Stages 2-3 over a Stage-1 AIRR and write the clonotype table + merged report.

Called by both run() (single node) and reduce() (after a sharded Stage 1).

Parameters:

• airr (str | Path) – Stage-1 mapped-reads AIRR TSV.

• out_dir (str | Path) – where the outputs land (see OUTPUTS).

• out_prefix (str) – where the outputs land (see OUTPUTS).

• map_report (dict | None) – the Stage-1 report to embed; for a sharded run, the merged one.

• echo – optional print-like callback for progress lines.

• organism (str)

• threads (int)

• assemble (bool)

• complete_only (bool)

• map_d (bool)

• d_max_evalue (float | None)

• ec_mode (str)

• min_junction_q (int | None)

• clonotype_key (str)

• call_level (str)

• isotype (bool)

• write_qc (bool)

Returns:

The merged report dict, as written to <prefix>.arda.json.

Return type:

dict

arda.rnaseq.pipeline.run(r1, out_dir, out_prefix, *, r2=None, organism='human', threads=0, reconstruct=False, min_score=75.0, kmer=12, assemble=True, complete_only=True, map_d=True, d_max_evalue=None, limit=None, two_pass=False, adaptive=False, fast_segments=False, prefilter=False, segment_only_v=False, indel_rescue=False, ec_mode='fast', min_junction_q=None, clonotype_key='full', call_level='allele', isotype=True, shm='framework', complete_junction_nt=0, echo=None)

Single-node map -> assemble -> correct.

Parameters:

• r1 (str | Path)

• out_dir (str | Path)

• out_prefix (str)

• r2 (str | Path | None)

• organism (str)

• threads (int)

• reconstruct (bool)

• min_score (float)

• kmer (int | None)

• assemble (bool)

• complete_only (bool)

• map_d (bool)

• d_max_evalue (float | None)

• limit (int | None)

• two_pass (bool)

• adaptive (bool)

• fast_segments (bool)

• prefilter (bool)

• segment_only_v (bool)

• indel_rescue (bool)

• ec_mode (str)

• min_junction_q (int | None)

• clonotype_key (str)

• call_level (str)

• isotype (bool)

• shm (str)

• complete_junction_nt (int)

Return type:

dict

arda.rnaseq.pipeline.reduce(shard_dir, out_dir, out_prefix, *, organism='human', threads=0, assemble=True, complete_only=True, map_d=True, d_max_evalue=None, echo=None)

Merge sharded Stage-1 output, then run Stages 2-3 once over the whole thing.

The shard AIRRs are merged from an explicit sorted list, not a bare *.tsv glob: shard names are zero-padded so sorted() is numeric (shard_10 must not precede shard_2 – concatenation order is read order here), and naming the glob means reduce can never swallow its own clones.tsv if someone points --out-dir at --shard-dir.

Parameters:

• shard_dir (str | Path)

• out_dir (str | Path)

• out_prefix (str)

• organism (str)

• threads (int)

• assemble (bool)

• complete_only (bool)

• map_d (bool)

• d_max_evalue (float | None)

Return type:

dict

arda.rnaseq.pipeline.write_stats_for(out_dir, out_prefix, *, organism='human', say=None)

Write <prefix>.stats.tsv from the run’s own artifacts. Returns the row count.

⛔ Written unconditionally, not behind a flag. It reads only files that already exist and costs one pass over each; the alternative is that the numbers an operator needs to decide whether a sample is usable exist only if they knew to ask for them BEFORE the run.

⛔ Called AFTER the report JSON is final. run() rewrites it with the whole-run wall time once Stages 2-3 return, so collecting inside finish() would put the Stage-2/3 time in the run scope under the name wall_seconds – a wrong number that looks like a right one.

Parameters:

• out_dir (str | Path)

• out_prefix (str)

• organism (str)

Return type:

int

arda.rnaseq.pipeline.OUTPUTS = {'airr': '{prefix}.airr.tsv', 'assembled_airr': '{prefix}.assembled.airr.tsv', 'clones': '{prefix}.clones.tsv', 'report': '{prefix}.arda.json', 'stats': '{prefix}.stats.tsv'}

Output basenames, relative to out_dir and given prefix.

Wall time and peak RSS for one pipeline stage.

Only map used to report resources, which left the expensive stage unmeasured: mapping is flat (~300-650 MB at any depth) while Stage 3 holds the clone set, and on a B-cell-rich tumour (28,444 clonotypes from 105M reads) correct peaked at 2,071.7 MB. Anyone sizing a SLURM –mem or Nextflow memory directive from the mapping number alone would be OOM-killed.

What the numbers mean, exactly – because resource.getrusage offers no way to be more precise, and a vaguer definition here would be a lie rather than a simplification:

peak_rss_mb

The whole process (plus reaped children) high-water mark as of the end of this stage. Monotone across stages: it can only rise. getrusage reports high-water marks only – RUSAGE_SELF since process start, RUSAGE_CHILDREN cumulative over children – and there is no per-stage reset, so a stage cannot be charged its own peak in isolation when all three run in one process (arda rnaseq / arda amplicon).

rss_gain_mb

How much this stage raised that mark; 0 if it stayed under an earlier stage’s peak.

The monotone number is the one an operator actually needs: it is what the process required at that point, which is what a memory directive has to cover. For per-stage attribution, run the stage in its own process (arda map|assemble|correct separately) – then peak_rss_mb is that stage alone.

arda.rnaseq._res.peak_rss_mb()

Peak RSS of this process AND its children, in MB; 0.0 where resource is absent.

RUSAGE_SELF alone is wrong and was: 92 % of a map run’s wall time is spent inside the mmseqs subprocess, whose nucleotide prefilter allocates a 4**k index table that dominates the footprint. Reporting only the Python process understated peak RSS by roughly an order of magnitude. ru_maxrss is bytes on macOS, KB on Linux.

⛔ Lives HERE, in the module with no arda imports, and is re-exported by arda.rnaseq._res. The other direction is an import cycle: arda.rnaseq.__init__ imports map, which needs Throttle.

Return type:

float

class arda.rnaseq._res.Stage

Bases: object

Times a stage and records its resource footprint onto a report dataclass.

Not a context manager: correct_airr and assemble_contigs write their JSON report inside the function, so measurement has to close before that write, not after the block.

stage = Stage() … stage.finish(report) # sets wall_seconds / peak_rss_mb / rss_gain_mb if report_path: …

property wall_seconds: float

finish(report=None)

Stamp the elapsed time and RSS onto report (if it has those fields).

Run QC and logging

Run QC: one long-format TSV describing a library, its reads and its clonotypes.

arda stats reads the artifacts a run already wrote – the Stage-1 AIRR, the clonotype table, the .arda.json report – and emits every number an operator needs to decide whether a sample is usable, without re-reading the FASTQ. It adds no alignment and no reference lookup beyond the germline gene list.

Long format, four columns, scope / key / metric / value. The scopes, and what keys them:

run

keyed by stage (map / correct / assemble) – the run report, flattened verbatim

sample

unkeyed – library-wide totals, junction lengths and quality, SHM rate, gene coverage

chain

keyed by locus (TRB, IGH, …) – reads AND clonotypes

v_gene / j_gene

keyed by gene (TRBV19) – reads and clonotypes per germline gene

allele_candidate

keyed by allele:mutation – a recurrent, high-quality V mutation

⛔ Long, not wide, and deliberately: the metric set differs per scope (a gene has no junction length, a chain has no allele frequency), so a wide table would be mostly empty cells, and the one thing a QC table must support is grep / join / a per-metric plot across samples. One value per cell, one row per fact – no 134/62 hybrids, and integers stay integers.

⛔ The chimera, non-functional and stop-codon counts are FLAGS, not filters. Nothing here removes a row from any output; stats only reads. See correct --flag-chimeras for why the chimera signature cannot separate a true PCR artefact from two real clones sharing a prefix and a suffix.

⚠ Alleles vs SHM is a heuristic, and it is reported as one. A mutation seen in most of an allele’s reads at high Phred is far more likely a germline the reference does not carry than somatic hypermutation or a miscall – but arda does not genotype, and allele_candidate is a shortlist to look at, never a call. The thresholds are exposed (--allele-min-frac, --allele-min-reads) precisely so the number can be re-derived rather than trusted.

arda.stats.collect(*, airr=None, clones=None, report=None, r1=None, r2=None, organism='human', allele_min_frac=0.5, allele_min_reads=10)

Every statistic arda can derive from a finished run, as (scope, key, metric, value).

Every input is optional and each contributes its own scopes, so this works on a bare arda annotate output as well as on a full arda rnaseq run directory.

Parameters:

• airr (str | Path | None) – Stage-1 (or annotate) AIRR TSV -> the chain read rows, *_gene reads, and allele_candidate.

• clones (str | Path | None) – clonotype table -> the chain clonotype rows and *_gene clonotypes.

• report (str | Path | None) – <prefix>.arda.json (or a bare --report JSON) -> the run scope, which is where total/mapped reads, threads, wall time and peak RSS come from.

• r1 (str | Path | None) – the input FASTQs. Used ONLY for their size on disk and for whether the library is paired – neither is recoverable from the AIRR, which holds the mapped subset.

• r2 (str | Path | None) – the input FASTQs. Used ONLY for their size on disk and for whether the library is paired – neither is recoverable from the AIRR, which holds the mapped subset.

• organism (str)

• allele_min_frac (float)

• allele_min_reads (int)

Return type:

list[tuple]

arda.stats.write_stats(rows, output)

Write rows as the QC TSV. Returns the row count.

Parameters:

• rows (list[tuple])

• output (str | Path)

Return type:

int

arda.stats.ALLELE_MIN_FRAC = 0.5

A V mutation is a candidate allele when it is carried by at least this fraction of the reads calling that allele AND by at least ALLELE_MIN_READS of them. The fraction is what separates germline from SHM (hypermutation is per-clone, so it does not reach half an allele’s reads); the count is what keeps a 2-read allele from producing a candidate off one read.

Logging, progress and resource reporting.

One arda logger, configured once from the CLI callback, that every module already feeds: cdr3fix, prefilter, segmap, refbuild.* all call logging.getLogger(__name__), so they are children of arda and inherit whatever setup installs. Nothing else has to know about verbosity.

⛔ Progress goes to stderr, results go to stdout. The stage lines used to be typer.echo on stdout, which means arda export-ref piped into a file interleaved a progress line with the data. Everything informational is a log record now; only paths and the export payload stay on stdout.

--log-file is always DEBUG whatever the console level is, and its format carries a timestamp and the process peak RSS. That is the artifact a cluster job leaves behind, and re-running a 10-hour bulk sample because the console was at the default level is not a thing anyone should have to do.

arda._log.setup(verbosity=0, quiet=False, log_file=None)

Configure the arda logger. Idempotent – handlers are replaced, never stacked.

Parameters:

• verbosity (int) – 0 = INFO (the stage/progress lines), 1+ = DEBUG with level and module names.

• quiet (bool) – WARNING and above only. Loses to --log-file, which stays at DEBUG.

• log_file (str | Path | None) – also write DEBUG records here, with timestamps and peak RSS per line.

Return type:

None

class arda._log.Throttle(seconds=30.0)

Bases: object

One-every-seconds gate for a progress line inside a hot loop.

A bulk map run flushes hundreds of chunks; logging each at INFO turns a 100 M-read sample into 500 lines of noise, and logging none leaves a multi-hour job with no sign of life. Time is the right axis, not chunk count – chunk wall time varies ~50x with the receptor fraction.

Parameters:

seconds (float)

ready()

Return type:

bool

arda._log.peak_rss_mb()

Peak RSS of this process AND its children, in MB; 0.0 where resource is absent.

RUSAGE_SELF alone is wrong and was: 92 % of a map run’s wall time is spent inside the mmseqs subprocess, whose nucleotide prefilter allocates a 4**k index table that dominates the footprint. Reporting only the Python process understated peak RSS by roughly an order of magnitude. ru_maxrss is bytes on macOS, KB on Linux.

⛔ Lives HERE, in the module with no arda imports, and is re-exported by arda.rnaseq._res. The other direction is an import cycle: arda.rnaseq.__init__ imports map, which needs Throttle.

Return type:

float

Cluster sharding

Multi-node (SLURM) sharding for very large inputs.

Single-node runs already stream in bounded memory (see annotate.mapper). For cluster scale we split the input once into N shards, annotate each shard as an independent SLURM array task, then concatenate the per-shard AIRR TSVs:

arda cluster split-fasta big.fastq work/shards –shards 50 # SLURM array task i: arda annotate -i work/shards/shard_<i>.fasta -o work/out/out_<i>.tsv … arda cluster merge work/out big.airr.tsv

arda slurm writes (and optionally submits) a single submit.sh that chains all three with SLURM job dependencies — see render_submit_script.

arda.cluster.split(input, out_dir, shards, *, prefix='shard')

Round-robin split a FASTA/FASTQ into shards FASTA files (one pass).

Round-robin (record k → shard k % shards) balances load even when record sizes vary, and every record lands in exactly one shard.

Parameters:

• input (str | Path)

• out_dir (str | Path)

• shards (int)

• prefix (str)

Return type:

list[Path]

arda.cluster.split_pairs(r1, out_dir, *, shards, r2=None, prefix='shard')

Split FASTQ into shards CONTIGUOUS blocks of read pairs, byte for byte.

Not split(). That one writes FASTA — dropping the quality string merge_pair’s per-base tie-break needs under --reconstruct — and round-robins records, which puts a fragment’s two mates in different shards. Mate separation is not a hypothetical defect here: it produced a published false discovery in this project’s own data (a spurious “R2-only blind spot”) that had to be retracted.

Contiguous, not round-robin, and that is load-bearing. Concatenating the per-shard Stage-1 AIRR in shard order then reproduces the single-node row order exactly, so Stage 2 and Stage 3 see byte-identical input and the sharded result is byte-identical to a single-node run. Round-robin would only give a permutation, and the clonotype fold is not permutation-invariant (correct collapses error children onto the parent it meets first, and coverage assignment is first-with-longest-overlap-wins).

R1 and R2 are cut at the same record boundaries, so mate k always lands in the same shard as mate k; no read ids are parsed to achieve it. map re-checks the pairing on every shard anyway.

Parameters:

• r1 (str | Path) – FASTQ (optionally gzipped). Single-end if r2 is None.

• out_dir (str | Path) – written as <prefix>_00000_R1.fastq (+ _R2 when paired).

• shards (int) – number of blocks; a shard receiving no records is not written.

• r2 (str | Path | None) – the mate file.

• prefix (str) – shard file stem.

Returns:

(r1_path, r2_path | None) per non-empty shard, in shard order.

Raises:

ValueError – on FASTA input, a bad shard count, an empty input, or mates of different lengths (a truncated R2, caught here rather than after N wasted tasks).

Return type:

list[tuple[Path, Path | None]]

arda.cluster.merge(shard_outputs, output)

Concatenate per-shard AIRR TSVs into one, keeping a single header.

shard_outputs may be a directory (its *.tsv are merged in sorted order) or an explicit list of files.

Parameters:

• shard_outputs (str | Path | list[Path])

• output (str | Path)

Return type:

Path

arda.cluster.render_submit_script(input, output, work_dir, *, shards, organism='human', seqtype='nt', threads=8, strand='both', map_d=True, partition=None, time='04:00:00', mem='8G', arda_mmseqs=None)

Render a submit.sh that chains split → array-annotate → merge on SLURM.

Uses sbatch --array with --wrap and an afterok dependency so the merge runs only once every shard succeeds. arda_mmseqs (if given) is exported so array tasks find the binary.

Parameters:

• input (str | Path)

• output (str | Path)

• work_dir (str | Path)

• shards (int)

• organism (str)

• seqtype (str)

• threads (int)

• strand (str)

• map_d (bool)

• partition (str | None)

• time (str)

• mem (str)

• arda_mmseqs (str | None)

Return type:

str

arda.cluster.render_rnaseq_submit_script(r1, out_prefix, work_dir, *, shards, r2=None, out_dir='.', organism='human', threads=8, kmer=12, min_score=75.0, reconstruct=False, assemble=True, complete_only=True, map_d=True, partition=None, time='04:00:00', mem='8G', reduce_time='08:00:00', reduce_mem='16G', arda_mmseqs=None)

Render a submit.sh chaining split → array-map → reduce for paired RNA-seq.

A sibling of render_submit_script() rather than a generalisation of it. That one chains split → array-annotate → merge over a single FASTA, and its last step is a pure concatenation. This chain differs in every step: the shard unit is a read pair, the files are FASTQ with quality, and the last step is a reduce — merge, then assemble and correct once over the whole merged AIRR. Folding both into one renderer would mean a parameter deciding which of two unrelated pipelines you get.

Only Stage 1 is distributed. correct collapses error variants and counts distinct fragments globally, and assemble grows contigs across reads, so sharding either would double-count clones and silently drop exactly the long-CDR3 contigs Stage 3 exists to build.

Two details in the array body that are not cosmetic:

• printf "%05d" — shard names are zero-padded so sorted() is numeric. The merge concatenates in name order and that order is read order.

• [ -s "$f" ] || exit 0 — a shard with no reads must not fail its task, or the afterok dependency drops the whole reduce step. (split_pairs does not write empty shards; this is the belt to its braces, for a resubmitted or hand-edited array range.)

Parameters:

• r1 (str | Path)

• out_prefix (str)

• work_dir (str | Path)

• shards (int)

• r2 (str | Path | None)

• out_dir (str | Path)

• organism (str)

• threads (int)

• kmer (int)

• min_score (float)

• reconstruct (bool)

• assemble (bool)

• complete_only (bool)

• map_d (bool)

• partition (str | None)

• time (str)

• mem (str)

• reduce_time (str)

• reduce_mem (str)

• arda_mmseqs (str | None)

Return type:

str

Reference build

Orchestrate the per-species reference database build.

For each locus: enumerate deduplicated V-J scaffolds, annotate them with IgBLAST, keep those with complete FR1-FR4 + CDR1-3 markup, translate to protein, and derive protein markup. Writes the committed artifacts under database/vdj/<organism>/ plus a comprehensive build.log.

arda.refbuild.build.build(organism='all', *, one_allele_per_gene=False, allow_chimeras=False)

Build one organism or "all" supported organisms.

Parameters:

• organism (str)

• one_allele_per_gene (bool)

• allow_chimeras (bool)

Return type:

None

arda.refbuild.build.build_species(organism, *, one_allele_per_gene=False, allow_chimeras=False)

Build the reference DB for one organism. Returns the output directory.

one_allele_per_gene builds scaffolds from a single representative allele per gene (*01 where it exists, else the lowest-numbered) – roughly a 4x smaller reference with no allele-level ambiguity. Off by default.

allow_chimeras additionally builds the TRDV × TRAJ scaffolds the default refuses – see loci_for() for the measurement that makes this a live question rather than a settled one. Off by default.

Parameters:

• organism (str)

• one_allele_per_gene (bool)

• allow_chimeras (bool)

Return type:

Path

Enumerate in-frame V-J reference scaffolds.

For markup transfer the FR/CDR region coordinates are fully determined by the V gene (FR1-3, CDR1-2, CDR3 start at the conserved Cys104) and the J gene (CDR3 end, FR4). The D segment lies inside the hypervariable CDR3 — query-specific at runtime — so we enumerate V×J scaffolds for every locus and, for VDJ loci, insert a short frame-neutral N spacer where D would sit so IgBLAST still annotates a plausible CDR3 + FR4.

Each scaffold is V + N*pad + J where pad keeps the J coding frame aligned to V’s reading frame (jframe from the IgBLAST aux file). Byte-identical scaffolds are deduplicated: one DB entry, with all contributing (V,J) allele pairs recorded.

class arda.refbuild.combinations.Scaffold(scaffold_id, locus, sequence, v_calls=<factory>, j_calls=<factory>, n_pad=0)

Bases: object

A deduplicated V-J reference scaffold.

Fields: scaffold_id (stable "{locus}_{index}"), locus, sequence (assembled V + N*pad + J), v_calls / j_calls (all alleles producing this scaffold), and n_pad (N nucleotides between V and J).

Parameters:

• scaffold_id (str)

• locus (str)

• sequence (str)

• v_calls (list[str])

• j_calls (list[str])

• n_pad (int)

scaffold_id: str

locus: str

sequence: str

v_calls: list[str]

j_calls: list[str]

n_pad: int = 0

arda.refbuild.combinations.load_j_frames(organism)

Parse bin/optional_file/<organism>_gl.aux -> {J allele: frame}.

Frame is the 0-based “first coding frame start position” (column 2).

Parameters:

organism (str)

Return type:

dict[str, int]

arda.refbuild.combinations.load_j_fr4_offsets(organism)

Parse the same aux file -> {J allele: (cdr3_stop, extra_bp)}, both 0-based nt counts.

IgBLAST’s aux carries five columns: allele, coding-frame start, chain type, CDR3 stop, and extra bp beyond the J coding end. arda only ever read column 2. Columns 4 and 5 pin FR4 inside the J exactly:

fwr4 = j_seq[cdr3_stop + 1 : len(j_seq) - extra_bp]

That is how a J + C scaffold gets an FR4 at all: igblastn cannot annotate a V-less sequence, so those scaffolds are not routed through it (see refbuild.build).

Verified against every V-J scaffold arda builds: the string this yields is byte-identical to IgBLAST’s own fwr4 on all 125 human J alleles where both exist – including IgBLAST’s own non-multiple-of-3 cases (IGHJ6*02 has extra_bp = 0 and a 34 nt FR4; 726 V-J scaffolds already carry one). Reproducing IgBLAST, quirks included, is the requirement here: the two scaffold kinds must agree, or a J->C read and a V-J read of the same clone disagree on FR4.

Pseudogene J entries carry only three columns and are skipped – they have no FR4 to report.

Parameters:

organism (str)

Return type:

dict[str, tuple[int, int]]

arda.refbuild.combinations.load_v_fwr3_stops(organism)

Parse bin/internal_data/<organism>/<organism>.ndm.imgt -> {V allele: FWR3 stop}.

IgBLAST’s own IMGT annotation of its V germlines. Column 11 is the 1-based FWR3 stop, and FR3-IMGT ends at position 104 – the conserved 2nd-CYS – so the Cys104 codon starts at fwr3_stop - 3 (0-based). This is the authoritative V junction anchor, and it is the same IgBLAST metadata the rest of the build already trusts.

It does NOT cover every IMGT allele (IgBLAST ships a subset), hence the motif fallback in refbuild.build._v_anchor.

Parameters:

organism (str)

Return type:

dict[str, int]

arda.refbuild.combinations.build_locus_scaffolds(locus, v_alleles, j_alleles, j_frames, *, d_spacer=None)

Build deduplicated V×J scaffolds for one locus.

Parameters:

• locus (Locus) – The locus definition.

• v_alleles (dict[str, str]) – {allele: ungapped V sequence}.

• j_alleles (dict[str, str]) – {allele: ungapped J sequence}.

• j_frames (dict[str, int]) – {J allele: 0-based coding frame} from the aux file.

• d_spacer (int | None) – N spacer length for VDJ loci (default DEFAULT_D_SPACER_NT); forced to 0 for VJ loci.

Returns:

Scaffolds, one per unique assembled sequence.

Return type:

list[Scaffold]

Build a SEGMENT reference: V, J and J+C as separate targets, not a V×J product.

The shipped reference enumerates every V×J combination — 15,069 scaffolds from 775 V alleles and 124 J alleles. That costs, and it costs twice:

• Speed. A read covering only V aligns against every scaffold carrying that V — a median of 13, and 67 for TRA. Measured on a TRA amplicon: 277 gapped alignments per hitting read, one kept. Aligning the same reads against 1,244 segment targets instead is 6.9× faster (TRA), 7.7× (TRB), 2.5× (bulk RNA-seq).

• Accuracy. 81 % of hitting amplicon reads sit at the –max-seqs 300 cap, so the true scaffold is sometimes never even a candidate; and the V call is decided by a whole-scaffold bit score whose J half is arbitrary. Scored against IgBLAST truth on TRA, the segment reference takes V-gene concordance 99.00 % → 99.99 % (95 errors → 1) and J-gene 98.47 % → 99.90 %.

So this is not a speed/accuracy trade — the product reference was losing on both.

Derived from the built reference (markup.tsv + alleles.fasta), not from IMGT, so it needs no download and is reproducible from any checkout that can already map.

Coordinates carry over almost for free, which is why this is cheap to build correctly:

• a V target is scaffold[:v_sequence_end], and fwr1/cdr1/fwr2/cdr2/fwr3 are already scaffold-relative, so they transfer unchanged. cdr3 is truncated at the V end (the read only ever sees the V-side stub of the junction here).

• a J target is scaffold[j_sequence_start-1:vj_end], so every coordinate shifts by j_sequence_start - 1. It carries fwr4 and the J-side stub of cdr3.

• a C target is scaffold[vj_end:] of a J+C scaffold — the constant region alone, one per distinct C allele. It carries no regions (a constant region has none of fwr1..fwr4) and no V or J call, only c_call.

The C side was the same cross-product, and this module used to leave it in place. The 345 J+C scaffolds were copied through verbatim, and they are a J×C product (IGH 14 J × 11 C, IGL 9 × 7, TRB 16 × 2) in which every scaffold of a locus ends in the same constant sequence. So a read reaching C was aligned against all of them to learn one c_call, at a redundancy factor equal to the locus’ J-allele count — 69× on TRA. Measured on a TRA amplicon: 345 of 1,244 targets (27.7 %) produced 76.4 % of all segment alignments, 4,977 alignments per target against 603 for a V target. Splitting them into the existing J| targets plus 25 C| targets takes the segment search to 1.89× and the alignment count to 4.23× fewer, with the V-and-J fast path down 0.81 %.

⚠ What this does NOT do. A read that spans the junction hits a V target and a J target separately, and something must merge the two into one AIRR record. That is 85.6 % of mapped amplicon reads but only 7.2 % of bulk RNA-seq reads, which is why the RNA-seq path is nearly free and the amplicon path is not. The merge lives in the mapper, not here; this module only builds the targets.

arda.refbuild.segments.build_segment_reference(organism='human', *, out_dir=None)

Write segments.fasta + segments.markup.tsv beside the reference.

One target per distinct (segment, allele): the longest scaffold instance of that allele is used as the donor, so a V target is never accidentally truncated by whichever J it happened to be paired with.

Returns:

SegmentStats — target counts and the reduction factor vs the V×J reference.

Parameters:

• organism (str)

• out_dir (Path | None)

Return type:

SegmentStats

class arda.refbuild.segments.SegmentStats(v_targets=0, j_targets=0, c_targets=0, source_scaffolds=0)

Bases: object

What the build produced, for the report and for tests to assert on.

Parameters:

• v_targets (int)

• j_targets (int)

• c_targets (int)

• source_scaffolds (int)

v_targets: int = 0

j_targets: int = 0

c_targets: int = 0

source_scaffolds: int = 0

property total: int

property reduction: float

as_dict()

Return type:

dict

Native nucleotide translation and reading-frame utilities.

No BioPython. The hot functions (translate, detect_coding_frame, reverse_complement, back_translate) are implemented in the C++ extension arda._markup and re-exported here; a pure-Python fallback keeps the module importable if the extension is unavailable. These mirror mirpy’s mirseq API so mirpy can later import arda and reuse them.

arda.refbuild.translate.translate(nt, frame=0)

Translate a nucleotide string from frame (0/1/2).

Parameters:

• nt (str)

• frame (int)

Return type:

str

arda.refbuild.translate.detect_coding_frame(nt)

Return the reading frame (0/1/2) with the fewest stop codons.

Parameters:

nt (str)

Return type:

int

arda.refbuild.translate.reverse_complement(nt)

Reverse-complement a nucleotide string (non-ACGT -> N).

Parameters:

nt (str)

Return type:

str

arda.refbuild.translate.back_translate(aa, unknown='NNN')

Mock back-translation via most-frequent human codons.

Parameters:

• aa (str)

• unknown (str)

Return type:

str

arda.refbuild.translate.aa_coords_from_nt(nt_start, nt_end, coding_start)

Map a 1-based closed nt interval to 1-based closed aa coordinates.

Parameters:

• nt_start (int) – 1-based start of the region in the nucleotide sequence.

• nt_end (int) – 1-based end (closed).

• coding_start (int) – 1-based nt position where translation begins (frame origin).

Returns:

(aa_start, aa_end) 1-based closed, in the translated protein.

Return type:

tuple[int, int]

IMGT/V-QUEST germline reference download, parsing, and ungapping.

The IMGT V-QUEST reference directory ships gapped germline FASTAs laid out as <Species>/<IG|TR>/<GENE>.fasta (e.g. Homo_sapiens/IG/IGHV.fasta). Sequences carry IMGT-numbering gap dots; IgBLAST’s edit_imgt_file.pl ungaps them and rewrites headers to bare allele names (what makeblastdb wants).

This module:

• downloads & extracts the reference zip into data/imgt (gitignored),

• parses the original gapped FASTA headers for per-allele functionality,

• ungaps a gene file via edit_imgt_file.pl into data/imgt/ungapped.

class arda.refbuild.imgt.ImgtAllele(allele, functionality, sequence)

Bases: object

A germline allele parsed from an IMGT FASTA header + sequence.

Parameters:

• allele (str)

• functionality (str)

• sequence (str)

allele: str

functionality: str

sequence: str

property is_functional: bool

arda.refbuild.imgt.download_reference(*, force=False)

Download and extract the IMGT V-QUEST reference directory.

Returns the extraction root (containing the per-species directories). Idempotent unless force.

Parameters:

force (bool)

Return type:

Path

arda.refbuild.imgt.gene_fasta_path(species_dir, group, gene_stem)

Path to a gene-type FASTA, e.g. Homo_sapiens/IG/IGHV.fasta.

Handles the occasional top-level wrapper directory inside the zip.

Parameters:

• species_dir (str)

• group (str)

• gene_stem (str)

Return type:

Path

arda.refbuild.imgt.parse_functionality(path)

Map allele name -> normalized functionality from gapped IMGT headers.

IMGT header: accession|allele|species|functionality|region|.... The functionality field may be wrapped, e.g. (F) / [F] for inferred.

Parameters:

path (Path)

Return type:

dict[str, str]

arda.refbuild.imgt.ungap_gene(species_dir, group, gene_stem)

Ungap a gene file with edit_imgt_file.pl; return the ungapped path.

Parameters:

• species_dir (str)

• group (str)

• gene_stem (str)

Return type:

Path

arda.refbuild.imgt.read_fasta(path)

Read a FASTA file into (header, sequence) pairs (sequence joined).

Parameters:

path (Path)

Return type:

list[tuple[str, str]]

External tool wrappers

Thin wrapper around the mmseqs binary.

Inspired by pymmseqs (MIT) but deliberately dependency-free: we only need binary discovery, a subprocess runner, and the createdb / search / convertalis (and easy-search) pipeline used by the annotator.

Discovery order for the binary: $ARDA_MMSEQS → the optional arda-mmseqs companion wheel → <project>/bin/mmseqs → mmseqs on PATH. Candidates after the explicit override are version-matched against the precompiled indexes in database/: an index is only reusable by the mmseqs release that built it, so accepting an arbitrary PATH binary silently discards the shipped index and rebuilds a private cache. If nothing matches, a known-good static binary is auto-fetched into <project>/bin/mmseqs (one-time, transparent) unless $ARDA_NO_AUTO_FETCH is set — so neither pip nor conda users need to install mmseqs manually – the binary is fetched on first use.

exception arda.mmseqs.MMseqsError

Bases: RuntimeError

Raised when an mmseqs invocation exits non-zero.

arda.mmseqs.mmseqs_binary()

Locate an mmseqs executable that can actually use the shipped indexes.

Resolution: $ARDA_MMSEQS → the arda-mmseqs companion wheel → <project>/bin/mmseqs → mmseqs on PATH → auto-fetched static build.

Version-matched, not merely present. Taking whatever mmseqs happened to be on PATH was a silent correctness and performance bug: an index is only reusable by the release it was compiled with, so a mismatched binary makes every run reject database/’s precompiled DBs and rebuild a private cache instead — no error, just a slow start and, if the two releases align differently, results that are not comparable with anyone else’s. Found in the wild: a cluster with a bare-git-hash build ahead of conda’s on PATH.

$ARDA_MMSEQS is never version-checked — an explicit override is the user’s call. If nothing matches, the known-good static build is fetched (unless $ARDA_NO_AUTO_FETCH); if that also fails we fall back to the best candidate and warn, naming the consequence.

Return type:

str

arda.mmseqs.run(args, *, check=True)

Run mmseqs <args> capturing stdout/stderr.

Parameters:

• args (list[str])

• check (bool)

Return type:

CompletedProcess

arda.mmseqs.version()

Return the mmseqs version string (e.g. 18-8cc5c).

Return type:

str

arda.mmseqs.version_key(v)

Canonical form of a version string, for comparing an index marker to the running mmseqs.

mmseqs version prints the same release+commit with different punctuation across builds: the official static binary says 18-8cc5c, the bioconda build 18.8cc5c. They are the same mmseqs and produce byte-compatible indexes; only the separator differs. Comparing the raw strings with == therefore rejected every committed index the moment the toolchain moved from conda’s mmseqs to the static one – so the precompiled DBs shipped in database/ were never used and every run rebuilt a private cache instead.

Fold each run of separator characters to a single - and lowercase. This bridges the cosmetic difference while still distinguishing genuinely different versions (17-b804f != 18-8cc5c), so an incompatible index is never accepted.

Not sufficient on its own – see versions_compatible(), which is what callers should use.

Parameters:

v (str)

Return type:

str

arda.mmseqs.versions_compatible(a, b)

Do two mmseqs version strings denote builds with interchangeable index formats?

Punctuation is not the only way the same build spells itself. The official static release asset prints its full 40-char commit hash (8cc5ce367b5638c4306c2d7cfc652dd099a4643f) while the bioconda build and the committed index marker print release+short-commit (18.8cc5c / 18-8cc5c). Release 18 is commit 8cc5c..., so those are one build – but no amount of separator folding makes the strings equal.

That mattered concretely: arda’s own auto-fetched binary is the static asset, so a pure version_key() comparison rejected the index arda itself ships, on every macOS install.

So: compare commit hashes when both carry one, accepting a prefix match in either direction (short vs full form). Fall back to version_key() when one has no hash at all. A genuinely different build has a different commit (76da68ad... is not 8cc5c...) and is still rejected.

Parameters:

• a (str)

• b (str)

Return type:

bool

arda.mmseqs.createdb(fasta, db, *, dbtype=None)

Create an mmseqs sequence DB from a FASTA file.

dbtype: None auto-detect, 1 amino-acid, 2 nucleotide.

Parameters:

• fasta (str | Path)

• db (str | Path)

• dbtype (int | None)

Return type:

Path

arda.mmseqs.search(query_db, target_db, result_db, tmp_dir, *, search_type=0, sensitivity=5.7, evalue=0.001, max_seqs=300, threads=1, kmer=None, extra=None)

Run mmseqs search with backtrace enabled (-a).

Parameters:

• kmer (int | None) – MMseqs2 -k. This is the memory knob. The nucleotide prefilter allocates a k-mer index table of 4**k entries, so the default k=15 costs 4**15 * 8 B ~ 8.6 GB regardless of database size, thread count or chunk size. k=13 costs ~0.7 GB. Leave None for MMseqs2’s own default.

• query_db (str | Path)

• target_db (str | Path)

• result_db (str | Path)

• tmp_dir (str | Path)

• search_type (int)

• sensitivity (float)

• evalue (float)

• max_seqs (int)

• threads (int)

• extra (list[str] | None)

Return type:

Path

arda.mmseqs.convertalis(query_db, target_db, result_db, out_tsv, *, format_output='query,target,qstart,qend,tstart,tend,qlen,tlen,alnlen,mismatch,gapopen,cigar,qaln,taln,evalue,bits,pident', threads=1, search_type=None)

Convert an alignment result DB to a TSV with the requested columns.

search_type must be passed for nucleotide results (3) so convertalis can interpret the alignment; otherwise mmseqs cannot tell nt from translated.

Parameters:

• query_db (str | Path)

• target_db (str | Path)

• result_db (str | Path)

• out_tsv (str | Path)

• format_output (str)

• threads (int)

• search_type (int | None)

Return type:

Path

arda.mmseqs.top_hit(result_db, out_db)

Reduce an alignment DB to the single best-scoring hit per query.

MMseqs2 already stores each query’s results sorted by descending score, so taking the first line per entry is the best hit – this is the idiom mmseqs’ own filterdb usage message shows. Verified on 100 k reads against the human reference: 4,101 queries, identical target and identical bit score to a full sort-and-dedupe in polars, on every one.

Why it matters: with –max-seqs 300, 4,101 hitting queries produced 804,341 alignment rows (194 MB of TSV, each row carrying cigar/qaln/taln). Parsing that dominated arda’s peak RSS – 877 MB, against 284 MB for the mmseqs subprocess itself. Reducing before convertalis writes 1.0 MB instead, and costs 0.04 s.

Parameters:

• result_db (str | Path)

• out_db (str | Path)

Return type:

Path

arda.mmseqs.easy_search(query_fasta, target_fasta_or_db, out_tsv, tmp_dir, *, search_type=0, sensitivity=5.7, evalue=0.001, max_seqs=300, threads=1, format_output='query,target,qstart,qend,tstart,tend,qlen,tlen,alnlen,mismatch,gapopen,cigar,qaln,taln,evalue,bits,pident', strand=None, extra=None)

One-shot createdb+search+convertalis producing a TSV.

strand (nucleotide search only): 1 forward, 2 both strands; None lets mmseqs default (forward).

Parameters:

• query_fasta (str | Path)

• target_fasta_or_db (str | Path)

• out_tsv (str | Path)

• tmp_dir (str | Path)

• search_type (int)

• sensitivity (float)

• evalue (float)

• max_seqs (int)

• threads (int)

• format_output (str)

• strand (int | None)

• extra (list[str] | None)

Return type:

Path

Wrapper around an IgBLAST release, fetched on demand.

Used at build time (Phase 1) to construct the curated reference DB, and by arda igblast, which is how every gold-standard comparison in the benchmark is produced. The runtime annotator does not depend on IgBLAST.

The release is a flat directory – executables plus the internal_data and optional_file trees, with $IGDATA pointed at it:

igblastn  igblastp  makeblastdb  edit_imgt_file.pl
internal_data/   optional_file/

Resolved by igblast_root(), in order:

1. $ARDA_IGBLAST – an explicit directory, never fetched;

2. <project>/bin if it already holds an IgBLAST (what setup.sh produces in a checkout);

3. <project>/igblast, auto-fetched from NCBI on first use.

Step 3 is why this module changed shape. It used to resolve only through <project>/bin, so a plain pip install arda-mapper – which has no checkout and never runs setup.sh – could not run arda igblast at all. The failure was also misleading: it surfaced as IgBlastError: IgBLAST ships no internal annotation for organism 'human', which names a missing data file rather than a missing install, and so reads as a broken reference rather than as “IgBLAST was never installed here”.

exception arda.igblast.IgBlastError

Bases: RuntimeError

Raised when an IgBLAST tool invocation fails or is missing.

arda.igblast.igblast_root()

The directory holding the IgBLAST executables and internal_data.

Fetches a release on first use; $ARDA_NO_AUTO_FETCH turns that into an error instead.

Return type:

Path

arda.igblast.igblast_version()

The auto-fetched IgBLAST release version, or None if it was not fetched by arda.

Return type:

str | None

arda.igblast.igdata_env()

Environment with IGDATA pointing at the IgBLAST data root.

Return type:

dict[str, str]

arda.igblast.tool(name)

Resolve an IgBLAST executable, fetching the release if it is not installed yet.

Parameters:

name (str)

Return type:

Path

arda.igblast.edit_imgt_file(imgt_fasta, out_fasta)

Ungap an IMGT germline FASTA via edit_imgt_file.pl.

Parameters:

• imgt_fasta (str | Path)

• out_fasta (str | Path)

Return type:

Path

arda.igblast.makeblastdb(in_fasta, out_db, *, dbtype='nucl')

Build a germline BLAST database from an ungapped FASTA.

Parameters:

• in_fasta (str | Path)

• out_db (str | Path)

• dbtype (str)

Return type:

Path

arda.igblast.igblastn_airr(query_fasta, out_tsv, *, organism, germline_db_v, germline_db_j, germline_db_d=None, auxiliary_data=None, ig_seqtype='TCR', num_threads=1)

Run igblastn -outfmt 19 (AIRR rearrangement TSV).

Parameters:

• query_fasta (str | Path)

• out_tsv (str | Path)

• organism (str)

• germline_db_v (str | Path)

• germline_db_j (str | Path)

• germline_db_d (str | Path | None)

• auxiliary_data (str | Path | None)

• ig_seqtype (str)

• num_threads (int)

Return type:

Path

arda.igblast.auxiliary_data(organism)

optional_file/<organism>_gl.aux — IgBLAST’s J-gene coding frames.

⛔ Without this file IgBLAST silently emits no CDR3 and no junction. It is what tells igblastn each J allele’s reading frame, and with no frame there is nothing to place the Phe/Trp 118 anchor against. Everything else still works: V and J are called, v_score is normal, the process exits 0 — only cdr3*, junction and junction_aa come back empty, on every read.

It lives beside the executables, under igblast_root(). Both callers used to look under paths.bin_dir() and then fall back to passing nothing when the file was not there. Those two are THE SAME directory in a source checkout (setup.sh installs IgBLAST into <repo>/bin), so it worked everywhere it was developed and failed on every auto-fetched install, where the root is $XDG_CACHE_HOME/arda/igblast while bin_dir() is $XDG_CACHE_HOME/arda/bin. Measured cost: a 10,000-read amplicon IgBLAST truth carrying j_call on 9,070 of 9,300 reads and junction_aa on zero — written up as an IgBLAST limitation at 151 bp before it was traced to here.

Raises:

IgBlastError – if the file is absent. A missing frame table must not degrade to “no junctions”: that is indistinguishable from a truth set which genuinely has none, which is exactly how this went unnoticed.

Parameters:

organism (str)

Return type:

Path

Cache layout and reference fetch

Where arda keeps its reference and how a plain pip install acquires one.

Filesystem layout discovery for arda.

Two modes, resolved transparently:

• Source checkout ($ARDA_HOME or an editable/development install) — bin/, data/ and database/ live next to pyproject.toml, exactly as committed.

• PyPI install (pip install arda-mapper, no source tree) — there is no bundled database/ in the wheel (it is 50+ MB of curated references), so everything lives under a per-user cache ($XDG_CACHE_HOME/arda or ~/.cache/arda). The curated reference is auto-fetched once into the cache on first use (see arda._database_fetch), and mmseqs target DBs are built there on demand. No ARDA_HOME needed.

Set $ARDA_NO_AUTO_FETCH to disable the one-time reference download (e.g. air-gapped runs where the cache was pre-populated).

arda.paths.project_root()

Source checkout root if present, else the per-user cache root.

Unlike older arda, this never raises: a PyPI install with no source tree resolves to the cache dir, where the reference is auto-fetched and mmseqs DBs are built.

Return type:

Path

arda.paths.cache_root()

Per-user cache root ($XDG_CACHE_HOME/arda or ~/.cache/arda).

Return type:

Path

arda.paths.bin_dir()

Directory holding the downloaded mmseqs/IgBLAST binaries (gitignored / cached).

Return type:

Path

arda.paths.data_dir()

Writable scratch directory for downloads and built mmseqs DBs (gitignored / cached).

Return type:

Path

arda.paths.database_dir()

Curated reference database root.

Source checkout → the committed database/. PyPI install → <cache>/database, with the reference auto-fetched on first call unless $ARDA_NO_AUTO_FETCH is set.

Return type:

Path

arda.paths.vdj_dir(species=None)

database/vdj (or a single species subdirectory).

Parameters:

species (str | None)

Return type:

Path

Auto-fetch the curated arda reference database (a GitHub release asset).

The wheel ships code only; the curated vdj/ reference (allele FASTAs + region markup, per species, AA + NT; ~50 MB on disk, ~3 MB compressed) is published as the arda-reference-vdj.tar.gz asset on the matching vX.Y.Z GitHub release and downloaded once into the per-user cache on first use. Version-sensitive precompiled mmseqs DBs are not shipped — they are built on demand from the fetched FASTAs into <cache>/data.

arda._database_fetch.reference_url(version=None)

Release-asset URL for the reference tarball.

Defaults to _REFERENCE_TAG (the release that published the current reference), not to the running package version — see the note there.

Parameters:

version (str | None)

Return type:

str

arda._database_fetch.fetch_database(dest, *, force=False, version=None)

Download + extract the reference (vdj/<org>/...) into dest; return dest.

Skips the download if dest/vdj already exists (unless force).

``vdj/`` must never be visible in a partial state, because its mere existence is the gate every other arda process tests (paths.database_dir) — a half-populated tree is not a slow download, it is a silent wrong answer. Two things guarantee that:

• a build_lock(), because arda is routinely run concurrently against the same cache (one Nextflow process per sample, one SLURM task per array index), and on first use in a fresh environment all of them would otherwise fetch into the same path at once;

• extraction into dest itself, then a single os.replace(). The extraction directory has to be a sibling of the target for that to be a rename: the previous code extracted into /tmp and called shutil.move, which silently degrades to a recursive copy across filesystems (/tmp and ~/.cache usually are different ones) — populating vdj/ file by file, in full view of every other process.

Extraction still rejects symlinks/hardlinks and any path escaping the staging dir.

Parameters:

• dest (Path)

• force (bool)

• version (str | None)

Return type:

Path