v2.20.0 — run QC, per-mutation quality, and one logger

arda stats — run QC as one long-format TSV

Written by every mode run as <prefix>.stats.tsv, and buildable from any subset of a run's
artifacts (--airr / --clones / --report / --r1/--r2), so it also works on a bare
arda annotate output. Reads only what already exists — no re-read of the FASTQ, no alignment.

Four columns, scope / key / metric / value, one value per cell:

scope	key	what
run	map / correct / assemble	the report verbatim: total and mapped reads, FASTQ bytes, read length, paired, threads, wall time, peak RSS, per-locus
sample	—	library totals, junction length and quality, SHM rate, V/J gene coverage
chain	TRB, IGH, …	per locus, reads and clonotypes: productive / non-functional, stop codons, out-of-frame, truncated junctions, min/max/mean junction length, junction quality, SHM rate, chimeras
v_gene / j_gene	TRBV19	reads and clonotypes per germline gene
allele_candidate	TRBV19*01:G45A	a recurrent, high-quality V mutation, with frequency and mean Phred

arda stats -i SAMPLE.airr.tsv -c SAMPLE.clones.tsv -r SAMPLE.arda.json -o SAMPLE.stats.tsv

• ⛔ Long, not wide. The metric set differs per scope, so a wide table is mostly empty cells.
• ⛔ A metric with no input is omitted, never emitted as 0. A run without --junction-quality
  must not read as "mean quality 0".
• ⛔ Truncation, a stop codon and an out-of-frame junction are counted separately.
• ⚠ allele_candidate is a shortlist, not a genotype call — arda does not genotype. Thresholds
  are exposed (--allele-min-frac, --allele-min-reads). Chimera / non-functional / stop-codon
  counts are flags, never filters.

Gene coverage is measured against the shipped cdr3_anchors.tsv — the germline set arda actually
maps to.

map --mutation-quality

v_mutation_quality / j_mutation_quality: the Phred of the read base behind each entry of
v_mutations / j_mutations, comma-joined, one-for-one and in the same order.

⛔ Driven by the emitted list, not by re-deriving one. _markup.segment_cigars finds a
superset of what the columns carry — arda.shm drops the junction-internal entries afterwards,
25 of 242 V rows on this repo's own fixture — so a re-derived list lines up in length only by
accident while pairing entry i with a different base's score.

⛔ The two quality columns use different encodings. junction_quality is raw Phred+33
characters; v_mutation_quality is comma-joined integers.

Logging, verbosity, progress

Global -v / -q / --log-file, before the subcommand, on one arda logger every module already
feeds. --log-file is always DEBUG whatever the console level is, with a timestamp and process
peak RSS on every line; -q does not silence it. map emits a throttled progress line — reads
seen, reads mapped, reads/s, peak RSS, at most one every 30 s.

⚠ Progress moved from stdout to stderr. A mode run now prints its output paths, one per line,
and nothing else on stdout — so $(arda map ...) and arda export-ref ... > out.tsv are safe.

Also

• map --report records paired / input_bytes / read_length_min|max|mean.
• Nextflow module bumped 2.16.0 → 2.20.0 and emits stats.
• New docs/qc.rst; README, docs/usage.rst, docs/api.rst, skills/arda/SKILL.md updated.

Full notes: CHANGELOG.md.

2.19.0 — correct --flag-chimeras

What's Changed

• 2.19.0 — correct --flag-chimeras by @mikessh in #101

Full Changelog: v2.18.1...v2.19.0

2.18.1 — the coverage assignment's inner scan, bounded in C++

What's Changed

• 2.18.1 — the coverage assignment's inner scan, bounded in C++ by @mikessh in #100

Full Changelog: v2.18.0...v2.18.1

2.18.0 — the assembler stops counting once the answer is known

What's Changed

• 2.18.0 — the assembler stops counting once the answer is known by @mikessh in #99

Full Changelog: v2.17.0...v2.18.0

2.17.0 — a truncated junction can be finished from the germline

What's Changed

• 2.17.0 — a truncated junction can be finished from the germline by @mikessh in #98

Full Changelog: v2.16.0...v2.17.0

