Somatic hypermutation#

A B-cell repertoire is not a set of germline calls. After activation, AID mutates the rearranged V region at ~10⁻³ per base per division, and the resulting lineage — who descends from whom — is the signal most BCR analyses are actually after. arda records that per read, in the one coordinate frame a downstream tool can use directly.

What arda records#

field

what it is

v_mutations / j_mutations

The read’s substitutions against its called germline: G45A,C112T — germline base, 1-based position in that segment’s own germline allele, read base. Empty when the read is germline over that segment.

v_identity

The same information as a fraction (identity over the V alignment).

v_cigar / j_cigar / c_cigar

Per-segment AIRR CIGARs. Indels live here, as I/D.

sequence_alignment / germline_alignment

The full aligned strings. Everything above is derived from them — see the warning below about deriving it yourself.

v_germline_start / v_germline_end, j_germline_*

The germline window the read actually covered. Mutation positions lie inside it.

The mutation lists come out of the same single walk of the alignment that builds the CIGARs, so they are close to free: +36 ms per 100,000 mapped reads (measured on 35,825 real bulk IG alignments, 26.9 → 39.7 ms, CIGARs byte-identical before and after) and +2.25 % on the output TSV (21,461,770 → 21,944,646 bytes). Amplicon wall clock does not move.

Why not just diff the two alignment strings#

The information was never missing — sequence_alignment and germline_alignment carry every column, and the germline arda reports matches the shipped allele on 28,365 of 28,365 mapped reads of a real bulk IG library (66,526 V mismatches, zero disagreements). What it was not, was usable.

Warning

arda aligns to a V + N-pad + J [+ C] scaffold, not to a germline. A consumer that does the obvious AIRR thing — diff the two alignment strings — gets 100,091 mismatches on that library of which 20,140 (20.1 %) are N-pad or constant-region columns. It attributes junction positions to a germline. Recovering it correctly needs arda’s scaffold geometry (mmseqs2_tstart, mmseqs2_t_vend, mmseqs2_t_jstart, mmseqs2_t_vjend); the two strings are also 30.7 % of the TSV.

Use v_mutations / j_mutations in preference to the alignment strings.

Why the junction is excluded#

Scoping SHM to the framework is not a convenience. The V germline’s 3’ tail and the J germline’s 5’ head both lie inside the junction, so counting mutations by segment boundary rather than by the junction boundary reports exonuclease chew-back and non-templated N/P bases as substitutions against a germline that does not template them.

The size of that effect was measured on a TRA amplicon (SRR5233636, 500,000 reads). T-cell receptors do not somatically hypermutate, so every entry on that library is spurious by construction:

  • 1.046 V and 1.658 J “mutations” per read under segment scoping;

  • 86.2 % of the J entries sat at J germline position ≤ 10 — the 5’ head, inside the junction;

  • splitting the load by the frequency of each (allele, position, alt) across reads carrying that allele: 13.0 % of J entries below frequency 0.01 (sequencing error), 80.8 % between 0.01 and 0.5 (junction diversity), 6.2 % at 0.5 or above (a genuinely wrong allele call).

Frequency alone does not separate those three populations; frequency and position relative to the anchor does. The high-frequency entries are junction-internal too — TRAV8-6*01 positions 281/282 at 0.88 against its anchor at 270, TRAJ8*01 position 1 at 0.67 against 26 — that is, allele differences in the templated V/J tail rather than mutation.

Since 2.16.0 v_mutations, j_mutations and v_identity are scoped to the framework using the per-read germline anchors, in every mode.

Important

The scoping change moves every SHM figure arda produced before 2.16.0. To reproduce one, use --shm both, which also emits the junction-inclusive values as v_identity_full, v_mutations_full and j_mutations_full. --shm off emits no SHM fields at all.

A read whose called allele has no usable anchor is left unscoped rather than emptied: arda does not know where that germline’s junction starts, and saying so is better than guessing. The raw anchors ship in the output, so the scoping stays checkable rather than merely trusted.

IGH, IGK and IGL are where SHM is real. The scoping runs on every locus anyway, because the effect is not IG-specific — it was measured on TRA. On the TR loci the entries that survive are allele mismatches in the templated framework, not hypermutation.

