API reference

arda is usable as a library as well as a CLI. The single entry point most callers need is arda.annotate_sequences(), which takes a batch of sequences and returns AIRR record dicts. It invokes MMseqs2 and manages temporary search files internally:

import arda

records = arda.annotate_sequences(
    [("q1", "TGTGCCAGCAGCTTAGCGGGAGGGAACACCGGGGAGCTGTTTTTTGGA")],
    seqtype="nt",          # "nt" or "aa"
    organism="human",      # human | mouse | rat | rabbit | rhesus_monkey
    threads=0,             # automatic: allocated physical cores
    map_d=True,            # D segments for VDJ loci; works on amino-acid input too
)

Every input produces a record. Unmapped inputs retain sequence_id and sequence with empty annotation fields. Submit one batch instead of calling the mapper once per sequence. The batch is held in memory; use the streaming CLI for a whole FASTQ.

The fields each record carries are documented in Output files and columns. The modules below are the rest of the public surface, grouped as the package is.

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, *, threads=0, map_dd=True, major_only=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. Supports nucleotide and amino-acid input.

• map_dd (bool) – Search a tandem second D (default True); False retains single-D mapping.

• major_only (bool) – Search one major allele per gene (default); False searches all alleles.

• threads (int) – CPU budget; 0 uses all CPUs available to the process/job. An already initialized Polars pool exceeding this budget raises ValueError.

Returns:

A list of AIRR-style annotation record dicts, one per input sequence. Unmapped sequences retain their ID and sequence with empty annotation fields.

Reference provenance

arda.reference_index_identity(organism, seqtype="aa") resolves and validates an index without annotating queries. Its JSON-compatible result contains the resolved source_fasta and database_path, source_fasta_sha256, and a files_sha256 mapping from MMseqs payload basenames to SHA256 digests. The mapping includes split sequence/header files and indexes. schema_version, dbtype, mmseqs_version, organism and seqtype describe the identity.

Only ARDA’s own source module triggers automatic source-root discovery; explicit ARDA_HOME remains authoritative. Index reuse checks content, not timestamps. Legacy indexes remain untouched. New generations publish atomically and remain immutable, so concurrent readers retain coherent sequence/header/index files. Capture the API result before and after a batch and require equality to detect a reference change during that batch.

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, map_dd=True, d_max_evalue=None, shm='framework', complete_junction_nt=0, major_only=True)

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)

• map_dd (bool)

• d_max_evalue (float | None)

• shm (str)

• complete_junction_nt (int)

• major_only (bool)

Return type:

Path

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

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)

• map_dd (bool)

• complete_junction_nt (int)

• major_only (bool)

Return type:

list[dict]

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

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

Publishes immutable generations under database/vdj/<org>/mmseqs/<seqtype> with source and payload SHA256 identities. Reuses validated generations unless force; existing generations are preserved for concurrent readers.

Parameters:

• organism (str)

• force (bool)

Return type:

None

arda.annotate.mapper.reference_index_identity(organism, seqtype='aa', *, major_only=True)

Resolve and validate the exact searched reference index, returning SHA256 provenance.

The returned JSON-compatible mapping identifies the resolved FASTA, immutable database prefix and every sequence/header/index file, including split payloads. Resolving may build a missing or stale index; it never runs an annotation search.

Parameters:

• organism (str)

• seqtype (str)

• major_only (bool)

Return type:

dict

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, map_dd=True)

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 input-strand read, stored in the AIRR sequence field after mandatory sequencing DNA IUPAC-to-A normalization; 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).

Direct calls require a canonical DNA query/alignment consistent with the normalized submitted read and strand. Normalize before aligning; ambiguous precomputed alignments are rejected rather than reinterpreted.

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)

• map_dd (bool)

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, major_alleles=<factory>)

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)

• major_alleles (dict[str, str])

organism: str

seqtype: str

target_fasta: Path

entries: dict[str, RefEntry]

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

anchors: dict

major_alleles: dict[str, str]

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.

Never: 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.

Never: 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

load_segment_markup()

Pick up a segments.markup.tsv written AFTER this reference was loaded.

Never: The segment reference is GENERATED, not shipped, and it is generated on the map path (mapper._cached_segment_db) — i.e. after load_reference has already read entries. So on the FIRST run against a fresh reference the two-pass search found its segment targets and then resolved every one of them to None: segment_j_call handed the raw target name to the combination lookup (no_such_combination for every J→C read) and _cannot_reach_cys104 saw no CDR3 markup, so –v-only-on-segment routed nothing. Measured on tests/data/rnaseq_real: 185 of 313 AIRR rows differed from every later run, fast_fraction read 0.0771 against 0.1707, and the v_only_on_segment counter was absent from the report. Correct output, exit 0, none of the speed — and the first run is the one a new user, a fresh container and a fresh CI checkout all get.

Segment keys are re-read wholesale rather than merged: a pre-2.8.0 file already loaded into entries names its targets JC|<scaffold> where the regenerated one names them C|<allele>, so the stale keys must go or they outlive the file they came from. Base scaffold ids never contain | (that separator is what marks a segment target), so the two key spaces are exactly separable.

Return type:

None

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

Load major-only search templates; major_only=False loads the full resource.

Each call returns its own metadata object. Major derivatives have a separate content-addressed path, so full-reference metadata and index keys stay distinct.

Parameters:

• organism (str)

• seqtype (str)

• major_only (bool)

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, seqtype='nt')

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

Sequencing DNA IUPAC bases are always replaced with A. seqtype="aa" denotes protein input, which is preserved.

with_qual=True yields (id, sequence, qual) instead – the FASTQ Phred string, or None for FASTA (which has no quality). 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)

• seqtype (str)

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]

arda.annotate.io.normalize_nt(sequence, *, aligned=False)

Uppercase sequencing DNA and replace every ambiguous DNA IUPAC base with A.

This mandatory rule preserves positions and rejects invalid non-DNA symbols. aligned=True additionally permits alignment gaps (- and .), never raw read gaps. Never apply this sequencing-input rule to germlines.

Parameters:

• sequence (str)

• aligned (bool)

Return type:

str

arda.annotate.io.normalize_nt_frame(frame)

Normalize observed AIRR nucleotide columns, preserving references and metadata.

Re-translate existing region/junction amino acids only where their observed nucleotide sequence changed. Nulls, coordinates and quality strings survive.

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.

Never: 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.

Never: 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, v_end=None, j_start=None, major_only=True, *, map_dd=True)

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

map_dd=False skips the second-D search while retaining the best single D. The default preserves tandem-D calling.

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.

v_end / j_start supply the interior to search instead of deriving it here, and a caller that already knows the boundaries should pass them. v_end is the COUNT of leading junction nucleotides the V templates and j_start the 0-based index of the first J-templated one – arda.cdr3fix.Cdr3Markup.v_end_nt and j_start_nt exactly.

⚠ Deriving them here is an exact germline prefix/suffix match against the CALLED allele, which is right for a read but wrong for a junction whose nucleotides were inferred under a model keyed on a different allele of the same gene: the prefix then breaks at the first synonymous difference and the interior opens up inside the V, where a spurious D wins. CATSIRFTDTQYF placed a TRBD2 at nucleotide 6 that way.

Parameters:

• junction_nt (str)

• v_call (str)

• j_call (str)

• species (str)

• d_max_evalue (float | None)

• v_end (int | None)

• j_start (int | None)

• major_only (bool)

• map_dd (bool)

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, map_dd=True, d_max_evalue=None, major_only=True)

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)

• map_dd (bool)

• d_max_evalue (float | None)

• major_only (bool)

Return type:

list[dict]

arda.annotate.contig.merge_contig(contig, reference, *, map_d=True, map_dd=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)