v2.16.0 — three named modes, and SHM that stops at Cys104

⛔ Two breaking changes

1. arda rnaseq run is removed — the regime is the command name

It was the only pipeline entry point and was used for amplicon libraries as well as bulk, with the regime spelled out as four loose flags that do not compose. --two-pass, the one tuning flag it exposed for four releases, is a loss in both regimes (0.762× bulk, 0.87× IGH amplicon) — so the combination that was easiest to reach was the dominated one.

was	now
arda rnaseq run … --prefilter	arda rnaseq …
arda rnaseq run … --two-pass --fast-segments --v-only-on-segment	arda amplicon …
arda rnaseq run …	arda rnaseq … --exact
arda rnaseq map/correct/assemble	arda map/correct/assemble
arda rnaseq split/reduce/slurm	arda cluster split/reduce/submit
arda split/merge/slurm	arda cluster split-fasta/merge/submit-fasta

arda singlecell is reserved and exits non-zero — arda has no barcode or UMI concept, which is the real single-cell gap. Design in ROADMAP.md.

⚠ Mode defaults change behaviour: arda rnaseq enables --prefilter (~0.15 % of mapped reads, measured across 122 bulk datasets, up to 2.46 % on one), and each mode's --ec-mode defaults to its own regime. --exact --ec-mode fast reproduces the pre-2.16.0 output.

2. v_identity / v_mutations / j_mutations are scoped to the FRAMEWORK

This retracts a guarantee docs/shm.rst printed until 2.14.0 and moves every SHM number arda has published. The lists were scoped by segment, which is not the same as being outside the junction: a rearranged junction is V 3′ tail + N/P + J 5′ head, so both germlines' templated tails lie inside it and exonuclease chew-back was reported as substitution against a germline that does not template it.

Measured on a TRA amplicon (500 k reads) — TCRs do not hypermutate, so every entry was spurious by construction: 1.046 V and 1.658 J entries per read, 86.2 % of the J ones at germline position ≤ 10. On the committed example, TRBV28*02's v_identity was 0.8723 and is now 1.0000.

• --shm framework (default) · --shm both (also emits the old values as *_full) · --shm off
• arda shm -i in.airr.tsv -o out.airr.tsv rescopes an existing table with no reference and no re-map; it raises on a pre-2.14.0 file rather than copying it through with a success message.

Added

--call-level {allele,gene} — collapses allele-level call splits (Jurkat: TRGJ1*01 64 reads vs TRGJ1*02 140 on one junction, invisible to every error model because identical junctions have no discriminating base) and allele-only tie lists. --no-isotype.

Fixed — the install path

Found by running setup.sh end to end: arda.__version__ and pyproject.toml were two independent literals and they disagreed; --tests ran a pytest that was never installed and swallowed the failure; nothing verified the C++ extensions built; a stale build/ broke every rebuild; pybind11 was unbounded; docs/installation.rst still described a conda setup.sh.

Nextflow

params.regime selects the arda command. New arda_shm / arda_call_level. Pinned 2.16.0 — a hard minimum in both directions.

Full changelog: CHANGELOG.md

2.15.0 — the contig assembly audit

Fixed — contig assembly (⛔ read counts and contig sequences change)

A second audit went over SHM calling and contig assembly. Eight defects, each with a regression test
verified to fail without it.

⛔ Contigs were READ-ORDER DEPENDENT. The extension tie-break was len(ext) > len(best_ext) with
a strict >, so equal-length candidates were resolved by the order the posting list happened to be
in — which is AIRR row order, which comes from a threaded MMseqs2 search. The contig sequence, and
every junction derived from it, could differ between runs on the same input. Now ordered on
(length, sequence), a total order, in both the 3′ and 5′ passes.

A contig's junction is now only attributed to members that COVERED it. Membership is granted on
a min_overlap match, and once the extension passes have accumulated germline at the contig ends
that match can be pure germline V — the 5′ pass says so itself: "that region is shared germline, so
any V-read of the gene extends it correctly". Every incomplete member was nonetheless stamped with
the contig's junction, so a read of a different clone of the same V gene was credited to this
clonotype's duplicate_count on no clone-specific evidence. Members now carry their span in contig
coordinates and must cover the junction by ≥10 nt.

