arda 2.29.0 — a model of where hypermutation is expected
Added: arda shm-model — where a mutation is EXPECTED
arda.shm reports where a read differs from its germline. This fits the other half:
P(substitution | 5-mer germline context, region), estimated from what every mode already writes
under the default --shm framework.
arda shm-model -i mapped.airr.tsv -o shm.tsv --locus IGH⛔ Not a per-allele-per-position table, and the measurement is why. Across two donors a
per-allele-per-position profile transfers at Pearson r .2557–.5601 against 5-mer context's
r .7424–.7854; within one donor the positional table repeats at r ≈ .99, because it is a
portrait of that donor's expanded clones rather than a property of the allele. ⚠ Context is also
the only one of the two that reaches the V tail inside the junction — exactly the positions
arda.shm scopes out as unidentifiable, where a position has no estimable rate at all and a
5-mer takes its rate from every allele carrying it in clean framework sequence.
✅ Both parameter families are earned. AID's WRCY/RGYW motifs carry 4.66× / 5.00× / 4.77× /
4.69× the rate of every other covered position across four libraries, on an overall rate that
itself moves 1.6× between donors; and the CDR:FWR residual after conditioning on context
reproduces across donors — FWR1 0.630 / 0.626, CDR2 1.389 / 1.291.
⚠ The fitted scale is the sample's, not the model's, so it is written as provenance and never
applied.
Added: arda scenarios --shm-model — the consumer
scenarios.lattice bounded the templated V length with an exact common prefix, so on IGH one
substitution in the V tail forced the rest of that tail to be explained as insertion — and delV
and insVD are the two distributions the estimator fits. Those positions are now priced:
Π(1 − μ) over matches, μ / 3 over mismatches.
Measured on 26,619 real IGH junctions, with the model fitted on a different donor's library:
| parameter | exact | with --shm-model |
change |
|---|---|---|---|
delV mean |
4.148 | 2.708 | −1.440 |
delV P(no deletion) |
.1426 | .2388 | +.0962 |
insVD mean |
12.065 | 10.972 | −1.093 |
insDJ mean |
12.342 | 12.144 | −0.198 |
delJ mean |
12.247 | 12.253 | +0.006 |
The exact bound was charging 1.44 nt of V germline per IGH rearrangement to deletion, and
delV's mode moves from 1 to 0. ✅ The model reaches the V side only, so delJ at +0.006 against
delV at −1.440 is the result's own control. Cost: +4.8 % wall on a 3-iteration fit.
✅ The exact-match bound turns out to be the μ = 0 case of the new emission, not a separate
rule. shm=None remains the default and is byte-identical, so TR and unmutated IG are untouched.
Fixed: a V/J call naming a FAMILY with one functional gene now resolves
resolve_allele climbed exact → gene*01 → first allele of the gene, and none of those reaches a
call naming a family whose genes all carry a suffix: nothing is named TRBV20, only
TRBV20-1, so markup_cdr3 returned FailedBadSegment and dropped the V end entirely.
6,291 human beta chains in VDJdb carry such a call (TRBV20 1,839, TRBV3 1,039, TRBV24
291, TRBV29 70), and the germline places their CDR3s without a single edit once it resolves.
Never: the new rung refuses rather than guesses — the family must hold exactly one
functional gene. TRBV3 resolves because TRBV3-2 is a pseudogene; TRBV6, with five
functional genes, still returns "".
Fixed: a tie-group v_call silently found no germline
IGHV3-23*01 lives inside IGHV3-23*01,IGHV3-23D*01, so a lookup by a single allele found
nothing — it cost the SHM estimator reads and made the first SHM-aware lattice score exactly
like the exact one, which reads as a clean null rather than a missing key. 325 human IGH keys
become 349.
Full detail in CHANGELOG.md.