• map_dd (bool)

Return type:

dict

arda.annotate.contig.merge_contigs(contigs, organism='human', seqtype='nt', *, reference=None, map_d=True, map_dd=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)

• map_dd (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.

Placement is one gapless local alignment per side (_markup.d_local_align, arda’s C++ D caller), and the outcome follows from where the hit starts in each sequence – VDJdb’s Cdr3Fixer table, unchanged:

start_in_segment  start_in_cdr3
0                 0              NoFixNeeded
0                 > 0            FixTrim    -- framework outside the junction
> 0               0              FixAdd     -- germline the submission cut
> 0               <= max_replace FixReplace -- substitute at the anchor
> 0               > max_replace  FailedReplace

Searching every offset is the point: CAMYLCASSLFGSPLHF against TRBV9 (CASSV) has a spurious CA at offset 0 and the real CASS at offset 5, and only the longest hit anywhere finds the second. An engine anchored at offset 0 scores the first and calls the record clean.

When it is not certain, flag it – do not change it. That is the whole policy, and it admits exactly four edits, in this order:

1. Re-call the V or J, when another allele of the locus explains 2 more residues contiguously from its own anchor (_RECALL_GAIN, anchor_depth()) and does not break an anchor the submission already holds. Confirming the submitted call is as much the job as replacing it.

2. Trim framework sitting past the conserved anchor of the allele that won. There is no N region beyond an anchor, so those residues are not junction.

3. Add back germline the submission was cut inside of, when the surviving run clears _MIN_HIT.

4. Substitute the first or last residue only, to the anchor the winning allele encodes, on a solid match.

Nothing else is ever rewritten. A residue that disagrees with germline inside the templated run is reported (mismatch) and left exactly as submitted, because a curation error and an allele IMGT does not record are indistinguishable from one junction – and rewriting them collapsed distinct records onto one string: CAAAETSYDKV{M,R,T,V}F on TRAJ50 are four junctions that all became CAAAETSYDKVIF. Four residues observed at one templated position is a suspect call or an unrecorded allele, not four independent typos.

A contradicted call changes the allele, not the sequence (guess_allele()), which is step 1 above. Measured over VDJdb’s corpus, 74.7 % of anchor-adjacent substitutions into an already-canonical junction were that case, against 1.1 % of untouched records.

Repair always targets a canonical junction. cdr3_repaired is only accepted when it opens with Cys104 and closes with Phe/Trp118 (_canonicalise); otherwise the submission is returned untouched. Legacy had no such rule. So good implies canonical, by construction.

The verdict is a set of flags per side (FLAGS), because one worst-wins label cannot say both “I trimmed a flank” and “I found a substitution I will not touch”. v_fix/j_fix keep VDJdb’s names so its cdr3fix JSON stays comparable key-for-key.

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, an edit further from the anchor than max_replace, 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_end_nt=-1, j_start_nt=-1, v_fix='FailedBadSegment', j_fix='FailedBadSegment', v_flags=('impossible', ), j_flags=('impossible', ), v_alts=(), j_alts=(), proposed=(), 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_end_nt (int)

• j_start_nt (int)

• v_fix (str)

• j_fix (str)

• v_flags (tuple[str, ...])

• j_flags (tuple[str, ...])

• v_alts (tuple[str, ...])

• j_alts (tuple[str, ...])

• proposed (tuple[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_end_nt: int = -1

j_start_nt: int = -1

v_fix: str = 'FailedBadSegment'

j_fix: str = 'FailedBadSegment'

v_flags: tuple[str, ...] = ('impossible',)

j_flags: tuple[str, ...] = ('impossible',)

v_alts: tuple[str, ...] = ()

Every allele that explains this end exactly as well as v_call / j_call does, the chosen one first. An amino-acid junction frequently cannot separate the alleles of a gene, so this is where that ambiguity is reported instead of being resolved by name order. Hand it to vdjtools.model.infer_nt_batch, which scores a LIST per row and can separate them.

j_alts: tuple[str, ...] = ()

proposed: tuple[str, ...] = ()

Sides whose call the SUBMISSION did not give and this proposed from the junction alone – ("V",), ("J",) or both. A consumer deciding whether to trust a call needs to know it was never curated, which is a different fact from allele (the submission named a different allele of the same locus).

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

Neither side is impossible, 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).

The anchor is guaranteed by _fix_side, which refuses any repair that would not land on the called allele’s own anchor residue – so this needs no separate canonicality conjunct, and must not use one: v_canonical/j_canonical are VDJdb’s literal C / [FW] tests, and TRBJ2-7*02’s anchor decodes V, which would make a perfectly canonical junction for that allele read as broken.

Read off the FLAGS, not off v_fix/j_fix. The fix-type names are VDJdb’s worst-wins vocabulary, and a declined repair projects onto NoFixNeeded there because that is what VDJdb’s scanner returned for it – its scan stopped at the anchor and never looked deeper in. arda does look deeper, so it knows the junction disagrees with the germline, and good has to mean what it says.

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='', loci=(), imgt_completeness='')

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)

• loci (tuple[str, ...])

• imgt_completeness (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 = ''

loci: tuple[str, ...] = ()

imgt_completeness: str = ''

arda.cdr3fix.load_anchors(organism)

{(segment, allele): Anchor}; Anchor.loci retains shared locus membership.

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.boundary_nt(junction, residues, germline_nt, side)

Where the germline stops in NUCLEOTIDES, from an amino-acid boundary.

residues is the amino-acid answer – Cdr3Markup.v_end (count of V-templated residues) or .j_start (index of the first J-templated residue). Returns the junction-nucleotide boundary and the probability of that value: for "V" a count, i.e. a half-open end, and for "J" the index of the first templated nucleotide. Both match the VEnd + 1 / JStart an aligner reports off the observed sequence.

The amino acid cannot say this and the nucleotide can. A germline run ends wherever the exonuclease stopped, which is not a codon boundary, so the last residue it touches is usually part germline and part N region – and an alignment on the protein can only round that to a whole residue. What the protein DOES fix is which nucleotides are admissible: 17 of the 20 residues have the same first base in every one of their codons, so on the V side the answer is frequently forced. GGA (Gly) against an observed Glu can only be GAA/GAG, both of which open with the germline’s G, so the germline demonstrably reaches one nucleotide further – whether it was templated or an insertion reproduced it, an aligner reads it as germline either way. Where the codons disagree, _BOUNDARY_W weighs them.

Measured against isalgo/airr_control’s human.trb.ntvj, on the 8,133 VDJdb human TRB junctions whose boundary every control observation agrees on. VDJdb’s v.end counts residues the V templates at least 2 nucleotides of, (nt + 1) // 3:

v.end	exact	within 1
k-mer scanner / 2.30.1	5,839 (71.8 %)	8,089 (99.5 %)
this	7,556 (92.9 %)	8,040 (98.9 %)

In nucleotides the same comparison is 28.96 % -> 80.39 % on V and 30.27 % -> 74.47 % on J. The two sides are limited by different things, and the numbers say which: the V-side residue count is right on 83.99 % of records and the extension on 95.71 % of those, so V is bounded by the protein alignment; the J-side count is right on 97.97 % and the extension on 76.02 %, because the nucleotide a J boundary turns on is the codon’s third, the one the genetic code leaves free.

j.start does not move: VDJdb defines it as the first FULLY J-templated residue, ceil(nt / 3), which is the same residue whether the germline reaches 1 or 2 nucleotides into the one before it. 97.97 % before and after, by construction.

Parameters:

• junction (str)

• residues (int)

• germline_nt (str)

• side (str)

Return type:

tuple[int, float]

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.guess_allele(cdr3, segment, locus, anchors, called='', min_gain=1, organism='')

The single best allele, or "" when nothing beats the submitted call. See guess_alleles() for the ones it ties with.

Parameters:

• cdr3 (str)

• segment (str)

• locus (str)

• anchors (dict)

• called (str)

• min_gain (int)

• organism (str)

Return type:

str

arda.cdr3fix.guess_alleles(cdr3, segment, locus, anchors, called='', min_gain=1, organism='')

Every allele that explains this junction end exactly as well as the best one does.

guess_allele() has to answer with one allele, so it resolves a tie by functionality and then by name – and a tie between two alleles is real information that the name order throws away. A junction’s first few residues often cannot separate the alleles of a locus at all, which is not a defect: it is what an amino-acid junction can say.

⚠ This is what the nucleotide stage should be given. vdjtools.model.infer_nt_batch takes a LIST of alleles per row and scores the reconstruction under each, so handing it the tie set lets codon plausibility and the model’s own usage settle what the junction could not – whereas handing it one allele makes arda’s name-order tie-break bind an answer it has no evidence for.

The chosen allele is first; the rest follow in candidate order (functional first, then by name).

Parameters:

• cdr3 (str)

• segment (str)

• locus (str)

• anchors (dict)

• called (str)

• min_gain (int)

• organism (str)

Return type:

tuple[str, …]

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

Note

arda.dpost was removed in 2.33.0. A posterior over the D gene is a recombination-model question, not a germline-reference one, so it lives in vdjtools as vdjtools.model.annotate_junctions (arda#144). The module was ported, not rewritten: its answer is identical field-for-field, and it gained the batch entry point arda#142 asked for. arda still ships the prior table it reads (database/vdj/<org>/d_prior.tsv) and still fits one (arda.scenarios.estimate()).

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, cell_from='', cell_regex=None, major_only=True, map_dd=True)

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.

  Never: 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)