⚠ This reduces reported abundance. On the committed 1,035-read example: rescued members 162 →
145, reads counted 384 → 367, clonotypes 19 → 18. Those 17 reads are genuinely uncounted, not
moved — a germline-only read matches no junction, so coverage assignment cannot re-place it. It is a
deliberate precision-over-recall trade, unlike the error-correction path where read conservation is
an invariant. ✅ The clonotype that went is the case for the change: CAASMAGGGNKLTF under
TRAV13-1*01 with 2 reads against the same junction under TRAV13-1*02 with 28 — two alleles of one
gene, i.e. de-duplication rather than loss.

A rejected contig now releases its reads. used was set as reads were recruited but a contig
dropped for having <2 members never gave them back, so a seed that failed to extend was permanently
consumed and could not join a later contig even as an ordinary extension member. Seeds are tried
longest-CDR3-tail first, so this stranded exactly the short-tailed reads that most need a contig.

The assembly k-mer index is bounded at insert. It posted every k-mer position of every mapped
read of the locus with no cap; scan_cap bounded only how many postings were consumed, and
--assemble is on by default. ⚠ Proved equivalent before shipping, since a cap that changes the
candidate set is a behaviour change and not an optimisation: both consumers already slice
[:scan_cap]. Measured on 20,000 reads sharing a 60 nt germline prefix — candidate sets identical
for every k-mer, postings 1,600,000 → 784,000, ~112 MB → ~82 MB.

Fixed — clonotype reporting

Isotype is one vote per FRAGMENT, at last. _dominant_ccall deduplicated its read list down to
fragments and then re-expanded to one entry per row, so a fragment whose two mates both carried a
c_class voted twice, and an assembly-rescued fragment voted again. A one-fragment minority could
outvote a two-fragment majority. ⛔ The tally was also order-dependent — Counter.most_common(1)
breaks ties by insertion order — so a two-way isotype tie could report a different class run to run.
Now lexicographic.

Under --all-junctions, an assembled row outranks the read's own truncated junction, and a
clonotype left with no reads is not emitted. Coverage pass 1 walks the concatenated frame mapped-rows
first, so the read's own truncated junction won the race and the contig's clonotype was emitted with
duplicate_count 0 — measured on a 2-read fixture: CAKGALQ dup=2 beside CAKGALQKW dup=0. Neither
half touches read conservation: the reads move, and a dropped row has none by definition.

reads_with_junction no longer double-counts assembly-rescued reads (once for the incomplete
Stage-1 row, once for the assembled row) while being documented as a Stage-1 statistic.
reads_from_assembly and empty_clonotypes report the two new quantities separately.

Fixed — arda resolve-ties degraded silently without a reference

It caught OSError per locus, continued with an empty germline set, and returned output
byte-identical to its input with expanded: 0 — indistinguishable from a library that genuinely had
no ties. The same "raise, never degrade" rule --min-junction-q and the quality rescue already
follow. ⚠ This had kept CI red since tie lists were added.

Added — the junction boundary in GERMLINE coordinates

v_anchor_nt and j_anchor_nt (0-based Cys104 / [FW]118 offsets in the called allele) are emitted
per read, appended last. They existed only inside the reference before, so a consumer holding the TSV
could not tell a framework mutation from a junction one. Now:

• a v_mutations entry at 1-based p is junction-internal iff p > v_anchor_nt;
• a j_mutations entry at p is junction-internal iff p <= j_anchor_nt + 3.

This is what makes allele re-assignment and SHM correction a downstream job — see
docs/shm.rst for the framework-only identity recipe and the measured split.

Changed

