arda

name	arda
description	Fast TCR/BCR FR/CDR region annotation (Antigen Receptor Domain Annotation). Use whenever the user wants to annotate immune-receptor sequences with framework/CDR regions, V/D/J gene calls, constant-region isotype (c_call/c_class), junction/CDR3 boundaries, or AIRR output — for TCR (TRA/TRB/TRG/TRD) or BCR (IGH/IGK/IGL), nucleotide or amino acid, single sequences or large FASTA/FASTQ. Also use for: marking up, validating or REPAIRING a bare (CDR3 amino acid, V, J) record that has no read behind it — a VDJdb-style row — including locating where it disagrees with germline and restoring a missing Cys104 / Phe118 anchor; calling the D gene (and tandem D-D) on a bare junction, or inferring it from junction length when the protein shows none of it; extracting the receptor repertoire (clonotypes, isotype usage) from bulk RNA-seq; getting germline FR1-FR3/CDR1-CDR2 (V) or FR4 (J) subsequences for individual alleles; building or rebuilding the reference database from IMGT germlines via IgBLAST; or diagnosing mmseqs2 setup. arda runs MMseqs2 + a C++ coordinate projection (IgBLAST is offline/build-time only), so it is much faster than IgBLAST at annotation time. Load references/ files for the detailed API, region/junction semantics, reference-build pipeline, or mmseqs install/troubleshooting.
license	GPL-3.0
compatibility	Python 3.10+; `pip install arda-mapper` (>=2.5.0 for junction markup / repair, the aa D posterior, and D on protein input) — no source checkout and **no `ARDA_HOME`**: the curated `vdj/` reference auto-fetches into `~/.cache/arda` on first use (set `ARDA_NO_AUTO_FETCH` for air-gapped runs), and the `mmseqs` binary auto-fetches a static build into the cache if missing — so a bare `pip install` annotates out of the box. Bulk RNA-seq needs nothing extra: `seqtree` is a core dependency since 2.5.5. (Before that it lived in an optional `[rnaseq]` extra, and a plain install would map and assemble a whole sample and only then die, before writing any clonotype table.) A source checkout / `$ARDA_HOME` still uses the committed `database/`. Shell is fish — use fish syntax in terminal commands.
metadata	{"repo":"https://github.com/antigenomics/arda"}

arda Skills Guide

arda annotates the framework (FR1–FR4) and complementarity-determining (CDR1–CDR3) regions of TCR/BCR sequences. The expensive IgBLAST markup is done once, offline, when the reference database is built; at annotation time arda only runs an MMseqs2 search + a C++ routine that projects the reference region coordinates onto each query. That makes it embeddable and ~4–8× faster than IgBLAST, with 98–99.7% region concordance on real GenBank mRNA across all five organisms.

It also handles records with no read behind them — a CDR3 amino acid plus a V and J call, as in VDJdb — marking up which residues each germline templates, repairing the junction, and inferring the D gene from the junction's length.

Core API

import arda

records = arda.annotate_sequences(
    ["GACGTGCAG...", ("clone7", "CAGGTG...")],  # raw strings or (id, seq) pairs
    seqtype="nt",          # "nt" or "aa"
    organism="human",      # human | mouse | rat | rabbit | rhesus_monkey
    map_d=True,            # map D segments for VDJ loci — works on aa input too
)
# -> list of AIRR record dicts (one per query)

For explicit control of strand / sensitivity / in-memory vs file streaming, use the mapper directly:

from arda.annotate.mapper import annotate_records, annotate_file

recs = annotate_records(queries, organism="human", seqtype="nt",
                        strand="forward", map_d=False, sensitivity=7.0)
annotate_file("reads.fastq.gz", "out.airr.tsv", organism="human")  # streamed, memory-flat

Each record dict carries (1-based closed coords, query space): locus, v_call/d_call/d2_call/j_call, the constant-region c_call/c_class (isotype), productive/stop_codon/vj_in_frame, rev_comp, v_identity, sequence_alignment/germline_alignment, {v,j,c,d}_cigar, *_germline_start/end, v_sequence_end, j_sequence_start, np1/np2/np3, d_support/d2_support, junction(_aa), and per region in (fwr1, cdr1, fwr2, cdr2, fwr3, cdr3, fwr4): {r}_start, {r}_end, {r}, {r}_aa. Ambiguous D and C calls are comma-joined allele lists (as V/J are). The TSV is a spec-valid AIRR Rearrangement file (passes airr.schema validation).