• cell_from (str)

• cell_regex (str | None)

• major_only (bool)

• map_dd (bool)

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 "".

Never: 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"].

Never: 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.

Never: 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>, unmapped=<factory>, major_only=True)

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)

• unmapped (dict)

• major_only (bool)

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

unmapped: dict

major_only: bool = True

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, major_only=True, map_dd=True)

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)

• major_only (bool)

• map_dd (bool)

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', threads=0, 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, cell_from='', cell_regex=None, major_only=True, map_dd=True)

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).

• threads (int) – CPU budget for batched neighbour searches; 0 uses allocated CPUs. Use 1 when distributing samples across external workers.

• 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 | None) – 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).

• cell_from (str) – identifier dialect for arda.cell.make_parser() (migec, cellranger, prefix, auto). Turns on the AIRR umi_count column: distinct UMIs, per the spec, not distinct records. Empty disables it, and the column is then OMITTED rather than written as a column of 1s.

• cell_regex (str | None) – regex escape hatch for cell_from, as on arda.rnaseq.map.map().

• clonotype_key (str)

• call_level (str)

• flag_chimeras (bool)

• isotype (bool)

• report_path (str | Path | None)

• major_only (bool)

• map_dd (bool)

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, error_rate: 'float' = 0.0, error_rate_source: 'str' = 'given', error_rate_parents: '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)

• error_rate (float)

• error_rate_source (str)

• error_rate_parents (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.

Type:

Never

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

error_rate: float = 0.0

The substitution rate the run actually used, and – under –error-rate auto – what it was measured from. Never omitted under auto: a resolved parameter that is not written down makes the run unreproducible, and “it was measured” is not a number.

error_rate_source: str = 'given'

error_rate_parents: 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.

A sample delivered as SEVERAL FASTQ pairs – one per Illumina lane, or per in-house chunk – is the same shape arriving ready-made: whoever wrote the files already sharded it. So run() takes read groups rather than a pair, maps each one, concatenates in declared order, and calls finish() once. Nothing below Stage 1 knows read groups exist, which is deliberate: the clonotype key and the contig layout are where this project’s determinism guarantees live, and neither has to change for a sample to arrive in four files.

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, cell_from='', cell_regex=None, major_only=True, echo=None, map_dd=True)

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)

• cell_from (str)

• cell_regex (str | None)

• major_only (bool)

• map_dd (bool)

Returns:

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

Return type:

dict

arda.rnaseq.pipeline.run(pairs, out_dir, out_prefix, *, 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, cell_from='', cell_regex=None, major_only=True, echo=None, map_dd=True)

Single-node map -> assemble -> correct, over ONE sample’s read groups.

Parameters:

• pairs – [(r1, r2 | None), ...] for this sample, in the order they were declared. One pair is the ordinary case; several are the lanes or chunks one sample was delivered in, and they are concatenated after Stage 1, never before.

• out_dir (str | Path)

• out_prefix (str)

• 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)

• cell_from (str)

• cell_regex (str | None)

• major_only (bool)

• map_dd (bool)

Return type:

dict

Stages 2-3 stay global over the sample, exactly as they are for a single file. Nothing below Stage 1 knows read groups exist.

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, ec_mode='fast', min_junction_q=None, clonotype_key='full', call_level='allele', isotype=True, cell_from='', cell_regex=None, major_only=True, echo=None, map_dd=True)

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

Never: the Stage-2/3 knobs are reachable here for one reason – a sharded run must be able to reproduce what the mode command would have done. arda rnaseq defaults ec_mode to "rnaseq" and arda amplicon to "amplicon", while this function defaults to "fast" like arda correct does. A caller that forwards neither gets a different clonotype table from the same reads and nothing says so. arda.cluster.regime_flags() is what fills them in for the generated scripts.

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)

• ec_mode (str)

• min_junction_q (int | None)

• clonotype_key (str)

• call_level (str)

• isotype (bool)

• cell_from (str)

• cell_regex (str | None)

• major_only (bool)

• map_dd (bool)

Return type:

dict

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

Write <prefix>.stats.tsv and .stats.json from the run’s own artifacts.

Returns the row count.

Never: 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.

Never: 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', 'stats_json': '{prefix}.stats.json'}

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.

Never: 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).

Single cell

Lift a cell barcode out of a sequence identifier that already carries one.

arda has no barcode or UMI concept and does not need one: the tools that produce single-cell input put the barcode in the record NAME, and the name is what arda reads – map yields rec.id, the first whitespace-delimited token. (dnaio itself keeps the FASTQ comment on .comment; arda is what drops it. See project/design-singlecell.md.) This module is the whole of the single-cell front end – a pure parser, no I/O.

Three dialects ship, and a regex escape hatch for the rest:

cellranger	AAACCTGAGAAACCAT-1_contig_2 – 10x, with a gem-group suffix
migec	PBMC.AAACCTGCAGCCTGTT.CGTTTTTATC and its .c<k>/.<m> contig and split variants
prefix	<barcode>_<anything> / <barcode>.<anything>

Named dialects rather than regex alone, because each carries a validator a regex cannot express – see parse_migec below, which is a right-to-left walk and not a pattern at all.

class arda.cell.CellKey(cell=None, umi=None, sample=None, contig=1, molecule=1)

Bases: object

What a single-cell identifier resolves to.

contig numbers the fragments of ONE molecule; molecule numbers DIFFERENT molecules that shared a barcode. Both default to 1, which is what an identifier carrying neither means.

Parameters:

• cell (str | None)

• umi (str | None)

• sample (str | None)

• contig (int)

• molecule (int)

cell: str | None

umi: str | None

sample: str | None

contig: int

molecule: int

arda.cell.parse(sequence_id, dialect)

Parameters:

• sequence_id (str)

• dialect (str)

Return type:

CellKey | None

arda.cell.sniff(sequence_ids, *, min_match=0.95)

