2.33.0 — the repair policy, and the D posterior moves to vdjtools - #145
Merged
Merged
Conversation
`arda.cdr3fix` replaced VDJdb's `Cdr3Fixer`. On one class of input it detects the defect, names it exactly, and does not apply the repair, so a malformed junction reaches the consumer where the retired fixer corrected it. Measured over VDJdb's corpus, 190,902 distinct (species, cdr3, v, j) keys: - 382 keys end up with a junction that does not run Cys104 to the anchor its own segment encodes (267 J, 115 V); - 2,165 keys are returned `NoFixNeeded` with a single residue of germline agreement, where legacy required a k-mer hit of `min_hit_size = 2` at offset zero in both sequences. The asymmetry is the clearest symptom: `YFCASSQSPGGVAFFGQG` on TRBV14/TRBJ1-1 gives `v_fix=FixTrim` and `j_fix=FailedNoAlignment`, while `YFCAVVGTGLGYTFGSG` on TRBV9/TRBJ1-2 gives exactly the opposite. Same defect at both ends, decision flips per end. The cause is the reference, not the algorithm. `Anchor.templated_aa` runs Cys104 through [FW]118 inclusive and stops, so framework past an anchor has nothing to align to and the boundary can only be inferred from the junction side. The legacy fixer sliced `sequence[reference_point - 3:]` for V and `sequence[:reference_point + 4]` for J, putting one FR3 codon before Cys104 and the [FW]GXG motif after 118 - with the flanks present the repair is positional and needs no alignment scoring. arda already ships what that needs: `alleles.fasta` holds full allele sequences and `cdr3_anchors.tsv` carries `anchor_nt` as the reference point. Files: - `ISSUE-cdr3fix-kmer-regression.md` - the writeup, per class, with V and J calls; - `cdr3fix-noncanonical-junctions.tsv` - the 382, classified; - `cdr3fix-nofixneeded-without-alignment.tsv` - the 2,165; - `reference-legacy-cdr3fixer/` - the original `Cdr3Fixer`, `KmerScanner`, `Utils` and `FixerDataModels`, vendored from vdjdb-db before #658 retired them; - `src/arda/cdr3kmer.py` - scanner and translation ported verbatim, four-outcome table documented. **Incomplete**: no flank-extended segment construction from `alleles.fasta` + `anchor_nt`, no `fix()`, no batch API, no tests. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The 2.16.0 semi-global Needleman-Wunsch engine is gone. One gapless local alignment
per side (`_markup.d_local_align`, the same C++ the D caller uses) places the germline's
templated run anywhere in the junction, and legacy's positional table decides the
outcome from the two offsets.
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 an engine anchored at
offset 0 scored the first and called the record clean. The anchored placement still
wins a tie, or `CASSQQQQQQQQQF` gains an F it already had.
Measured against the authoritative VDJdb 2026-06-03 release, whose `cdr3fix` column
carries the retired fixer's own cdr3_old -> cdr3. Coverage first: 184,765 of 189,596
distinct (species, junction, V, J) curation keys join (97.45 %).
2.31.0 2.32.0
repaired junction agrees 180,374 (97.6235 %) 181,892 (98.4451 %)
vEnd agrees 98.2013 % 98.3571 %
jStart agrees 97.8316 % 99.7242 %
release's own repairs 584 of 779 623 of 779
repairs the release did not 4,196 2,717
good beside an unrepaired 19,436 0
throughput, one process 12,900 keys/s 20,400 keys/s
The verdict is now a SET of flags per side (`v_flags`, `j_flags`): ok, allele, sub,
add, trim, shallow, impossible. One worst-wins label could not say both "I trimmed a
flank" and "I found a substitution I will not touch", and collapsing them is how a
declined repair came to read as a clean record — `CAISGEFGSGA` reported `V sub@2 I>S
d=2` and still returned NoFixNeeded and good, as did 19,436 keys. `good` is now
exactly "neither side impossible"; `v_fix`/`j_fix` stay VDJdb's names for the JSON.
A contradicted call changes the ALLELE, not the sequence (`guess_allele`, legacy's
`guess_id` rebuilt on the same scan). Of the anchor-adjacent substitutions 2.31.0
wrote into an already-canonical junction, 74.7 % were records whose own 3' end matched
a different J allele better than the called one, against 1.1 % of untouched records —
a 68x enrichment. `CASSLRGAATDTQYF` is a clean TRBJ2-3 junction called TRBJ2-1 that
2.31.0 rewrote to `CASSLRGAATDTQFF`, a string no germline supports.
arda does not force a junction canonical on a thin match: writing the conserved anchor
takes a hit clearing `_MIN_HIT`, and `_canonicalise` scores against the called allele's
OWN anchor residue, not the [FW] motif — TRBJ2-7*02 templates `SYEQYV`.
A call is also checked against every gene name in the vocabulary, which was the largest
single failure class: 3,067 V and 152 J keys that came back FailedBadSegment. TRAV14
(999) is filed as TRAV14/DV4, TRBV19-1*01 (246) names a suffix IMGT does not use,
TRBJ2.1 (5) is dot nomenclature. A spelling must land on exactly one gene or it is
refused.
Removed `_MAX_TRIM` and `_MAX_FIX`, two bounds layered on top of `max_replace`: the
release trims framework `_MAX_TRIM` refused, and two knobs on one decision meant the
outcome depended on which bit first.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The repair policy, stated by the author on 2026-09-30 and now the module's rule. Four
edits are admissible, in this order, and nothing else is ever rewritten:
1. re-call the V or J when another allele explains 2 more residues CONTIGUOUSLY from
its own anchor (`anchor_depth`, `_RECALL_GAIN`), and does not break an anchor the
submission already holds;
2. trim framework past that allele's anchor;
3. add back germline the submission was cut inside of;
4. substitute the FIRST OR LAST residue only, to the anchor, on a solid match.
A residue that disagrees with germline INSIDE the templated run is reported (new
`mismatch` flag) and left exactly as submitted, because a curation error and an allele
IMGT does not record are indistinguishable from one junction. 2.32.0 rewrote them:
`CAAAETSYDKV{M,R,T,V}F` on TRAJ50 are four junctions that all became `CAAAETSYDKVIF`,
which is four residues observed at one templated position — a suspect J call or an
unrecorded allele, not four independent typos.
Re-calling is now step ONE for every record rather than a rescan triggered by a
substitution, which it can be because `_candidates` indexes the anchor table by
(segment, locus) and dedupes by templated run: the whole-table Python filter cost 19x
the markup itself, 110 s against 18 s over 187,488 keys.
`_RECALL_GAIN = 2` is measured against `isalgo/airr_control`'s nucleotide-called V/J
(junctions seen in >= 3 donors with one V and one J throughout), on the 3,059 human TRA
and 426 human TRB curation keys that appear there. It is the smallest margin at which no
re-call moves a call AWAY from the nucleotide answer, and every larger margin only
leaves errors standing:
TRA V TRA J TRB V TRB J
no re-call 88.17 % 99.18 % 94.13 % 98.12 %
gain 1 90.68 % (14 wrong) 99.87 % (0 wrong) 96.71 % (4 wrong) 99.77 % (0 wrong)
gain 2 89.47 % (0) 99.80 % (0) 96.71 % (0) 99.77 % (0)
gain 3 88.92 % 99.80 % 94.60 % 99.77 %
Three guards, each because the alternative was measured to re-call a record wrongly.
The hit must start at the germline's OWN anchor (`CSAR` is TRBV20-1's run and it does sit
inside `CGGSARSGELFF`, at offset 3, as N region). It may not place that anchor further
inside the junction than the call does (TRBJ2-2P's `LRGAAG` matches seven residues in and
trimmed `CASRPGAAGGRPELYF` to `CASRPGAAG`). And it may not break an anchor the submission
already explains: `CAARLGNNYKLIW` is a clean TRAJ33 junction whose `WIL` is contiguous
from its own W, which TRAJ12's five residues outscored only by skipping the anchor —
whereupon step 4 turned the W into an F. `_canonicalise` cannot catch that, because BOTH
anchors are canonical. The converse stays allowed: `CASSKRGGYEQYV` on TRBJ2-7*01 does not
hold its anchor, so *02 takes the call at equal depth.
A junction that lost only its conserved anchor is now repaired. Neither legacy nor 2.32.0
did: the survivors agree from position 1 on, so the best placement is a one-residue
coincidence below `_MIN_HIT`. Prepending the anchor is admissible exactly when doing so
makes the germline agree from the anchor outward — `AQGLLTGGGNKLTF` on TRAV29/DV5 (`CAAS`)
becomes `CAQGLLTGGGNKLTF`, where the 2026-06-03 release ships it unrepaired and not good,
while `CALRPA` on TRAJ17 stays refused because prepending its Phe still leaves one residue.
`good` is "neither side impossible", and `mismatch` keeps it: the junction is well formed
and every residue in it is the curator's. `impossible` is now exactly an ANCHOR the rules
cannot restore. `max_replace` is legacy's `max_replace_size` again — junction residues
standing where the anchor should be — not 2.32.0's distance from the anchor.
Measured against the authoritative 2026-06-03 release, all 187,488 distinct
(species, cdr3_old, vId, jId) keys, both engines given the same four inputs
(`scripts/compare_vdjdb_release.py`, new):
2.32.0 this
repaired junction agrees 98.3604 % 99.7840 %
the release's own repairs made 4,267/4,499 4,389/4,499
rewrites the release never ships 2,842 295
of those, already canonical 2,744 186
of those, made canonical 69 98
non-canonical junctions shipped 709 708 (release: 715)
good 170,218 186,139 (release: 174,407)
throughput, one process 32,400 keys/s 10,100 keys/s
vEnd concordance moves 99.86 % -> 98.19 % and jStart 99.65 % -> 98.88 %, entirely on
records whose call arda re-assigns. On those the release's own boundary is wrong 100 % of
the time against nucleotide truth and arda's is right 25-100 % (TRA vEnd 0 % -> 25.42 %,
TRA jStart 0 % -> 78.26 %, TRB jStart 0 % -> 100 %); on confirmed calls arda matches or
beats it everywhere.
Also restores CLAUDE.md, which ae112a8 deleted.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… 3.1x, and better
Where the germline run begins is a property of the JUNCTION: a framework flank in front of
Cys104 is there whatever allele is named. Scoring each candidate at a placement of its own
let one win by putting its anchor DEEPER inside the junction, where the trim that follows
deletes everything before it — `TRBJ2-2P`'s `LRGAAG` matches `GAAG` seven residues into
`CASRPGAAGGRPELYF` and trimmed the record to `CASRPGAAG`. Every candidate is now scored at
the called allele's own offset, which is one `anchor_depth` walk instead of a `scan` (two
`_extend` sweeps plus a C++ alignment).
`_extend` was 1.4 M calls over 20,000 records, 36 % of the markup, and a candidate index on
the offset-1 residue — the one position every route to `anchor_depth >= 2` must match — cuts
the ~60 alleles of a locus to three to five. Together: 9,200 -> 28,900 keys/s over 187,488
keys, within 15 % of doing no re-call at all.
Measured against the authoritative 2026-06-03 release (187,488 keys), and every number
moves the right way:
2.32.0 prev commit this
repaired junction agrees 98.3604 % 99.7840 % 99.8352 %
rewrites the release never ships 2,842 295 141
of those, already canonical 2,744 186 35
of those, made canonical 69 98 95
non-canonical junctions shipped 709 708 677 (release: 715)
arda now ships FEWER malformed junctions than the release it is measured against. The
re-call's accuracy against `isalgo/airr_control`'s nucleotide-called V/J is unchanged to
four places on all four columns (TRA V 89.47 %, TRA J 99.80 %, TRB V 96.71 %, TRB J
99.77 %, zero calls moved the wrong way).
The release's own repairs reproduced goes 4,389 -> 4,331 of 4,499: 58 records where its
repair followed from a placement this no longer considers. Deliberate, and the same trade
that removed 154 rewrites it does not ship.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ion cannot resolve Closes #144, #142 and antigenomics/vdjtools#183. `arda.dpost` is gone: a posterior over the D gene marginalises the generative model's insertion-length and D-trimming distributions and multiplies in P(D | J), which are model marginals, not germline facts. It is now `vdjtools.model.posterior_d` plus `posterior_d_batch`, the batch entry point #142 asked for. Ported, not rewritten — over 3,000 real human TRB junctions the new module agrees with this one on every field of every row, so the measured accuracy travels intact. `arda markup --d-posterior` / `--d-prior` are gone with it, and a removed flag now FAILS rather than being accepted and ignored. What stays here, because arda is the germline reference and the markup against it: the prior table `database/vdj/<org>/d_prior.tsv` (its other two consumers are `arda.hmm` and `arda.scenarios`, so a copy in vdjtools would be a second copy of a fitted artifact), the fitter, and now its READER — `arda.scenarios.load_prior_table`, beside the writer, so the format has exactly one parser. vdjtools reads the table through it. Two additions the junction pipeline needs, both of them arda's own job: `v_alts` / `j_alts` — every allele that explains an end exactly as well as the chosen one. `CAISE` is the run of TRBV10-3*01, *02 AND *03, so an amino-acid junction cannot separate them, and resolving that by functionality then by name binds a choice with no evidence behind it. The set travels instead, chosen allele first, and the stage that CAN separate them is `vdjtools.model.infer_nt_batch`, which scores a list per row. `map_d_junction(v_end=, j_start=)` — the interior to search, in nucleotides, instead of re-deriving it by exact germline prefix match. That derivation is right for a read and wrong for a junction whose nucleotides were INFERRED under a model keyed on a different allele of the same gene: the prefix 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. Deprecated: `arda.hmm`. Nothing consumes it — not the annotation path (two measured negatives in ROADMAP.md), not another repository, not the new pipeline — and the question it answers, the most likely nucleotide reading of a junction, is `infer_nt_batch`: batched, threaded and C++. `arda.scenarios`, which fits the model, is NOT deprecated: it is measured to still be needed, since the length-and-prior posterior it parameterises answers 1,753 of 4,000 real human TRB rearrangements that the nucleotide alignment declines, and is right on 53 % of them where the alignment reaches 86 % on the rest. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`arda.hmm`'s removal and the `d_prior.tsv` location travel together (loop 1 gates loop 2); the 141 keys where arda edits a junction the authoritative release does not have never been read record by record; and the consumer migration is tracked in the consumer (antigenomics/vdjdb-db#713). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The asset still ships `d_prior.tsv` and still must — `arda.hmm` and `arda.scenarios` read it, and vdjtools reads this copy rather than vendoring a second copy of a fitted artifact. Only the POSTERIOR moved (#144), so the check is unchanged and the comment beside it was stale. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #141, #142, #143, #144.
The repair policy: when it is not certain, flag it — do not change it
cdr3fixadmits exactly four edits, in this order, and rewrites nothing else:A residue that disagrees with germline inside the templated run is reported (
mismatch) and left as submitted — a curation error and an allele IMGT does not record are indistinguishable from one junction. 2.32.0 rewrote them, which turnedCAAAETSYDKV{M,R,T,V}Fon TRAJ50 into four copies ofCAAAETSYDKVIF._RECALL_GAIN = 2is measured againstisalgo/airr_control's nucleotide-called V/J on 3,059 human TRA and 426 human TRB curation keys: the smallest margin at which no re-call moves a call away from the nucleotide answer (at 1 it breaks 20).Against the authoritative 2026-06-03 release, all 187,488 keys
arda ships fewer malformed junctions than the release it is measured against.
vEnd/jStartconcordance falls only on records whose call arda re-assigns, where the release's own boundary is wrong 100 % of the time against nucleotide truth. Reproduce withscripts/compare_vdjdb_release.py.arda.dpost→vdjtools.model.posterior_d(#144, #142)A posterior over the D gene marginalises model distributions, so it belongs with the recombination model. Ported, not rewritten: it agrees field-for-field on 3,000 real human TRB junctions, and it gained the batch entry point #142 asked for.
arda markup --d-posterior/--d-priorare gone. The prior table, its fitter and now its reader (arda.scenarios.load_prior_table) stay here.#143 —
_align, the pure-Python DPGone with the 2.32.0 port; there is no Python alignment DP in
cdr3fixany more. Throughput over the same 187,488 keys is 28,900 keys/s including the new per-record re-call.Also
v_alts/j_alts: the alleles a junction cannot separate travel with the answer, for vdjtools' nucleotide stage to settle.map_d_junction(v_end=, j_start=): take the interior instead of re-deriving it from an inferred sequence.arda.hmmdeprecated — nothing consumes it, andvdjtools.model.infer_nt_batchanswers the same question batched and in C++.1,312 tests pass.
🤖 Generated with Claude Code