Rescoping a table written before 2.16.0#

arda shm rescopes an existing AIRR TSV. It needs no reference and no re-map — the anchors and the alignment strings are already in the file, and have shipped since 2.14.0:

arda shm -i mapped.airr.tsv -o rescoped.airr.tsv              # framework (default)
arda shm -i mapped.airr.tsv -o both.airr.tsv --mode both      # + the old values as *_full

Note

On a file older than 2.14.0, which has no v_anchor_nt column, this raises. It does not copy the input through with a success message — a command that reports success over an unchanged file is worse than one that fails.

Recomputing what you actually want#

arda applies the framework scope itself, but every read still carries the junction boundary in germline coordinates, so the split stays checkable — and re-derivable — from the TSV alone, with no reference:

v_anchor_nt

0-based offset of the Cys104 codon in the called V allele. A v_mutations entry at 1-based position p is junction-internal iff p > v_anchor_nt.

j_anchor_nt

0-based offset of the [FW]118 codon in the called J allele. A j_mutations entry at p is junction-internal iff p <= j_anchor_nt + 3 (the anchor codon is 3 nt, and the junction includes it — junction ≠ CDR3).

Framework-only V identity — which is what v_identity now is, and this is how to verify it:

fw_hi  = min(v_germline_end, v_anchor_nt)          # framework stops at Cys104
span   = max(0, fw_hi - v_germline_start + 1)
fw_mut = sum(1 for p in v_mutation_positions if p <= v_anchor_nt)
identity = 1 - fw_mut / span

Measured on the five records of examples/example.airr.tsv. The middle column is what arda emitted before 2.16.0; the one beside it is what it emits now:

v_call

v_identity (pre-2.16.0)

v_identity (2.16.0+, framework)

fw muts

junction muts

TRBV28*02

0.8723

1.0000

0

4

TRAV12-2*02

0.9964

1.0000

0

0

IGLV1-40*01

0.9833

0.9925

2

0

IGHV3-9*01

0.9262

0.9298

20

2

IGKV1-33*01

0.6920

0.6880

78

5

The first two rows are TR loci, where somatic hypermutation does not occur, so the junction-inclusive v_identity of 0.8723 for TRBV28*02 was measuring junction diversity outright. The bias is worst where the aligned span is mostly junction: that record covers 47 nt of germline, only 30 of it framework. The IG rows move much less, because their alignments are mostly framework.

Important

A mutation inside the junction is not attributable to any germline, and arda should not claim one.

V(D)J recombination trims the V, D and J ends by a variable amount and inserts non-templated N/P bases in their place, so the V-end / NDN / J-start partition of a junction frequently is not identifiable from the sequence at all. A “mutation” placed inside the NDN is a statement about a boundary that has no ground truth — see D segments and tandem D-D for the same rule on the D side. That is the design intent, and since 2.16.0 the implementation honours it.

Two more properties worth knowing:

  • Substitutions only. An SHM indel is one event of unbounded length, not a per-position call; it stays in the CIGAR as I/D. Germline coordinates on the far side of an indel are still correct, because the walk tracks the target position across gap columns.

  • Sequencing IUPAC bases become A at input, before SHM comparison. A column contributes only when both the normalized read and reference bases are ACGT; reference ambiguity is preserved.

  • Positions are on the coding strand. For rev_comp = T the sequence field is the read in its input orientation after normalization, so a read-side lookup (a Phred quality, say) must be made against the reverse complement of it.

Accuracy against IgBLAST#

3,000 hypermutated IGH reads (v_identity < .97), IgBLAST truth at v_score >= 70, per-mutation precision and recall over the whole mutation set:

arm

reads

precision

recall

identical mutation set

shares a V allele with IgBLAST

1,947

.98991

.98926

1,880 / 1,947 = .9656

+ germline frame verified

1,882

.99902

.99834

1,876 / 1,882 = .9968