The dialect that parses at least min_match of sequence_ids to a consistent cell.

Never: an alphabet check alone fabricates cells. A bulk sample named TCGA, CAG or any donor code over those four letters parses as a one-cell library, and nothing downstream flags it. So a dialect is only accepted when every parsed cell has the SAME LENGTH – a real 10x cell field is exactly 16 nt in every record, and a sample id that happens to look like DNA is not.

The caller is responsible for handing this a REPRESENTATIVE sample. migec writes in barcode order, so the first N identifiers of a file share their leading bases and are not one.

Parameters:

min_match (float)

Return type:

str | None

arda.cell.make_parser(dialect='auto', regex=None, *, sample=(), with_sample=False)

A sequence_id -> cell or None callable.

regex wins when given and must define a cell group. auto runs sniff() over sample and refuses rather than guessing when nothing fits.

with_sample returns (sample, cell) instead, so a caller can check that one file holds one sample. A 10x barcode comes from a fixed 737,280-entry whitelist, so the same cell across two samples is guaranteed, not unlucky – keying on the barcode alone merges two cells and nothing downstream can tell. Only the migec dialect carries a sample; the others return None for it.

Parameters:

• dialect (str)

• regex (str | None)

• with_sample (bool)

Reference-free per-cell contig assembly.

A droplet V(D)J library gives one barcode per cell and, under it, many UMI-tagged molecules. An upstream tool (migec, fgbio, UMI-tools) collapses each molecule’s reads into one consensus; what arrives here is therefore a cell worth of short, error-corrected fragments. Reads of one molecule are co-terminal in this chemistry, so a molecule’s consensus is a pile at one position and covers one window of the transcript however deep it is – but different molecules start at different positions, so the cell’s molecules tile the whole transcript. Measured on 10x’s sc5p_v2_hs_PBMC_1k VDJ-T: every one of Cell Ranger’s 943 filtered contigs is >= 99.9% covered by the 25-mers of its own cell’s molecules. The full-length contig is already in the data; assembling it is what this module does.

The method is overlap-layout with no germline reference anywhere: index every 25-mer of a cell’s molecules, propose an offset wherever two molecules share one, verify the implied overlap, and union-find the survivors into components. Each component’s weighted column consensus is a contig.

Two things are load-bearing and both were measured by removing them (CDR3 nt recovered exactly inside a contig, 244 Cell Ranger chains over the first 120 cells):

trimmed + depth-weighted	0.988
trimmed, every molecule weighted 1	0.967
untrimmed + depth-weighted	0.898
untrimmed, every molecule weighted 1	0.873

• Never: verify the overlap before joining, and trim the read-through adapter first. Every molecule that ran off the end of its insert carries the same Illumina adapter, so an adapter suffix is a shared k-mer and a 100%-identical overlap: verification alone cannot refuse it. Before verification existed, one component swallowed a cell’s TRA and TRB alike. With it, the damage is quieter and still 9 points – the largest component of a cell grows from 70 to 78 molecules (median over 40 cells) and the adapter columns corrupt the consensus where it does merge. Both failures have the same cure, and neither substitutes for the other.

• Never: weight the consensus by the molecule’s read depth. A one-read molecule is one read; giving it the same column vote as an 800-read consensus is what the 0.967 row costs.

Assembly is per cell and cells are independent, so nothing here scales with the library. Note: it is single-threaded Python at ~54 ms per cell (479 cells in 26 s). A process pool over cells is the obvious upgrade and was tried and reverted – see the comment in assemble_cells(); on macOS spawn re-imports the caller’s __main__ in every worker and forkserver will not start, so a pool here crashes library callers to save minutes on a 10k-cell library. The upgrade path is the assembler in C++, not a pool.

class arda.singlecell.Contig(cell, index, sequence, molecules, reads)

Bases: object

One overlap component of a cell, and what supports it.

Parameters:

• cell (str)

• index (int)

• sequence (str)

• molecules (int)

• reads (int)

cell: str

index: int

sequence: str

molecules: int

reads: int

property contig_id: str

<cell>_contig_<n>, which arda.cell.parse() reads back as the cellranger dialect – so an annotation run over these contigs recovers cell_id for free.

class arda.singlecell.CellAssemblyReport(cells: 'int' = 0, molecules_in: 'int' = 0, molecules_placed: 'int' = 0, contigs: 'int' = 0, contigs_dropped_small: 'int' = 0, contig_lengths: 'list[int]' = <factory>, molecules_by_cell: 'dict[str, int]'=<factory>)

Bases: object

Parameters:

• cells (int)

• molecules_in (int)

• molecules_placed (int)

• contigs (int)

• contigs_dropped_small (int)

• contig_lengths (list[int])

• molecules_by_cell (dict[str, int])

cells: int = 0

molecules_in: int = 0

molecules_placed: int = 0

contigs: int = 0

contigs_dropped_small: int = 0

contig_lengths: list[int]

molecules_by_cell: dict[str, int]

cell -> molecules READ for it, before assembly. Never derive this by summing contig membership: a molecule that joined no contig is missing from that sum, and one shared between two phased haplotypes is counted twice – so a doublet reads as the deepest cell in the library and the knee plot’s y axis stops being unique molecules.

as_dict()

Return type:

dict

arda.singlecell.assemble_cell(cell, sequences, weights=None, *, k=25, min_overlap=30, min_identity=0.95, min_molecules=2, split=True, min_minor_fraction=0.2, min_minor_molecules=2, min_informative_columns=2)

Assemble one cell’s molecule consensuses into contigs, largest component first.

weights is each molecule’s read depth (cD:i: in a migec consensus); None weighs every molecule 1, which costs about two points of junction recovery – see the module table.

Never: a component must be phased before it is consensed, or a doublet is invisible. Two chains of the SAME locus in one cell share their constant region and most of their V, so the overlap layout puts them in one component and the column consensus averages their junctions into a sequence that is neither. Measured on two TRB receptors tiled into one synthetic cell: one 918 nt contig, no junction called at all, both true junctions lost. With the split on, both come back. This is the same failure a de Bruijn graph has on shared germline, and it is why split defaults to on: the layout finds the component, the phasing finds the chain.

Never: the split gate is on the minor base’s WEIGHT SHARE and its MOLECULE count, both. A share alone promotes one deep molecule’s error to an allele; a count alone promotes a column that a hundred molecules cover and two disagree on. min_informative_columns requires the disagreement to co-segregate across positions rather than stand on one.

Parameters:

• cell (str)

• sequences (list[str])

• weights (list[int] | None)

• k (int)

• min_overlap (int)

• min_identity (float)

• min_molecules (int)

• split (bool)

• min_minor_fraction (float)

• min_minor_molecules (int)

• min_informative_columns (int)

Return type:

list[Contig]

arda.singlecell.assemble_cells(reads, output, *, cell_from='auto', cell_regex=None, k=25, min_overlap=30, min_identity=0.95, min_molecules=2, min_length=0, trim_adapter=True, split=True, min_minor_fraction=0.2, min_minor_molecules=2, min_informative_columns=2, cells=None)

Assemble every cell in a UMI-consensus FASTQ and write the contigs as FASTA.

Parameters:

• reads (str | Path) – a per-molecule consensus FASTQ/FASTA whose record names carry the cell barcode – what migec assemble writes, or anything arda.cell can parse.

• output (str | Path) – contigs FASTA, named <cell>_contig_<n>.

• cell_from (str) – a arda.cell dialect, "auto" to sniff, or "regex" with cell_regex.

• min_molecules (int) – molecules a component needs before it is emitted. Never: 1 is not a useful setting – on the 10x library, an extra chain supported by exactly one molecule is contamination 96% of the time (see pair_chains()).