d_support is the Karlin–Altschul E-value the D call was gated on (accepted at <= 0.2 for nt, <= 0.05 for aa). It ships so a consumer can re-threshold: keeping rows with d_support <= x for x < 0.2 reproduces exactly what a stricter arda would have called. A missing d_call on a VDJ locus usually means the best hit did not clear the gate, not that mapping was skipped.

Read references/annotation.md for the full field list, parameter semantics (strand/sensitivity/threads/chunking), AIRR column order, the D E-value gate + genomic-order constraint, and performance notes.

Batch annotation — never loop (use mmseqs2's own parallelism)

Always gather every sequence first, make ONE annotate_sequences call, then do downstream analysis on the batch output. Each annotate_* call pays a fixed ~825ms mmseqs2 process+index-load cost; a batch of 300 sequences costs the same ~930ms total because mmseqs2 parallelises internally across threads. So:

# RIGHT — one batched call, mmseqs2 threads internally
recs = arda.annotate_sequences([(cid, seq) for cid, seq in all_chains], organism="human")
by_id = {r["sequence_id"]: r for r in recs}     # then map back per-item, downstream

Do not wrap per-item annotate_* in a Python ProcessPoolExecutor/ThreadPoolExecutor or a loop: a process pool that forks after mmseqs2/BLAS have spawned threads deadlocks, a thread pool just serialises on the same overhead, and either way you pay the fixed cost N times instead of once. mmseqs2 is the parallel layer — Python orchestration is single-call.

Region & junction semantics

• Region coordinates are projected through the MMseqs2 alignment, so they are correct even for truncated, mutated, or reverse-strand queries.
• There is no coverage filter: a partial read (or a bare germline V or J) maps to its scaffold and returns only the regions inside its coverage. A bare V → fwr1..fwr3; a bare J → fwr4. This is how callers get per-allele germline FR/CDR subsequences without synthesising a rearrangement.
• junction spans Cys104 through the [FW]118 that opens FR4; cdr3 is J-anchored. Out-of-frame junctions are reported with an N-bridge (_).

Read references/region-segments.md for the bare-germline recipe, junction/CDR3 details, and coordinate round-trip rules.

Bare records — a CDR3 amino acid, a V call, a J call, no read

VDJdb-style rows have nothing to align. The V and J germlines still template a known run of residues into each end of the junction, and arda ships those per allele.

from arda.cdr3fix import markup_cdr3, markup_records   # markup_records: a whole polars frame
from arda.annotate.dmap import map_d_junction          # D (+ tandem D-D) on a bare nt junction
from arda.dpost import posterior_d                     # D gene + position from junction LENGTH

mk = markup_cdr3("CAIRDDKII", "TRAV12-3*01", "TRAJ30*01", "HomoSapiens")
mk.cdr3_repaired             # 'CAIRDDKIIF'  -- the Phe118 anchor restored
mk.v_end, mk.j_start         # residues templated by V / index of the first J residue
[str(e) for e in mk.errors]  # ["J del@8 missing 'F' d=0"]
mk.good                      # both sides repaired AND both anchors present
mk.to_cdr3fix()              # VDJdb's `cdr3fix` JSON object, key for key

CLI: arda markup -i vdjdb.txt -o marked.tsv --vdjdb --report - [--d-posterior].

The single biggest correctness trap. These coordinates are junction space: Cys104 through Phe/Trp118, both anchors included. That is what VDJdb's cdr3 column holds. It is not arda's cdr3 field, which excludes both — junction_aa is two residues longer than cdr3_aa. Conflating them silently corrupts every coordinate, and downstream corrupts Pgen, clustering and matching.

Repair is deliberately conservative and its two decisions are separate:

• Every germline disagreement is reported (side, kind, position, extent, distance from the anchor). Only edits adjacent to a conserved anchor are applied; deeper ones are left alone, because there a mismatch is as likely to be the real V/N boundary as a typo. Cdr3Error.applied is true only when the edit reached cdr3_repaired.
• A repair always lands on a canonical junction. If the result would not open with Cys104 and close with Phe/Trp118, it is refused and the submission returned untouched. So good implies canonical. An allele with no derivable anchor gives FailedBadSegment — flagged, never guessed.

posterior_d infers the D gene and where it sits from the junction's nucleotide length, which pins insVD + |D surviving| + insDJ. Shipped for human IGH/TRB/TRD and mouse TRB only (the pairs with a published generative model); every other pair returns None rather than guessing — do not substitute a human proxy.

Organisms & loci

Organism	Loci with full markup
human, mouse	TRA, TRB, TRG, TRD, IGH, IGK, IGL
rat, rabbit, rhesus_monkey	IGH, IGK, IGL (IG only)

VDJ loci (D segments mapped): IGH, TRB, TRD. D-D fusions sought in all three. D mapping runs on protein input too, against each D germline's three translated frames — useful for IGH (a call on ~36% of real records, agreeing with the nucleotide call on 98% of them), mostly silent for the TR loci, whose D is too short to survive trimming into protein. On aa input d_germline_* and d_cigar stay empty: those offsets index a reading frame, not the germline.

Genomic order constrains the call. TRBD2 lies 3′ of the entire TRBJ1 cluster, and V(D)J joining deletes the intervening DNA, so a TRBJ1 rearrangement is never assigned TRBD2 — in any species with that architecture. IGH and TRD place every D 5′ of every J, so nothing is excluded there.

Constant-region J + C scaffolds (isotype c_call/c_class) are built for every locus with a CH1 exon in the bundle.

CLI

arda info                                   # versions + available references
arda annotate -i reads.fastq.gz -o out.airr.tsv --organism human --seqtype nt
arda annotate -i prot.fasta -o out.tsv --seqtype aa --no-map-d
arda markup -i vdjdb.txt -o marked.tsv --vdjdb --report -   # bare (CDR3aa, V, J) records
arda rnaseq map --r1 R1.fq.gz --r2 R2.fq.gz -o mapped.airr.tsv   # bulk RNA-seq: filter receptor reads
arda rnaseq assemble -i mapped.airr.tsv -o assembled.airr.tsv    # rescue CDR3s no read spans
arda rnaseq correct -i mapped.airr.tsv -o clones.tsv            # collapse CDR3 errors into clonotypes
arda rnaseq run --r1 R1.fq.gz --r2 R2.fq.gz -p SAMPLE -d out/   # one-shot map+assemble+correct
arda igblast -i reads.fastq -o truth.airr.tsv                   # gold-standard IgBLAST (all loci)
arda build-db --organism all                # offline reference build (needs IgBLAST)
arda build-index --organism all             # rebuild mmseqs indexes for local mmseqs version
arda slurm -i big.fastq -o big.airr.tsv --shards 50   # multi-node: split → array → merge

Bulk RNA-seq mode (arda rnaseq)

For libraries where only 1–5% of reads are receptor-derived. Three stages, run separately or in one shot with rnaseq run (which does all three by default). Needs the rnaseq extra: pip install arda-mapper.

map — streams paired FASTQ (--r1/--r2), keeps only reads mapping to a receptor scaffold, writes them as AIRR. Recall-first, with --min-score/--kmer/--max-seqs around one default preset.

• The reference includes J + C constant-region scaffolds, so a read spanning the J→C splice (no V, hence no junction) still maps and carries c_call/c_class. In paired mode a CDR3-bearing read gets its isotype from its constant-region mate.
• --reconstruct merges overlapping mates into one fragment, giving a short read the mate's V/J context; overlap mismatches resolve to the higher-Phred base. FASTQ quality is read only on this path, so the default stays fast.

assemble (Stage 3) — recovers clonotypes whose CDR3 no single 100–150 bp read spans (V(DD)J ultralong, ~20–40 aa), by anchored greedy overlap-extension over Stage-1's per-read cdr3_start. It carries the contig's D call onto every member read: an ultralong CDR3 is where a tandem D-D is both most likely and least visible to one read.

annotate.contig gives an assembled contig its AIRR cigars two ways, producing the same record: reannotate_contigs (re-align it — what assemble uses) and merge_contig (stitch the reads' existing alignments via C++ _markup.merge_alignment). Merge is ~9× faster at ~10⁵ contigs/sample (scRNA-seq) and is the intended default once the assembler emits read layouts.

correct — collapses sequencing-error CDR3 variants into clonotypes keyed by (locus, v_call, j_call, junction).

• Abundance is the AIRR duplicate_count (every read encompassing the junction), with consensus_count for distinct fragment consensuses. There is no count column.
• A neighbour is an error child when count[parent] * p_sub**n_subs * p_ind**n_indel >= count[child]. Knobs: --max-subs, --max-indel, --error-rate, --indel-rate (per-BASE, length-scaled), --require-vj, --error-method (simple|binom|betabinom), --complete-only (on by default).
• Row order is deterministic — abundance ties break on (junction, v_call, j_call).
• Each clonotype's D is mapped once into its corrected junction (d_call/d2_call/ d_support), not voted over reads: D is a function of the junction, and a read's copy of it carries sequencing error.

mmseqs2 (auto-installed)

Annotation needs the mmseqs binary. Resolution order: $ARDA_MMSEQS → <project>/bin/mmseqs → mmseqs on PATH → auto-fetch a static binary into bin/mmseqs. So neither conda nor pip users must install it manually. The conda env (environment.yml) also ships mmseqs2 from bioconda.

Read references/install-mmseqs.md for env vars (ARDA_MMSEQS, ARDA_MMSEQS_ASSET, ARDA_NO_AUTO_FETCH), the shipped/precompiled indexes, and version-mismatch handling.

Rebuilding the reference

Most users never build anything — database/vdj/<organism>/ ships with precompiled markup and MMseqs2 indexes. Rebuild only when adding/refreshing an organism (needs IgBLAST, fetched by setup.sh into bin/).

Every build writes loci_manifest.tsv — one row per defined locus (V-J / J+C scaffold counts, D germlines, unreachable-D count, ok/EMPTY status) — and warns at build end on any EMPTY locus or unreachable D germlines. This is what makes an absent reference visible rather than silent: rat/rabbit/rhesus have no TR loci in IMGT, so their TCR loci build EMPTY (the IG-only limitation in the table above).

Read references/reference-build.md for the arda.refbuild pipeline (IMGT germlines → V×J scaffolds → IgBLAST → markup TSVs) and build-db / build-index.

Sequence primitives

arda.refbuild.translate exposes fast C++-backed helpers, mirpy-API-compatible: translate(nt, frame=0), detect_coding_frame(nt), reverse_complement(nt), back_translate(aa), aa_coords_from_nt(nt_start, nt_end, coding_start).

Gotchas

• junction is not cdr3. junction/junction_aa include both conserved anchors; cdr3/cdr3_aa exclude both, so cdr3_aa == junction_aa[1:-1] always. Everything in arda.cdr3fix / dmap / dpost works in junction space, matching VDJdb's cdr3 column. Mixing the two conventions is the most expensive mistake available here.
• An empty d_call is a decision, not a gap. The call is gated on d_support (E-value ≤ 0.2 nt, ≤ 0.05 aa). Human TRB gets a D on only ~47% of junctions because an ordinary TRB interior is 11–21 nt and heavily-trimmed TRBD1 scores below the gate. VJ loci (TRA, TRG, IGK, IGL) have no D gene at all.
• aa input returns region *_aa directly with no frame bridging, so stop_codon and vj_in_frame stay empty — but productive and the D columns are populated.
• posterior_d returns None for organisms with no shipped generative model (rat, rabbit, rhesus) and for VJ loci. That is deliberate. Do not fall back to a human model.
• The shipped MMseqs2 indexes are used only when the local mmseqs version matches; otherwise arda rebuilds a private cache in data/ on first run. arda build-index (re)builds the shipped indexes for your version.
• map_d=True on synthetic/partial input with no real junction simply finds no D — harmless; pass map_d=False to skip the search.
• IgBLAST is needed only to build references, never at annotation time.
• arda rnaseq correct uses seqtree (a core dependency since 2.5.5); without it every correct test skips silently rather than failing.