The 65-read gap between the arms is a scoring artifact, not an error: IgBLAST answers with an allele tie list and aligns to one member of it, tied members differ in 5′ length, so a shared gene can be two different coordinate frames (33N vs 105N on IGHV4-30-2). The frame arm additionally requires IgBLAST’s own germline string at its own coordinate to be arda’s allele there. The 6 residual frame-verified disagreements are gap placement — 251 of the 3,000 reads carry a V indel.

Building a lineage tree#

A tree builder needs three things: a root (the germline), member sequences in one coordinate frame, and the clone membership that says which reads belong together. arda gives all three, but they come from two different stages.

arda map     --r1 R1.fq.gz --r2 R2.fq.gz -o mapped.airr.tsv
arda correct -i mapped.airr.tsv -o clones.tsv --read-map read_map.tsv

--read-map writes sequence_id -> corrected junction, which is the join key: it links each read’s v_mutations back to the error-corrected clonotype it was assigned to. Replay the mutations onto the clone’s common V germline window and you have a gapless alignment whose root is the germline — the input format dnapars, gctree and IQ-TREE read directly.

On a real bulk IG library that yields 66 IGH clones with ≥ 5 reads, 735 member sequences, 70 parsimony-informative sites, 0 malformed alignments, the largest clone 88 reads over 37 distinct mutation sets.

Warning

Clone membership requires a complete junction, and most hypermutated reads do not have one. On that library, 8,345 of 9,755 mutated IGH reads (85.5 %) carry SHM but join no clonotype — they are V-only fragments that never reach the Cys104 anchor, so --read-map cannot see them by construction. That is exactly the population with the most SHM signal per read.

If the question is the mutation spectrum (targeting, selection, replacement/silent ratios), score v_mutations over all mapped reads and do not go through the clonotype table. If the question is lineage, accept the coverage cost — a tree needs the junction to know what a clone is.

Two more cautions when consuming this:

  • An allele tie list shares one reference target named with the comma-joined list. Look the full v_call string up first; matching on the first member silently changes the coordinate frame.

  • Count over clonotypes, never reads, when inferring a donor’s germline variants from apparent mutations. One expanded clone otherwise makes its own SHM look germline.

Fitting a model of where mutation is expected#

Everything above says where a read differs. arda shm-model turns a pile of those reports into P(substitution | germline context, region) — a 5-mer context table, a per-region multiplier, and the sample’s own overall rate:

arda shm-model -i mapped.airr.tsv -o shm.tsv --locus IGH

It reads exactly what map / amplicon / rnaseq already write under the default --shm framework, and takes under a second on a bulk library (4,661 observations, 0.76 s). The output is one # provenance line, then kind / key / value rows.

Why context and a region multiplier, and not a per-allele-per-position table. Measured on four bulk IGH libraries from two donors, a per-allele-per-position profile transfers between donors at Pearson r .26–.56, against r .74–.79 for 5-mer context. Inside one donor the positional table repeats at r ≈ .99 — it is a portrait of that donor’s expanded clones, which is also why it must not be read as a property of the allele. The region term is earned separately: context carries about half the CDR-over-FWR contrast (4.65× raw, 2.2× residual) and the leftover reproduces across donors, with the shipped estimator giving FWR1 ×0.683 and ×0.675 on the two.

Important

The scale is the sample’s, not the model’s. Two donors measured here differ 1.6× in overall substitution rate (.031 against .050) with the same shape, so the fitted rate is written as provenance and a consumer rescales to the library in front of it.

Two properties of the table:

  • A region with no observed substitution is omitted, never given a multiplier of 0 — zero would read as “a substitution here is impossible” and silently zero every probability in that region. An omitted region means “no evidence”, i.e. the context term alone.

  • The denominator is capped at v_anchor_nt, because the numerator is framework-scoped. An uncapped count divides scoped mutations by unscoped coverage and understates the rate by however far the read runs into the junction.

Nothing in the annotation path reads this model yet, exactly as arda scenarios shipped before arda markup --d-prior existed to consume it. Fitting a model and adopting one are separate decisions; project/design-shm.md records what the consumers need next.

See also#

  • Usage — the full AIRR column list.

  • Error correction — a multi-base indel survives correction on purpose; it is the SHM signature, not an instrument error.

  • D segments and tandem D-D — the same non-identifiability rule, on the D side.