• cells (set[str] | None) – restrict to these barcodes – the called cells. Never: cell CALLING is the upstream tool’s job (migec’s refine writes <sample>.cells.tsv); assembling every empty droplet costs the run and calls chains in ambient RNA.

• cell_regex (str | None)

• k (int)

• min_overlap (int)

• min_identity (float)

• min_length (int)

• trim_adapter (bool)

• split (bool)

• min_minor_fraction (float)

• min_minor_molecules (int)

• min_informative_columns (int)

Returns:

The contigs, and a CellAssemblyReport.

Return type:

tuple[list[Contig], CellAssemblyReport]

arda.singlecell.read_molecules(path, *, cell_from='auto', cell_regex=None, sniff_sample=10000)

Group a consensus FASTQ/FASTA by cell barcode, keeping each record’s read depth.

Never: the dialect is sniffed over a RESERVOIR sample of the names, never the first N. migec assemble writes in barcode order, so the head of the file shares its leading bases and any statistic taken from it describes one corner of the barcode space.

Never: TWO SAMPLES IN ONE FILE ARE REFUSED. The grouping key is the cell barcode alone, and a droplet barcode comes from a fixed whitelist – concatenating two samples’ consensus FASTQs merges every barcode they share, and the result reads as an ordinary doublet rather than as an error. migec writes one consensus file per sample; run this once per file.

Parameters:

• path (str | Path)

• cell_from (str)

• cell_regex (str | None)

• sniff_sample (int)

Return type:

dict[str, list[tuple[str, int]]]

arda.singlecell.trim(seq, *, adapter='AGATCGGAAGAGC', min_stub=8)

Cut a read-through adapter and any homopolymer tail off a molecule consensus.

Both orientations are tried, and a partial stub counts: a consensus that ended eight bases into the adapter still shares those eight bases with every other such consensus.

Parameters:

• seq (str)

• adapter (str)

• min_stub (int)

Return type:

str

arda.singlecell.pair_chains(rows, *, min_extra_molecules=3, min_extra_fraction=0.2, require_productive=True)

Turn annotated contigs into one row per (cell, locus, junction), with a doublet flag.

rows need cell_id, locus, junction_aa and molecules; productive, reads, sequence and the gene calls are carried when present. Within a cell and locus the chains are ranked by supporting molecules. Rank 1 is primary. A lower-ranked chain is believed only if it is productive and clears BOTH support gates, and then a second HEAVY chain is doublet_candidate while a second LIGHT chain is secondary – allelic inclusion at TRA and IGK/IGL is a known, real population, not a doublet. Everything else is extra, kept in the table with the columns that say why.

Never: an extra chain must be PRODUCTIVE, and that is the discriminator, not the count. Measured on sc5p_v2_hs_PBMC_1k against Cell Ranger’s own per-cell calls: of the extra chains it agrees with, 60/60 are productive; of those it does not, 20/128 are. Molecule support separates the same two groups far less sharply (median 10 against 4).

Never: the threshold is on MOLECULES, not reads, and 1 is not a threshold. Sweeping the extra chain’s molecule support against Cell Ranger’s calls, before the productivity gate:

locus	min	precision	recall
TRA	1	0.360	1.000
TRA	2	0.946	0.972
TRA	3	1.000	0.917
TRB	1	0.034	1.000
TRB	3	0.727	1.000
TRB	5	0.889	1.000

An extra chain supported by exactly one molecule is contamination in 96-97% of cases and a real second chain in none of the TRB cases; the distribution is bimodal, not a tail. That is why assemble_cells() defaults min_molecules to 2 and never to 1.

Parameters:

• rows (list[dict])

• min_extra_molecules (int)

• min_extra_fraction (float)

• require_productive (bool)

Return type:

list[dict]

arda.singlecell.cell_summary(pairs, molecules_per_cell=None)

One row per cell: what chains it holds and what that makes it.

pairs is pair_chains() output. The row carries the top two HEAVY and top two LIGHT chains by molecule support, which is what the doublet scatter is drawn from – a plot of heavy_molecules_2 against heavy_molecules_1 separates the ambient cloud along the bottom from the multiplets off the axis.

cell_status is one of:

value	meaning
paired	one heavy and at least one light chain
heavy_only	a heavy chain and no light one
light_only	a light chain and no heavy one
doublet_candidate	two heavy chains both clearing the support threshold
no_chain	no chain was called at all

Parameters:

• pairs (list[dict])

• molecules_per_cell (dict[str, int] | None)

Return type:

list[dict]

arda.singlecell.cell_rank(molecules_per_cell)

The knee curve: cells ranked by how many molecules they carry.

Never: the y axis is MOLECULES, never reads. One over-amplified molecule otherwise puts an empty droplet high on the curve, which is the artefact the plot exists to show.

Parameters:

molecules_per_cell (dict[str, int])

Return type:

list[dict]

arda.singlecell.threshold_sweep(pairs, reference, *, thresholds=(1, 2, 3, 5, 8, 12))

Precision and recall of the extra-chain filter, per locus, against a reference call set.

reference is {(cell_id, locus, junction_aa)} from a tool you are willing to score against – Cell Ranger’s filtered_contig_annotations.csv, a hashtag demultiplex, a simulation’s truth. Only NON-primary chains are swept: the top chain of a cell is kept at every threshold, so including it would put a constant in both columns and flatten the curve. Each locus is swept twice, with the productivity gate off and on, so the two levers are separable rather than confounded.

This is the table that picks min_extra_molecules and require_productive; both are emitted rather than baked in, because the right setting is a property of the library’s ambient load and not of the code.

Parameters:

• pairs (list[dict])

• reference (set[tuple[str, str, str]])

• thresholds (tuple[int, ...])

Return type:

list[dict]

arda.singlecell.HEAVY = {'IGH', 'TRB', 'TRD'}

Which locus of a class carries the rearrangement that is unique per cell. Two of these in one cell is the doublet signal; two light chains are not – allelic inclusion at TRA and IGK/IGL is a known, real population.

arda.singlecell.read_called_cells(path)

The called barcodes, from a bare one-per-line list or from migec’s <sample>.cells.tsv.

A header naming a cell (or cell_id/barcode) column switches to table mode, and a called column is then respected: migec writes every barcode it saw with a boolean, so taking column 1 unfiltered would restore the empty droplets this option exists to remove. A -1 gem-group suffix is stripped.

Parameters:

path (str | Path)

Return type:

set[str]

arda.singlecell.read_reference_calls(path)

{(cell_id, locus, junction_aa)} from a comparator’s own call table.

Accepts Cell Ranger’s filtered_contig_annotations.csv (barcode/chain/cdr3, gem-group suffix stripped, is_cell respected when present) and any TSV/CSV carrying cell_id/locus/junction_aa. This is only ever a scoring input – assemble_cells() never reads a reference of any kind.

Parameters:

path (str | Path)

Return type:

set[tuple[str, str, str]]

arda.singlecell.read_records(path)

Yield (name, sequence, depth) from a consensus FASTQ or FASTA.

Never: the read depth is in the FASTQ comment (cD:i:), which the shared reader drops with the rest of the header. It is read here rather than by widening read_sequences(), whose callers are the bulk hot path. A record with no cD:i: weighs 1.

Parameters:

path (str | Path)

arda.singlecell.run(reads, prefix, *, organism='human', cell_from='auto', cell_regex=None, k=25, min_overlap=30, min_identity=0.95, min_molecules=2, min_length=100, trim_adapter=True, split=True, min_minor_fraction=0.2, min_minor_molecules=2, min_informative_columns=2, min_extra_molecules=3, min_extra_fraction=0.2, require_productive=True, cells_file=None, reference=None, threads=0, map_d=True, plot=None, gnuplot=None, major_only=True, map_dd=True)