scripts/pack_reference.sh rejects an empty reference file, not just a missing one. A
build-db run without IgBLAST emits 0-byte alleles.*/markup.*, which passed the presence check,
shipped, and crashed every fresh install with NoDataError: empty CSV — v2.5.7 shipped exactly that.
The smallest genuine file is 24,372 B against a 1,000 B floor (#83).

Correction to the 2.14.0 notes: the coverage k-mer cap's under-assignment is library-dependent,
not a general property. It tracks how much germline the junctions share: 1.6 % on a TRA amplicon
(short, near-germline junctions over 36,741 roots) against 0.013 % on a TRB amplicon of comparable
depth (72,339 roots, D+N-bearing junctions that are far more specific).

Known defects, documented not fixed

v_mutations / j_mutations and v_identity include junction-internal positions — the V
germline's 3′ tail and the J germline's 5′ head lie inside the junction, and both are scoped by
segment rather than by the junction boundary. Measured on a TRA amplicon, where TCRs cannot
hypermutate so every entry is spurious by construction: 1.046 V and 1.658 J entries per read, and
v_identity reporting 0.8723 on a TCR whose framework-only identity is 1.0000. ⚠ Frequency
alone does not separate error from allele from junction diversity — frequency and position against
the anchor does. Fixing it inside arda changes every published SHM number, so it gets its own round;
docs/shm.rst carries the recipe and the measurements meanwhile.

2.14.0 — the coverage read leak, fixed and measured

Fixed — a full-depth read leak in coverage assignment (⛔ read counts change)

--ec-mode accurate|amplicon|rnaseq lost reads out of the clonotype table, and the amount
scaled with library depth, so no local test saw it. Measured across the 16-sample golden set at full
depth: 7 samples leaked, −107,440 reads in total, worst on the TRA amplicons (SRR5233636
−25,468 of 1,838,213, i.e. 1.39 %; SRR5233635 −7,383) and present on both MIGEC published-truth
libraries (−755 IGH, −1,351 TCR). The cell lines all gained (Jurkat +44, Raji +46,
MouseSpleen +5), which is exactly why a monoclonal QC panel never caught it.

The cause is not the error model. _assign_coverage bounds each k-mer's postings at cap and fills
them in descending abundance; an alias — a junction the quality gate vacated, kept in the index
only so partial reads that covered it still reach the parent — is ordered by its parent's count,
which is high by construction. Aliases therefore sorted to the front and evicted genuine
low-abundance roots, and every partial read whose only home was such a root went unassigned. On
SRR5233636, every arm emits an identical 36,587-row clonotype table, so this is purely read
assignment:

index	reads assigned	vs fast
cap 64, aliases ordered by abundance	1,812,740	−25,473
cap 64, aliases off	1,838,181	−32
cap 64, roots first (this release)	1,841,624	+3,411
cap 1024, aliases by abundance	1,869,556	+31,343

The alias mechanism was added to rescue 5 reads of 9,208 on Ramos and was costing 25,441 here. It is
worth keeping — the same aliases gain 31,343 once the index can hold both — it simply must not
outrank the roots. Every mode now gains rather than loses: on SRR5233636 accurate +3,411,
amplicon +2,827, rnaseq +3,393.

CorrectReport.reads_assigned is new and is the reason this was invisible. reads is the
spanning-read count taken before correction, so it cannot move; comparing it across --ec-mode
returns 0 on every sample and reads exactly like conservation holding. reads_assigned is
sum(duplicate_count) over the emitted table — the quantity the invariant is actually defined on.

The quality rescue no longer merges across loci. It grouped candidate parents by junction
length alone, and --ec-mode amplicon opens the radius to 12 substitutions, so on SRR5233636
3 of 9,025 rescues were 1-read TRB clonotypes absorbed into abundant TRA clonotypes at 11–12
substitutions — misattributed expression, not lost reads. The search now runs once per locus, which
also bounds it: it is quadratic within a length bin and single-threaded, and a 397 k-clonotype
library was still running after 39 minutes of CPU against ~5 for a 49 k one.

⚠ V and J stay ignored by the rescue, deliberately. The abundance model defaults to
--require-vj because a true error keeps the germline call — true for the 1–3 substitution
neighbours it collapses. The rescue targets the opposite class, where the whole junction window is
unreliable (median mean Phred 16.5–20.1), and a read that bad has an unreliable V/J call for the
same reason. 50.9 % of amplicon rescues cross V. What protects a genuine clone is the two gates
that are trustworthy: its reads must be measurably bad, and the parent must be
lowq_min_ratio times more abundant.

The quality gate's parent must now be able to have PRODUCED the read. Its only test was
count(parent) > count(child) — one extra read made anything within --max-subs a parent, and at
3 substitutions nothing supports that. The plausibility test is computed from this read's own
Phred: the product of 10^(−Q/10) over the discriminating positions, times the parent's count, must
reach the child's. ⚠ Using the global --error-rate instead was tried and is wrong the other way —
it makes the gate a strict subset of the abundance model, which is what the gate exists to reach
past (MIGEC at 1e-5 went 1,633 → 1,633 error clonotypes, from 158).

--ec-mode amplicon|rnaseq now raises without junction_quality instead of silently skipping
the rescue and reporting counters that are indistinguishable from a clean library — the same
"raise, never degrade" rule --min-junction-q already followed.

The Nextflow module emitted a flag arda rnaseq run did not have. main.nf appended
--clonotype-key to the run command line, where only correct accepted it, so
arda_clonotype_key = 'junction' failed the process outright. run now accepts it and plumbs it
through to correct.

Changed

Orphon V/J genes are excluded from the reference. IMGT /OR genes sit outside their locus
(TRBV20/OR9-2 is on chromosome 9, TRB is on 7) and cannot rearrange, so they were pure false-call
surface. ⚠ This affects arda build-db; an auto-fetched prebuilt reference needs regenerating to
pick it up.

D germline loading is cached. _d_germlines re-opened and re-parsed d_germlines.fasta on
every call and _clonotype_d calls it once per clonotype — 54.4 µs each, so 2.0 s wasted on a
36,741-clonotype amplicon and 21.6 s on a 397,305-clonotype library, on the default --map-d
path. Now 0.056 µs.

Added

Tie lists (--tie-lists, off by default) and arda resolve-ties. A read aligned over a span
that several germlines share carries no base that separates them, so naming one is a claim the data
does not support. Membership is decided per read from the span it already aligned — a string
comparison against the reference, not a new alignment — then ranked library-wide so the allele the
whole library supports leads, with only unambiguous reads voting.

⛔ Known defect, newly measured — v_mutations / j_mutations include junction positions

docs/shm.rst claimed the mutation lists were scoped to V and J structurally, so "a junction
position has no germline coordinate to be filed under and cannot enter the list by any code path."
That claim is false and is retracted in this release. The N-pad is excluded, but the V
germline's 3' tail and the J germline's 5' head lie inside the junction, and mutations are scoped
by segment (t <= t_vend, [t_jstart, t_vjend]) rather than by the junction boundary — so
exonuclease chew-back and non-templated N/P bases are reported as substitutions against a germline
that does not template them.

Measured on a TRA amplicon (SRR5233636, 500,000 reads), where T-cell receptors do not
hypermutate so every entry is spurious by construction: 1.046 V and 1.658 J mutations per read,
with 86.2 % of J entries at J germline position <= 10. Splitting the load by the frequency of
each (allele, position, alt) across reads carrying that allele separates three superimposed
populations:

frequency	share of J entries	share of V entries	what it is
< 0.01	13.0 %	39.0 %	sequencing error
0.01–0.5	80.8 %	59.2 %	junction diversity misattributed as SHM
>= 0.5	6.2 %	1.9 %	a genuinely wrong allele call

The high-frequency tail is real and separate — TRAV8-6*01 positions 281/282 at 0.88, TRAJ8*01
position 1 at 0.67 — i.e. the called allele is not the one the reads carry. v_identity has the
same scope defect: it runs to t_vend, so it is depressed by junction diversity rather than by
mutation load.

Not fixed in this release. The fix needs the scaffold's CDR3 boundary threaded into the C++
markup and the Python reference implementation together, and it changes every published SHM number,
so it gets its own round with its own measurement. Documented in docs/shm.rst with the workaround
in the meantime: separate the three populations by frequency.

Known limitation, measured

The coverage k-mer cap under-assigns ~1.6 % of reads in the default path, independent of the alias
fix. On SRR5233636 with --ec-mode fast (no aliases exist), the clonotype table is identical at
every cap:

cap	reads assigned	wall
64	1,838,213	207.4 s
128	1,844,359	312.1 s
256	1,850,917	485.9 s
512	1,867,867	740.7 s
1024	1,867,904	865.0 s

It saturates at 512: the whole prize is +29,654 reads and 1024 adds 37 more. Peak RSS is flat, so it
is CPU in the alignment inner loop. The default stays at 64 — 3.6× on this stage for 1.6 % of
coverage-based abundance is a per-workload judgement. ⛔ Guaranteeing every root a posting does
not help: built, measured, byte-identical on an amplicon and a bulk sample at two caps. No root
is ever fully unreachable; the reads are lost because a root's surviving postings sit at k-mer
positions the partial read does not cover.

arda 2.13.2 — Nextflow + cluster integration

Added — the denoising framework is reachable from Nextflow, and the cluster path is documented

Audit of "is everything shipped and usable" found two real gaps. Both are integration, not
behaviour: no arda output changes in this release.

The Nextflow module could not reach the framework. main.nf, nextflow.config and meta.yml
had zero mentions of --ec-mode, --clonotype-key or --junction-quality, so a pipeline user had
no way to it short of task.ext.args surgery. Now:

params {
    arda_ec_mode       = 'amplicon'   // fast (default) | accurate | amplicon | rnaseq
    arda_clonotype_key = 'junction'   // full (default) | junction
}

Both are validated against their allowed sets, and the module warns when arda_ec_mode and
regime disagree — the amplicon and rnaseq presets are tuned to opposite clonotype-size
distributions, so the mismatch is a real cost rather than a no-op. Declared in nextflow.config
with the measurement behind each, and read through params.getOrDefault so the module stays
correct when included without its config.

SLURM was implemented but undocumented. arda slurm / split / merge and
arda.cluster.render_submit_script have shipped for releases with no page describing them. New
docs/cluster.rst: the one-command chain, why split and split_pairs are not interchangeable
(one writes FASTA and round-robins records, splitting a fragment's mates across shards), why
Stage 2 must run once globally rather than per shard, aligner pinning, regime choice, and
resource sizing — including that prefilter threads saturate at 16 and regress at 32.

New docs/use_cases.rst: bulk RNA-seq, amplicon, monoclonal QC, negative controls, low-frequency
variants, SHM, and the rules for comparing arda against another tool (name the stage; benchmark
every tool at its best config; compare at gene level; give each call metric its own denominator).
It also states plainly that a 200 k subsample of a bulk library yields 1–3 clonotypes and nothing
computed on it means anything.

README gains the monoclonal-QC and per-regime rows, the read-conservation invariant, and why the
modes are off by default.

arda 2.13.1 — a Q1 base broke AIRR parsing

Fixed — a Q1 base in junction_quality made the whole AIRR unreadable

junction_quality is a Phred+33 string, and chr 34 is ", i.e. Q1 — a legitimate score that
any low-quality base produces. polars' CSV reader treats it as a quote character, so one such
base collapsed the parse of the entire file:

ComputeError: CSV malformed: expected 1 rows, actual 155 rows

Found on a real Raji run in the round-23 benchmark: exactly one row of the file contained a ",
and every correct leg for that sample died. It is not rare — it needs one Q1 base anywhere in one
junction, so any library with a low-quality tail hits it.

Underneath was a second problem: arda wrote the format two ways. The streaming writer
(_markup.format_rows, which produces every map output) emits raw fields and a truly empty
string for a missing value; polars' write_csv quotes, rendering an empty string as the two
characters "". Reading unquoted is right for the big files and would have turned every empty
field of an older polars-written one into a literal "".

So both sides are fixed: every AIRR reader now reads unquoted and normalises a literal "" back
to empty (an AIRR field is never legitimately that two-character string), and arda's polars
writers emit quote_style="never" so there is one dialect going forward.

⚠ examples/rnaseq/clones.tsv changes accordingly — 18 rows where "" becomes genuinely empty.
No value changes.

Pinned by test_a_q1_base_in_the_quality_string_does_not_break_the_parse.