v4.8.0 — the junction pipeline

The junction pipeline

model.annotate_junctions(cdr3_aas, v, j, species=) takes a bare amino-acid junction with its V and J calls — all a curated database has — and returns a repaired junction, confirmed or corrected gene calls, V/J boundaries in residues and nucleotides, the most plausible nucleotide junction, and a D with coordinates in both alphabets. One row out per row in, in input order.

Four stages, each an existing entry point: arda repairs the junction and places the boundaries, infer_nt_batch infers the nucleotides (the templated flanks are written from germline, so only the N region is inferred), the model names the D gene, and the aligner places it on those nucleotides — whose coordinates then fold onto the residues their codons touch. That last step is the point: a D contributing one amino acid is invisible in a translated junction but not in nucleotides, because the residues at each end of it are part germline and part N region.

Naming the D and placing it are separate questions

The model names the gene (d_call, with d_posterior) from its own scenario weights — no second model, no fitted prior table. The aligner then places that gene greedily and ungated (d_start_nt / d_end_nt, with np1 / np2 either side), because the gene is already chosen and refusing to say where it sits only leaves a row with nothing to draw.

On 4,000 real human TRB rearrangements whose D comes from the nucleotides: the gene is right on 74.35 % of all rows and 99.80 % have coordinates. Letting the aligner choose the gene as well reaches 47.93 % gated and 67.10 % ungated, and it is no better even where it is confident — 85.96 % against the model's 86.32 % on those rows. So there is no d_best, no d_posterior_call, no d_entropy and no arda.dpost port: a second D estimator was measured and is not needed.

The nucleotides are the authority

v_end_nt, j_start_nt, d_start_nt and d_end_nt are all read off the inferred nucleotide junction, and every amino-acid bound is recomputed from them, so a view showing both alphabets cannot draw them disagreeing.

model_source="auto" is a measured chain, not a preference

OLGA's bundled fit first, then arda's on what it left empty. Human TRB nucleotide-exact 14.40 / 17.32 / 17.32 % and D gene 72.58 / 74.08 / 74.35 % for arda / OLGA / the chain, and human TRA keeps 4,000 of 4,000 nucleotide junctions where OLGA alone declines 143.

Every species answers, not just the fitted ones

A fitted model exists for human (seven loci) and mouse (TRA/TRB) and for nothing else, so every other organism came back with no nucleotides, no D and no bounds — 1,457 rhesus keys in VDJdb answered zero times, invisible inside a single corpus-wide total. The chain now ends in a germline scaffold (from_arda(locus, organism=)): the templated flanks come from that organism's own germline and only the N region is a default.

species	locus	keys	nucleotide junction	D gene	D coordinates
human	TRB	115,829	115,819	115,819	113,232
human	TRA	58,281	58,203	—	—
mouse	TRB	9,015	8,805	8,805	8,694
mouse	TRA	8,140	8,140	—	—
rhesus	TRB	1,383	1,379	1,379	787
rhesus	TRA	74	73	—	—
human	TRD	4	4	4	4

192,423 of 192,726 keys (99.84 %) with a nucleotide junction, against 124,022 before, and no
species or locus at zero. TRA and TRD are VJ loci here — there is no D to find.

Speed

122 µs per key on the real corpus — VDJdb's 192,726 distinct (species, cdr3, V, J) keys in 23.6 s in one process. Every stage is batched (one markup_batch, then one native threaded infer_nt_batch and one best_aa_scenarios_batch per organism, locus and model set, on a model loaded once), so do not wrap the call in a pool of your own.

Blank calls

A blank V or J is handled, not refused: arda proposes that side from the junction and proposed says which side was never curated. 3,130 of the 192,726 keys leave a side out — 2,669 name one side, 461 name neither, and for those the locus is proposed too (agreeing with the column the record was filed under on 457 of 461, 99.13 %).

Fixed

• infer_nt_batch raised on a batch whose first row had no D. The native call returns its string columns as numpy object arrays and polars types one from its first element alone, so [None, "TRBD1*01"] came back Object and the declared Utf8 cast raised cannot cast 'Object' type — while ["TRBD1*01", None] was fine. The failure depended on row order.
• infer_nt_batch raised on a per-row list of alleles, one of the three call forms it documents, when the boundary columns were computed.
• A conditioned V could return nothing. A bundled model is an EM fit and an allele its training cohort never showed comes out at probability zero — 36 of 66 human TRB V alleles — so conditioning on one produced no scenario at all. The pipeline floors those before inferring.

⚠ Two things that must not be quoted wrongly

OLGA is not human-only; this bundle is. OLGA ships mouse_T_alpha, mouse_T_beta, mouse_B_heavy, mouse_B_kappa and mouse_B_lambda, and all five import cleanly through from_olga — they are simply not vendored here. On mouse TRB the vendored arda fit is the better one anyway (D gene 53.97 % against 52.33 %, 4,000 rearrangements), so what vendoring them would buy is mouse IGH/IGK/IGL, which no bundled set covers.

Every accuracy figure above is human TRB. The truth set (isalgo/airr_control) carries TRA and TRB and no immunoglobulin, so none of them may be quoted for IGH. The D call does reach IGH, but this pipeline has no somatic-hypermutation term, so a hypermutated V tail is priced as insertion, which moves v_end_nt and therefore where the D can sit. For B cells with real nucleotides that is arda.hmm(..., shm=ShmModel) — a complement to this, not an alternative.

Requires arda-mapper>=2.36.0.

Full changelog: https://github.com/antigenomics/vdjtools/blob/master/CHANGELOG.md