Assemble, annotate, pair and diagnose a droplet V(D)J library. Returns the run report.

Writes, all under prefix:

suffix	what is in it
.contigs.fasta	the assembled per-cell contigs
.contigs.airr.tsv	their AIRR annotation, cell_id included
.chains.tsv	one row per (cell, locus, junction), ranked, with a status
.cells.tsv	one row per cell, with the doublet scatter’s two columns
.cell_rank.tsv	the knee curve, molecules per cell
.contig_lengths.tsv	the contig length histogram
.sweep.tsv	precision/recall of the extra-chain filter (needs reference)
.partition.tsv	clustering agreement of the clonotypes; see arda.partition
.arda.json	the run report

Parameters:

• reads (str | Path)

• prefix (str | Path)

• organism (str)

• cell_from (str)

• cell_regex (str | None)

• k (int)

• min_overlap (int)

• min_identity (float)

• min_molecules (int)

• min_length (int)

• trim_adapter (bool)

• split (bool)

• min_minor_fraction (float)

• min_minor_molecules (int)

• min_informative_columns (int)

• min_extra_molecules (int)

• min_extra_fraction (float)

• require_productive (bool)

• cells_file (str | Path | None)

• reference (str | Path | None)

• threads (int)

• map_d (bool)

• plot (str | None)

• gnuplot (str | None)

• major_only (bool)

• map_dd (bool)

Return type:

dict

arda.singlecell.find_knee(sizes, *, min_ratio_to_mean=10.0)

Where a descending molecules-per-barcode curve breaks. (rank, molecules), 1-based.

sizes must be sorted descending. Returns (0, 0) when the curve has no knee worth reporting, which is a real answer and not a failure.

This is Kneedle (Satopaa 2011) taken at the global maximum of its difference curve. On a decreasing curve the distance to the chord joining the endpoints and Kneedle’s normalised |y_n - (1 - x_n)| are the same quantity up to a constant; what differs is the selection rule.

Never: Kneedle’s published rule – walk the LOCAL maxima, stop at the first that then falls by a sensitivity step – needs the smoothing spline the paper specifies, and without it the answer is degenerate. Measured on a 136,032-barcode 10x library: unsmoothed, the difference curve carries 287 local maxima and the walk stops at rank 15 (1,057 molecules), which would call fifteen cells. Smoothed over a 0.01 window in normalised log-rank it gives rank 358 (314 molecules) – the global maximum, rank 376 (308 molecules), to within 5% – and over-smoothing drifts it off again (rank 180 at a 0.1 window). The global maximum reaches the same place with no smoothing parameter to tune.

Never: the guard is not decoration. A global maximum always exists, so an ambient-only library returns a rank too, and a “knee” there is a threshold that calls every barcode a cell. The floor is min_ratio_to_mean times the mean molecules per barcode – one order of magnitude, which is the unit a log-log curve is read in. The mean alone is too weak and was tried: an ambient-only library of 1-3 molecules per barcode is a step curve whose corner sits at 3 against a mean of 2.0 and clears a mean-only test at 1.5x. On the real library the knee is 308 against a mean of 3.44, a factor of 89. Two orders apart, so the boundary is not tuned.

Parameters:

• sizes (list[int])

• min_ratio_to_mean (float)

Return type:

tuple[int, int]

External clustering validation: does our partition of the cells match a reference’s?

Recall and concordance answer “did we find this chain”. They do not answer “did we put the same cells in the same clone”, which is the question a repertoire analysis actually runs on, and which a per-chain metric cannot see: splitting one clone in two and merging two clones into one both leave every individual chain call correct.

The scores implemented here are the homogeneity-parsimony trade-off proposed in arXiv:2607.20799, whose reference implementation is qimmuno/clustereval. Three families, each a pair of one score that punishes merging and one that punishes splitting:

family	punishes merging	punishes splitting	entropy
information-theoretic	homogeneity	parsimony	Shannon
set-matching	normalized purity	normalized inverse purity	min-entropy
pair-counting	pair specificity	pair sensitivity	collision

Note: the formulas are implemented here rather than imported, so arda gains neither a dependency nor an optional-import branch for sixty lines of counting. They are transcribed from the reference implementation and tests.unit.test_partition checks them against its published worked values; the trade-off is that a change upstream will not follow automatically, so the version transcribed is recorded here: clustereval as of commit 9f6daba.

Never: parsimony is not completeness. Completeness normalises H(K|C) by H(K), the entropy of the clustering being judged, so a method can raise its own score by making its clustering more uniform. Parsimony normalises by log N - H(C), which depends only on the reference – so the denominator is the same for every method being compared, which is the whole point of comparing them.

arda.partition.homogeneity_score(labels_true, labels_pred)

1 when no cluster mixes two reference classes. Punishes MERGING.

Return type:

float

arda.partition.parsimony_score(labels_true, labels_pred)

1 when no reference class is split across clusters. Punishes SPLITTING.

Return type:

float

arda.partition.purity_score(labels_true, labels_pred)

Share of items in the majority reference class of their own cluster.

Return type:

float

arda.partition.normalized_purity_score(labels_true, labels_pred)

Purity rescaled so that the one-cluster clustering scores 0 rather than a class share.

Return type:

float

arda.partition.inverse_purity_score(labels_true, labels_pred)

Share of items in the majority cluster of their own reference class.

Return type:

float

arda.partition.normalized_inverse_purity_score(labels_true, labels_pred)

Inverse purity rescaled so that the singleton clustering scores 0.

Return type:

float

arda.partition.pair_specificity_score(labels_true, labels_pred)

Share of reference-different pairs we also kept apart. Punishes MERGING.

Return type:

float

arda.partition.pair_sensitivity_score(labels_true, labels_pred)

Share of reference-same pairs we also kept together. Punishes SPLITTING.

Return type:

float

arda.partition.compare_partitions(labels_true, labels_pred)

All six scores plus the sizes, as one flat dict ready for a TSV row.

Never: report BOTH scores of a family or neither. Either alone is trivially maximised – homogeneity by putting every item in its own cluster, parsimony by putting them all in one – which is exactly what the trade-off framing exists to prevent.

Never: read classes_singleton before reading the scores. A reference partition that is almost all singletons carries almost no clustering to agree about – an unexpanded T-cell repertoire is 476 clonotypes over 479 cells – and every score is then decided by a handful of items. It is reported beside them for that reason, not as decoration.

Return type:

dict[str, float]

Single-cell QC figures, drawn with gnuplot from the tables arda cells already wrote.

Nothing here computes anything. Every panel is a gnuplot script over a TSV that is written whether or not gnuplot exists, which is what lets a figure be redrawn months later by someone who no longer has the FASTQ. It also means arda gains no plotting dependency: the .gp scripts are written either way and gnuplot runs only if it is on PATH.

Colours are ColorBrewer Dark2 (qualitative, colour-blind safe), never hand-picked. Transparent background and one mid-grey ink, so a single SVG serves a light page, a dark page and print.

arda.scplot.draw(prefix, *, fmt='svg', gnuplot=None)

Write a .gp script per available table and render it if gnuplot is present.

gnuplot="" writes the scripts and renders nothing, which is also how the tests exercise the no-gnuplot path. Returns the figures that were actually rendered.

Parameters:

• prefix (str | Path)

• fmt (str)

• gnuplot (str | None)

Return type:

list[Path]

arda.scplot.PANELS: dict[str, tuple[str, str, str]] = {'.cell_rank.tsv': ('cell_rank', 'Molecules per cell, ranked', '\nset title "Cell barcode rank plot"\nset xlabel "cell barcode rank"\nset ylabel "molecules (unique UMIs)"\nset logscale xy\nset key inside bottom left box lw 0.5 lc rgb "#808080" textcolor rgb "#808080"\n# The knee, drawn from the same table rather than recomputed: the rank whose molecule count is\n# the largest in the file at or above the reported knee. `arda cells` writes the knee into its\n# JSON report and marks the row here, so the line and the number cannot drift apart.\nplot "{data}" using (column("rank")):(column("molecules")) \\\n        with lines lw 2 lc rgb "#1b9e77" title "barcodes", \\\n     "{data}" using (strcol("knee") eq "1" ? column("rank") : 1/0):(column("molecules")) \\\n        with points pt 7 ps 1.4 lc rgb "#d95f02" title "knee"\n'), '.cells.tsv': ('doublet_scatter', 'Second heavy chain against first', '\nset title "Doublet scatter: cells carrying two heavy chains"\nset xlabel "molecules on the cell\'s top heavy chain"\nset ylabel "molecules on its second heavy chain"\nset logscale xy\nset key inside top left box lw 0.5 lc rgb "#808080" textcolor rgb "#808080"\n# Never: only cells that HAVE a second heavy chain are drawn. A cell with one chain has y = 0,\n# which a log axis cannot show, and substituting a floor value would invent the datum the plot\n# is being read for. The one-chain cells are the knee plot\'s population, not this one.\nh1 = "heavy_molecules_1"; h2 = "heavy_molecules_2"\nplot "{data}" using (column(h2) > 0 ? column(h1) : 1/0):(column(h2)) \\\n        with points pt 7 ps 0.6 lc rgb "#666666" title "two heavy chains", \\\n     "{data}" using (strcol("cell_status") eq "doublet_candidate" ? column(h1) : 1/0):(column(h2)) \\\n        with points pt 7 ps 1.0 lc rgb "#d95f02" title "doublet candidate"\n'), '.chains.tsv': ('chain_support', 'Molecules per chain, by rank', '\nset title "Molecules supporting a cell\'s chains, by rank within the cell"\nset xlabel "molecules on the chain (powers of two)"\nset ylabel "chains"\nset logscale xy\nset key inside top right box lw 0.5 lc rgb "#808080" textcolor rgb "#808080"\nbin(x) = x < 1 ? 1 : 2**(floor(log(x)/log(2)))\nplot "{data}" using (column("chain_rank") == 1 ? bin(column("molecules")) : 1/0):(1) \\\n        smooth frequency with steps lw 2 lc rgb "#1b9e77" title "rank 1", \\\n     "{data}" using (column("chain_rank") > 1 ? bin(column("molecules")) : 1/0):(1) \\\n        smooth frequency with steps lw 2 lc rgb "#e7298a" title "rank 2+"\n'), '.contig_lengths.tsv': ('contig_lengths', 'Assembled contig lengths', '\nset title "Assembled contig length"\nset xlabel "contig length (nt, 25 nt bins)"\nset ylabel "contigs"\nset logscale y\nunset key\nset boxwidth 22\nplot "{data}" using (column("length_bin")):(column("contigs")) \\\n        with boxes lc rgb "#7570b3"\n'), '.sweep.tsv': ('filter_sweep', 'Extra-chain filter, precision against recall', '\nset title "Extra-chain filter: precision and recall against the reference calls"\nset xlabel "minimum molecules on the extra chain"\nset ylabel "rate"\nset yrange [0:1.05]\nset key inside bottom left box lw 0.5 lc rgb "#808080" textcolor rgb "#808080"\nprod(v) = strcol("require_productive") eq v ? column("min_molecules") : 1/0\nplot "{data}" using (prod("F")):(column("precision")) \\\n        with points pt 6 ps 0.9 lc rgb "#1b9e77" title "precision, any chain", \\\n     "{data}" using (prod("T")):(column("precision")) \\\n        with points pt 7 ps 0.9 lc rgb "#1b9e77" title "precision, productive only", \\\n     "{data}" using (prod("F")):(column("recall")) \\\n        with points pt 4 ps 0.9 lc rgb "#d95f02" title "recall, any chain", \\\n     "{data}" using (prod("T")):(column("recall")) \\\n        with points pt 5 ps 0.9 lc rgb "#d95f02" title "recall, productive only"\n')}

suffix -> (stem, title, body). The body is gnuplot with {data} for the TSV path.

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

junction_aa_len / read_len / clone_size / isotype

the distributions, all keyed locus:bucket – junction length in residues, aligned read length in 10-nt buckets, clonotype size in powers of two, constant-region class

Never: 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.

Never: 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 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. Never: arda genotype is a DIFFERENT object and this one does not feed it – that call is a likelihood ratio between diploid genotypes over clonotypes, restricted to alleles the reference already catalogues, and it never reads this shortlist.

arda.stats.collect(*, airr=None, clones=None, report=None, cells=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.

• cells (str | Path | None) – an arda cells output PREFIX (not a file) -> its report plus .chains.tsv, as the same sample / chain / *_gene / junction_aa_len scopes a bulk run produces, so one batch table holds both kinds of sample.

• 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.write_stats_json(rows, output, *, arda_version='')

The same rows as nested, TYPED JSON: {scope: {key: {metric: value}}}.

The TSV is what an operator greps and what a shell pipeline joins; this is what a program reads, and what arda qc report inlines into its HTML. Both are written from the SAME rows list rather than from two passes over the data, so they cannot disagree – the cost is re-parsing a number that was just formatted, which is nothing against the frames it came from.

Values are typed here and stringified there for the same reason the TSV drops blank rows: a consumer must be able to tell “no IGH read spanned both anchors” from “the minimum was 0”.

Parameters:

• rows (list[tuple])

• output (str | Path)

• arda_version (str)

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.

Never: 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.

Never: 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.arda_binary()

The arda executable a generated script should call, resolved at RENDER time.

Never: A generated script must name the arda it was generated BY, not whatever arda happens to resolve to on a compute node. Measured on a real cluster: the map array succeeded across every read group and the reduce array then died with No such option: --ec-mode, because the workers’ PATH held 2.19.0 – the map flags happened to exist in that version, the reduce flag did not. A sheet’s worth of Stage-1 AIRR, no clonotypes, and an error message that names a flag rather than a version. ARDA_MMSEQS is already pinned for exactly this reason; arda itself was the one tool in the script left to PATH.

Absolute but NOT symlink-resolved: an HPC module system points a stable name at the current build on purpose, and following that link would pin a path the site means to move.

Return type:

str

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

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)

• seqtype (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.plan(samples, work_dir)

Write the two manifests a multi-sample run is executed from; make the shard dirs.

arda’s unit of parallel work is the read group, not the sample. Stage 1 is per read and shards perfectly, so every read group of every sample is an independent job – a sheet of 5 samples x 4 lanes is 20 concurrent tasks, not 5. Stages 2-3 are global per sample (a clone split across read groups would be counted once per group, and a contig tiling across them would never be built), so they run once per sample, in arda.rnaseq.pipeline.reduce().

Any scheduler can therefore drive arda from these two files and nothing else: one arda map per row of readgroups.tsv, then one arda cluster reduce per row of samples.tsv.

Never: the part names are shard_%05d.airr.tsv and the index is the read group’s position in the sheet, not its filename. reduce merges the parts in NAME order and that order is read order, so both the 5-digit padding (shard_00010 must follow shard_00002) and declared order are load-bearing – together they are what makes the result byte-identical to the same reads in one file.

Parameters:

• samples – arda.samples.Sample objects, in the order they were declared.

• work_dir (str | Path) – scratch root; each sample gets work_dir/<id>/.

Returns:

(readgroups.tsv, samples.tsv).

Return type:

tuple[Path, Path]

arda.cluster.regime_flags(regime, *, organism='human', threads=8, kmer=12, min_score=75.0, reconstruct=False, map_d=True, assemble=True, complete_only=True, ec_mode=None, cell_from='', map_dd=True)

The arda map and arda cluster reduce flags one regime implies. Returns both.

Never: THE TWO HALVES MUST AGREE, and only this function knows how. arda rnaseq runs both stages in one call, so it can wire Stage 1 to Stage 2 itself – it turns on --junction-quality whenever the denoising preset will look for that column, because “asking for the gate has to imply producing its input”. A path that invokes map and reduce as separate jobs has no such link: the preset defaults differ between the two commands, the quality column is simply absent, and the gate then does nothing at all. Exit 0, plausible clonotypes, a quality-directed rescue that silently never ran.

Never: the two speed levers do NOT compose and each is a loss in the other’s regime, so the regime is named here ONCE, by library type, rather than spelled as loose flags per job.

Never: --cell-from is a BOTH-HALVES flag, for the same reason --junction-quality is. Stage 1 needs it to write cell_id; Stage 2 needs it to write umi_count, which it derives from sequence_id itself and cannot get from Stage 1’s output. Passed to one half only, a sharded run silently drops whichever column that half owns – exit 0, plausible clonotypes, a missing column.

Parameters:

• regime (str)

• organism (str)

• threads (int)

• kmer (int)

• min_score (float)

• reconstruct (bool)

• map_d (bool)

• assemble (bool)

• complete_only (bool)

• ec_mode (str | None)

• cell_from (str)

• map_dd (bool)

Return type:

tuple[str, str]

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, arda_bin=None, map_dd=True)

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)

• arda_bin (str | None)

• map_dd (bool)

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', regime='rnaseq', arda_mmseqs=None, arda_bin=None, map_dd=True)

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)

• regime (str)

• arda_mmseqs (str | None)

• arda_bin (str | None)

• map_dd (bool)

Return type:

str

arda.cluster.render_samples_submit_script(readgroups, samples_manifest, *, n_read_groups, n_samples, out_dir='.', log_dir=None, organism='human', threads=8, kmer=12, min_score=75.0, reconstruct=False, assemble=True, complete_only=True, map_d=True, regime='rnaseq', partition=None, time='04:00:00', mem='8G', reduce_time='08:00:00', reduce_mem='16G', arda_mmseqs=None, arda_bin=None, map_dd=True)

Render a submit.sh for a whole sample sheet: array-map per READ GROUP, then reduce.

A third sibling of render_submit_script() and render_rnaseq_submit_script(), not a generalisation of either – see the note on the second one about a parameter that decides which of two unrelated pipelines you get. What differs here is that there is no split step: a sample delivered as several files was already sharded by whoever wrote it, so the read groups are the shards and the array indexes them straight out of the manifest.

Two arrays, not one job per sample. The map array is sized by read groups across every sample, so the whole sheet is in flight at once. The reduce array is sized by samples and gated on the entire map array with afterok – deliberately not aftercorr, whose index-for-index pairing is exactly wrong here, since one sample consumes many map tasks.

Parameters:

• readgroups (str | Path)

• samples_manifest (str | Path)

• n_read_groups (int)

• n_samples (int)

• out_dir (str | Path)

• log_dir (str | Path | None)

• organism (str)

• threads (int)

• kmer (int)

• min_score (float)

• reconstruct (bool)

• assemble (bool)

• complete_only (bool)

• map_d (bool)

• regime (str)

• partition (str | None)

• time (str)

• mem (str)

• reduce_time (str)

• reduce_mem (str)

• arda_mmseqs (str | None)

• arda_bin (str | None)

• map_dd (bool)

Return type:

str

arda.cluster.READGROUP_COLUMNS = ('sample', 'read_group', 'r1', 'r2', 'part_airr', 'part_report')

Work-unit manifest written by plan(). One row per READ GROUP – the unit a scheduler allocates a worker to.

arda.cluster.SAMPLE_COLUMNS = ('sample', 'shard_dir', 'out_prefix')

The reduce manifest written beside it. One row per SAMPLE.

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. NCBI’s IMGT conversion ungaps them and rewrites headers to allele names (what makeblastdb wants); the packaged Python implementation follows edit_imgt_file.pl without requiring Perl.

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 using NCBI’s IMGT conversion 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 NCBI’s IMGT conversion; 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(*, fetch_missing=True)

Pin one complete runtime for execution and metadata; optional probes stay offline.

Parameters:

fetch_missing (bool)

Return type:

Path

arda.igblast.igblast_version()

Version of this process’s pinned root, independent of later refreshes.

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)

Apply NCBI’s edit_imgt_file.pl conversion without a Perl prerequisite.

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.

Never: 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 root, pinned to one complete generation per process.

Source checkout → the committed database/. PyPI install → a verified generation within <cache>/database, fetched once 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 versioned data release (independent of the package version) and downloaded once into the per-user cache on first use. Readers pin a verified immutable generation. 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.installed_database(dest, *, version=None)

Resolve one complete immutable generation, rejecting missing/corrupt cache files.

Parameters:

• dest (Path)

• version (str | None)

Return type:

Path | None

arda._database_fetch.legacy_database(dest, *, version=None)

Accept a flat legacy cache only when all canonical data hashes match.

Parameters:

• dest (Path)

• version (str | None)

Return type:

Path | None

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

Fetch a complete reference and return its immutable root (containing vdj/).

Publication replaces a small pointer atomically. Existing readers keep their old root; failed or interrupted refreshes cannot remove it. A warm valid cache needs no network.

Parameters:

• dest (Path)

• force (bool)

• version (str | None)

Return type:

Path

Preliminary segment evidence

arda.alignment_candidates(records, *, organism='human', threads=0, max_seqs=300, evalue=0.001, sensitivity=7.0, batch_size=10000, major_only=True)

Return independent nucleotide V/J/D alignment candidates for each input.

Major-only search is the default; major_only=False searches all alleles. The returned major_alleles map uses the shared gene-to-major convention. These are preliminary evidence, not assignments or recombination boundaries. Returns flat Polars queries and alignments frames plus limits. Calls are bounded to batch_size inputs; iterate batches for large cohorts. Coordinates are 1-based closed in the input query and segment germline. V/J uses one batched forward-strand MMseqs segment search; callers supply coding-strand sequences. max_seqs caps prefilter targets, not alleles. MMseqs returns one local alignment per target subject to evalue. D retains each allele’s best positive (+1 match/-1 mismatch) ungapped span per allele, without an E-value gate or greedy D/D2 selection. Local endpoints are evidence only; inference must reconsider trimming and N/P. threads budgets the sequential MMseqs and native D phases. D uses contiguous query spans; batches below 32 queries run serially. ARDA_MMSEQS_SEARCH_OPTS must be unset; search settings are explicit here. No claim of exhaustive alignments or recombination scenarios is made.

Major allele naming

arda.select_major_alleles(alleles)

Return gene -> major allele, preferring *01 then the lowest numeric allele.

Names without a star are valid (notably imported D gene references). This function imposes no functionality/completeness policy on the supplied names.

arda.major_allele_map(organism='human', locus=None)

Return usable runtime gene -> major allele for an organism or one locus.

Prefer ok F/ORF anchors among available V/J targets. Genes with only incomplete targets retain one preliminary alignment template, with its original anchor status; this never establishes complete emission eligibility. D names come from exported germlines. No raw IMGT installation is needed.

arda.d_germline_metadata(organism='human')

Return exact D-template metadata keyed by (allele, nucleotide sequence).

Preliminary alignment targets may be partial; complete emission consumers must reject status=partial. Functionality policy belongs to the consumer.