From d3cfc2aa6ae4e0cb3e3b78fcd97cbc93567d7a7a Mon Sep 17 00:00:00 2001 From: mikessh Date: Sat, 18 Jul 2026 22:22:52 +0300 Subject: [PATCH] Relocate memory/ and benchmarks/ out of the tool repo MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The design-rationale notes (`memory/`) and measured benchmark results (`benchmarks/RESULTS.md`) are dev/analysis artifacts, not part of the shipped package (neither was in pyproject/MANIFEST, so the wheel is unchanged). They now live in the separate arda-benchmark analysis workspace. Drops the README links that pointed at the moved dirs so nothing 404s. Two in-code comment breadcrumbs ("see memory/mmseqs-params.md", "see memory/markup-transfer.md") are left as-is — internal prose, not links. Co-Authored-By: Claude Opus 4.8 --- README.md | 10 +-- benchmarks/RESULTS.md | 81 -------------------- memory/README.md | 14 ---- memory/cluster-slurm.md | 22 ------ memory/d-mapping.md | 121 ------------------------------ memory/discordance-and-scaling.md | 84 --------------------- memory/igblast-gotchas.md | 29 ------- memory/junction-markup.md | 111 --------------------------- memory/markup-transfer.md | 79 ------------------- memory/mmseqs-params.md | 43 ----------- memory/reading-frames.md | 24 ------ memory/scaffold-enumeration.md | 28 ------- 12 files changed, 3 insertions(+), 643 deletions(-) delete mode 100644 benchmarks/RESULTS.md delete mode 100644 memory/README.md delete mode 100644 memory/cluster-slurm.md delete mode 100644 memory/d-mapping.md delete mode 100644 memory/discordance-and-scaling.md delete mode 100644 memory/igblast-gotchas.md delete mode 100644 memory/junction-markup.md delete mode 100644 memory/markup-transfer.md delete mode 100644 memory/mmseqs-params.md delete mode 100644 memory/reading-frames.md delete mode 100644 memory/scaffold-enumeration.md diff --git a/README.md b/README.md index bd60488..21d3c20 100644 --- a/README.md +++ b/README.md @@ -100,8 +100,7 @@ covering every junction-repair outcome, including one arda reports and refuses t it refuses outright; and a 1,035-read FASTQ that runs the whole bulk RNA-seq pipeline in ~6 s. Two tests re-run that script and fail if a committed artifact stops reproducing. -See [`CHANGELOG.md`](CHANGELOG.md) for what changed per release and -[`benchmarks/RESULTS.md`](benchmarks/RESULTS.md) for measured speed and accuracy. +See [`CHANGELOG.md`](CHANGELOG.md) for what changed per release. The reference database ships with **precompiled MMseqs2 indexes** (`database/vdj//mmseqs/`), so annotation runs out of the box with no @@ -258,7 +257,7 @@ gold-standard reference for benchmarking (see the `arda-benchmark` project). its errors located and conservatively repaired, and optionally given a D gene inferred from the junction *length* (`--d-posterior`). -See [`memory/`](memory/) for design rationale and gotchas. Fast sequence +Fast sequence primitives (`translate`, `detect_coding_frame`, `reverse_complement`, `back_translate`) live in the C++ extension and are re-exported from `arda.refbuild.translate` — mirpy-API-compatible, so mirpy can `import arda` and @@ -330,7 +329,4 @@ honest, `.[test]` for `airr` schema validation. Without them those tests **skip* `pip install -e '.[test]'` before reading a green suite as full coverage. Layout: `src/arda/{refbuild,annotate}`, C++ in `src/_markup/markup.cpp`, -references in `database/`, downloads in gitignored `bin/` + `data/`. Design rationale and -the gotchas that cost us the most live in [`memory/`](memory/) — read -[`memory/d-mapping.md`](memory/d-mapping.md) and -[`memory/junction-markup.md`](memory/junction-markup.md) before touching either. \ No newline at end of file +references in `database/`, downloads in gitignored `bin/` + `data/`. \ No newline at end of file diff --git a/benchmarks/RESULTS.md b/benchmarks/RESULTS.md deleted file mode 100644 index b4bba82..0000000 --- a/benchmarks/RESULTS.md +++ /dev/null @@ -1,81 +0,0 @@ -# arda — achieved performance & accuracy - -Measured on Apple M3 (16 threads), MMseqs2 18.x, IgBLAST 1.22.0, against the -committed reference DB and test fixtures. All numbers are reproducible with the -scripts/tests noted; this file is the on-record snapshot. - -## Accuracy vs IgBLAST (gold standard) - -Offline, on ~7.3k real GenBank mRNA across all five organisms -(`tests/assets/realworld/`, compared on IgBLAST-productive records). - -| organism | region concordance (productive) | -|---|---| -| human | 98.6% | -| mouse | 99.6% | -| rat | 98.0% | -| rabbit | 99.7% | -| rhesus_monkey | 99.4% | - -- **V gene** assignment agrees ~**100%**. -- **junction_aa / cdr3_aa** match IgBLAST ~**99–100%** on productive-canonical - records, and satisfy the AIRR invariant `cdr3_aa == junction_aa[1:-1]` exactly - for *every* emitted junction (including out-of-frame ones, rendered with a `_`). -- **D gene** (where both tools call a D): TRB/TRD ~**97%** gene agreement; IGH - ~**46–69%** — IGH D is inherently ambiguous (≈50 paralogous germlines + SHM), - consistent with inter-tool reports. -- Remaining region diffs are one-residue FR1/CDR3 boundary conventions, not errors. - -Reproduce: `pytest tests/realworld -s` (offline); raw vs-IgBLAST run on synthetic -IGH: `scripts/bench_vs_igblast.py` (98.9% region concordance at n=10k). - -## Speed vs IgBLAST - -Synthetic human IGH, 16 threads (`scripts/bench_vs_igblast.py`): - -| sequences | arda | arda rate | speedup vs IgBLAST | -|---:|---:|---:|---:| -| 10,000 | 5.5 s | ~1.8k/s | **4.4×** | -| 50,000 | 16 s | ~3.0k/s | **7.3×** | -| 100,000 | 30 s | ~3.3k/s | **7.9×** | - -IgBLAST runs at ~0.44k seq/s on the same data. arda's rate rises with input size -as fixed costs (DB load) amortize; for small batches that overhead dominates. - -## Bulk RNA-seq (prefilter) - -mmseqs k-mer prefilter rejects non-receptor reads before alignment, so realistic -bulk RNA-seq (~1–5% receptor) is several-fold faster than amplicon. 150 nt reads, -16 threads (`scripts/bench_prefilter.py`): - -| receptor content | throughput | -|---:|---:| -| 100% (amplicon) | ~5.7k reads/s | -| 10% | ~19k reads/s | -| 1% (blood RNA-seq) | ~25k reads/s | - -**32-core / 30M-read estimate:** ~10–20 min for a ~1%-receptor library — same order -of magnitude as a STAR genome-alignment pass (STAR is faster per read but indexes -the whole genome). Memory is flat (bounded-chunk streaming); multi-node via -`arda slurm`. - -## D-mapping overhead - -D mapping is a short per-hit gapless local alignment (C++ `d_local_align`), not a -second mmseqs pass. On a VDJ-only query set (worst case), 16 threads -(`tests/benchmark/test_d_mapping_overhead`): **~±1%** end-to-end, **~7 µs/seq**. -On by default; `--no-map-d` to disable. - -## How to reproduce - -```bash -# accuracy (offline, committed fixtures) -ARDA_MMSEQS=$(which mmseqs) pytest tests/realworld -s - -# guarded benchmarks (timing / memory / scaling / D-overhead) -RUN_BENCHMARK=1 ARDA_MMSEQS=$(which mmseqs) pytest tests/benchmark -s - -# headline speed + concordance vs IgBLAST -ARDA_MMSEQS=$(which mmseqs) python scripts/bench_vs_igblast.py -ARDA_MMSEQS=$(which mmseqs) python scripts/bench_prefilter.py -``` diff --git a/memory/README.md b/memory/README.md deleted file mode 100644 index 92a9ea0..0000000 --- a/memory/README.md +++ /dev/null @@ -1,14 +0,0 @@ -# arda design memory - -Running notes on non-obvious design decisions and gotchas, so future work -(and future sessions) don't re-derive them. One topic per file. - -- [scaffold-enumeration.md](scaffold-enumeration.md) — why V×J only (not V×D×J), dedup, N-spacer. -- [reading-frames.md](reading-frames.md) — V frame detection, J frame from aux, FR4 sanity check. -- [igblast-gotchas.md](igblast-gotchas.md) — dummy D db, IGDATA, aux files, AIRR coords. -- [markup-transfer.md](markup-transfer.md) — runtime projection, indel/strand semantics, D-segment TODO. -- [mmseqs-params.md](mmseqs-params.md) — tuned nt/aa params, all-loci single DB, caching, speedup. -- [discordance-and-scaling.md](discordance-and-scaling.md) — why arda≈IgBLAST, alignment-phase frame fix, bulk RNA-seq prefilter speed, 30M-read estimate. -- [cluster-slurm.md](cluster-slurm.md) — split/merge/slurm sharding design + how it's tested without a cluster. -- [d-mapping.md](d-mapping.md) — anchored interior, one E-value not four floors, TRBD2×TRBJ1 is impossible, aa D, D-D, dead ends. -- [junction-markup.md](junction-markup.md) — junction≠CDR3, per-allele anchors, detect-vs-repair, the aligner tie-break bug. diff --git a/memory/cluster-slurm.md b/memory/cluster-slurm.md deleted file mode 100644 index 52604a5..0000000 --- a/memory/cluster-slurm.md +++ /dev/null @@ -1,22 +0,0 @@ -# Multi-node (SLURM) sharding - -`src/arda/cluster.py` + CLI `split` / `merge` / `slurm`. - -- **split(input, out_dir, shards)**: one streaming pass, round-robin (record k → - shard k%shards) into `shard_.fasta`. Round-robin balances load and is the - reason shard sizes differ by ≤1. Writes FASTA (quality dropped; arda ignores it). -- **merge(dir|list, output)**: concatenate per-shard AIRR TSVs, header once. -- **render_submit_script(...)**: emits a `submit.sh` that runs `arda split`, then - `sbatch --array=0-(N-1) --wrap 'arda annotate ... shard_${SLURM_ARRAY_TASK_ID}'`, - then `sbatch --dependency=afterok:$ARRAY_JID --wrap 'arda merge'`. `ARDA_MMSEQS` - is exported into the script if set in the caller's env. -- CLI `arda slurm -i ... -o ... --shards N [--submit]` writes `arda_slurm/submit.sh` - (chmod +x) and optionally runs it. - -Design choice: pre-split once (low I/O) rather than each array task reading the -whole file (round-robin in `read_sequences` would mean N× read I/O on a huge file). - -Tested without a cluster: `tests/unit/test_cluster.py` (split partitions all records -disjointly + balanced; merge keeps one header; script chains the three steps with -the dependency). A live SLURM run is still TODO (user will provide a cluster). -Single-node path is unchanged (streaming `annotate_file`). diff --git a/memory/d-mapping.md b/memory/d-mapping.md deleted file mode 100644 index 689abe0..0000000 --- a/memory/d-mapping.md +++ /dev/null @@ -1,121 +0,0 @@ -# D-segment mapping: interior, gate, genomic order, alphabets - -Everything here was measured. Don't re-derive it; the ROADMAP records the dead ends too. - -## 1. The interior must come from the anchors, not the projection - -A scaffold is `V + 9 nt N-pad + J`. A real read has a 20–40 nt N-D-N region there, and -mmseqs — unable to align anything to a run of N — parks those bases against the flanking -V and J. Both projected boundaries then march inward. On the human fixtures the projected -`v_sequence_end` sits +6 nt too far and `j_sequence_start` −13 nt too early (IGH, medians), -leaving an 11 nt "interior" where the truth is 37 nt: **a window too small to hold an IGH D -at all.** - -`transfer._anchored_vj_bounds` fixes this from the per-allele germlines in -`cdr3_anchors.tsv` — the V templates the junction's 5′ end, the J its 3′ end, so the bounds -are the longest common prefix / suffix. vs IgBLAST: `v_sequence_end` within 2 nt for 85 % of -IGH (projection: 43 %), `j_sequence_start` for 99 % of TRB (projection: 49 %). - -**This — not "paralogous germlines + SHM" — was the cause of the old 46–69 % IGH D -concordance.** ROADMAP said the latter for months. It was wrong. - -SHM truncates the exact match early, which *widens* the interior. That is the safe -direction: a D is never clipped, only surrounded by a little more sequence. - -## 2. One E-value, not four score floors - -Four hand-tuned per-locus floors (9/7/7/6 for huIGH/huTRB/huTRD/moTRB) became one knob: - - E = K·m·n·e^(−λS) → S_min = ceil( ln(K·m·n/E) / λ ) - -with `λ = ln((1−p)/p)` for `p` = the chance two residues match. **p = 1/4 recovers λ = ln 3** -for nucleotides. For amino acids p is *not* 1/20: N-region inserts and D germlines are both -G/S/Y-rich, and the measured p over real middles × real D frames is **0.0613 → λ = 2.7285**. - -The absolute E is not calibrated (inserts are Markov, not uniform; K is borrowed). Treat it -as an m·n-corrected score and the threshold as an operating point. nt: `E ≤ 0.2`. aa: `E ≤ -0.05` plus a floor of 4, because a 22–38 residue aa database leaves the E-value badly -under-calibrated at small n (at 0.2 the IGH false-call rate reaches 13 %). - -`d_support` ships the E-value so a consumer can re-threshold. - -## 3. Genomic order forbids TRBD2 × TRBJ1 - -The TRB locus runs `TRBD1 – TRBJ1 cluster – TRBC1 – TRBD2 – TRBJ2 cluster – TRBC2`, and -V(D)J joining *deletes* the intervening DNA. TRBD2 therefore cannot reach any TRBJ1. This is -genomic order, **not** a usage preference, and it holds in human, mouse, rat and rhesus. -IGH and TRD place every D 5′ of every J: nothing is forbidden there (verified — 0 structural -zeros, and IGH's 29 sub-1e-3 cells are low-usage genes, not a block). - -Unenforced, TRBD2 (16 nt) simply outscores TRBD1 (12 nt) on noise: **17 % of real human TRB -J1-cluster D calls (18/104) and 12 % of mouse (7/59) were impossible**, at E-values (median -0.096) sitting in the chance band — against 0.014 for the producible TRBJ2 × TRBD2. - -Two traps found while fixing it: - -* **OLGA's human TRB model does not encode the constraint, and hides that behind an allele.** - `TRBD2*01` correctly falls to ~1e-5 on J1 rows, but `TRBD2*02` absorbs 21–27 % of every J1 - row. At *gene* level — which is what `dpost` consumes — the unmasked model claims - `P(TRBD2 | TRBJ1) ≈ 0.23`. The mouse model does learn it (~0.003). Mask, don't trust. -* **`dpost` backs off to the marginal when the J allele is outside the model, and the human - model has no `TRBJ1-6*01` row at all** — so the mask must be re-applied after the backoff - (`_mask_forbidden`), or TRBD2 walks back in through the door. - -**Mask-only, never mask-and-rescore.** `_D_MAX_EVALUE` is an operating point calibrated -against the *full* locus D set, so shrinking `n` to the masked set silently loosens the gate: -measured, it admitted 58 new calls at median E = 0.098, squarely in the chance band. Holding -`n` fixed, a J1 record can lose an impossible call but never gain a weak one. Result: 18 -human impossible calls gone, 3 recovering as TRBD1 at E ≤ 0.054, the 86 good D1 calls -untouched. - -Where the accuracy actually landed: **all of it in the nt caller.** The aa posterior changed -*zero* calls on real data — the unmasked prior already favoured TRBD1 3.3:1 and the aa -evidence agreed. What masking buys `dpost` is honest confidence (posterior 1.0, entropy 0, -instead of 0.92/0.40 for a certainty) and an unrepresentable rather than merely unlikely -wrong answer. - -## 4. aa input: three frames, one allele - -A trimmed D has no knowable reading frame, so the aa reference carries all three translations -of every allele (`reference._load_d_germlines_aa`). Two frames of one allele can tie on one -span — `_best_d` de-duplicates, because that is one allele, not an ambiguity. - -At E ≤ 0.05, floor 4, against a composition-preserving shuffled null (the "D trimmed away -entirely" null has n = 0..42 per locus and is useless): - -| locus | call rate | gene acc | false call | ambiguous | -|---|---|---|---|---| -| human IGH | 69 % | 99 % | 2.0 % | 5 % | -| human TRB | 8 % | 100 % | 0.3 % | 0 % | -| human TRD | 12 % | 100 % | 0.0 % | 0 % | -| mouse TRB | 11 % | 100 % | 0.7 % | 58 % | - -Out of model, on real fixtures: human IGH calls a D on 36 % of records where nt manages 68 %, -agreeing with it on 98 %. SHM and codon degeneracy cost half the recall, almost no precision. - -Mouse TRBD1/TRBD2 translate to near-identical poly-glycine — hence 58 % ambiguity. That is -the honest answer, and it is why `_best_d` reports lists rather than picking a winner. - -`d_germline_*` and `d_cigar` are **withheld** on aa: those offsets index a reading frame, not -the D germline. - -## 5. D-D is limited by surviving length, not orientation - -Order-respecting injection (TRBD1 5′ of TRBD2, the only producible orientation, and it can -only join a J2-cluster J): sensitivity 14.7 % human / 12.9 % mouse TRB, false `d2_call` on a -true single-D junction 0.1 % / 0.0 %. Stratified: **~0 % when either D survives under 6 nt, -57–72 % once both survive 7+.** Trimming usually leaves less. IGH reaches 61 %. - -Orientation was never the limit: only 3 of 223 human calls came back as the wrong pair/order. - -## 6. Dead ends (measured; see ROADMAP) - -* Re-ranking nt D candidates by the length prior (`λS + log P(insVD) + log P(dlen) + - log P(insDJ) + log P(D|J)`) changes **nothing**: IGH gene accuracy 98.9 → 97.8 %, huTRB - 94.2 → 94.5 %, identical call rate and `d_start` error. With 10–18 matched nt, `λS` is - 11–20 nats and the prior moves ±3. In protein the ratio inverts — which is exactly why - `dpost` needs the prior and `_map_d` does not. -* A present/absent Bayes-factor gate buys IGH ~+3 pp recall at matched FP (93.6 % vs 90.7 % - @ ~2 % FP) and nothing for TRD — but needs a **per-locus** threshold (BF > 6 for IGH, - BF > 10 for TRD at the same FP), reintroducing the knobs the E-value removed, and the - priors exist for only 4 of the 15 (organism, D-locus) pairs arda ships. diff --git a/memory/discordance-and-scaling.md b/memory/discordance-and-scaling.md deleted file mode 100644 index 5fc2a52..0000000 --- a/memory/discordance-and-scaling.md +++ /dev/null @@ -1,84 +0,0 @@ -# Discordance analysis & scaling - -## Why arda ever disagrees with IgBLAST (resolved to ~0) -On 100 real IGH mRNA, region concordance is 99.7%. Categorized -(`scripts/analyze_discordance.py`): -- **582 exact**, **18 one-residue boundary** (FR1 start / CDR3 termini — both tools - valid, just trim the partial terminal codon differently), **0 frameshift, 0 other**. - -The earlier ~1-2% "garbage" cases had ONE root cause: when the mmseqs alignment -starts **mid-codon** (e.g. q1↔t3), the projected FR1 start is 1 nt off the true -codon boundary. The first fix (stop-free `detect_coding_frame`) failed on a record -(PX613029.1) where BOTH frames were coincidentally stop-free over the short V. -**Real fix** (transfer.py): derive the query reading frame from the **alignment -phase** — project the first in-phase target codon boundary (`p = tstart + -((t0-tstart)%3)`) to the query; that query position IS `coding_start`. Exact, no -heuristic. detect_coding_frame remains only as a gap fallback. - -## V/J output -Yes — `v_call`/`j_call` come from the best-hit scaffold's allele set (comma-joined -ambiguous calls). V-gene agreement with IgBLAST was 100/100. - -## Prefilter / bulk RNA-seq speed (measured, 16 threads, 150nt) -mmseqs k-mer prefilter rejects non-receptor reads before alignment, so mostly-junk -input is far faster: 100% receptor ~5.7k/s, 10% ~19k/s, 1% ~25k/s -(`scripts/bench_prefilter.py`). This matches blood bulk RNA-seq (~1-5% receptor). - -## CDR3 / junction correctness (AIRR-critical) -CDR3 length is query-specific, so its end is **J-anchored** (= FR4 start − 1), NOT -taken from the fixed-length scaffold — otherwise long somatic CDR3s get truncated. -`cdr3_aa == junction_aa[1:-1]` holds **by construction** for every emitted junction. - -**Out-of-frame junctions are reported** (not dropped). `transfer.py _junction_nt`: -the nt junction is the real query slice (Cys104 codon .. [FW]118 codon). For -translation, when V and J are in different frames (phase = (fwr4_start − -coding_start) % 3 ≠ 0), insert k=(3−phase)%3 N **after the V germline end** -(`v_sequence_end`, clamped to stay inside CDR3 so Cys/[FW] flanks survive); the -codon(s) containing inserted N render as `_`. FR4 is translated in its own J frame -(`translate(Q[fwr4_start-1:])`) so it reads `WGQG…` even for non-productive reads. -`productive` = in-frame AND no stop in V-side/junction. - -**Extended scaffold markup**: `markup.tsv` now stores `v_sequence_end` / -`j_sequence_start` (scaffold nt); these are transferred to queries (point -projection) → AIRR `v_sequence_end`, `j_sequence_start`, locating the V/J split in -the junction. V and J can rarely overlap/cross — handled by clamping. - -**Concordance scoring**: non-canonical junctions (not C…[FW], or containing `_`) -are reported but **excluded** from junction/cdr3 concordance metrics. On the -committed fixtures, productive-canonical junction & cdr3 match IgBLAST ~100% per -species; productive region concordance 98–99.7%. - -## Trimmed inputs & region deletion -- V-only (FR1-FR3) fragments: annotate FR1-3 + CDR1-2 correctly, no CDR3/FR4. -- J-side (CDR3-FR4) fragments: coding frame is derived from the **alignment phase** - (not FR1), so CDR3/FR4 aa are produced even without V. -- Deleting an internal region (e.g. CDR2) collapses it to ~0 residues; flanks and - distal regions still round-trip exactly (nt and aa). Tests in tests/synthetic. - -## Committed test fixtures (offline) -`tests/assets/realworld/.fasta.gz` + `.igblast.airr.tsv.gz`: -balanced ~7.3k GenBank mRNA across all 5 organisms × loci (IG all; TR human/mouse; -NCBI lacks 500 for rare groups like TRG so totals < 10k), gzipped, with IgBLAST AIRR -reference. `tests/realworld` runs offline (mmseqs + DB only), parametrized per -organism. Rebuild via `scripts/build_test_fixtures.py`. - -## Multi-species concordance (productive records) -Region concordance vs IgBLAST: human 98.6%, mouse 99.6%, rat 98.0%, rabbit 99.7%, -rhesus 99.4%. junction/cdr3 vs IgBLAST ~99% all species. KEY: compare only on -IgBLAST-productive records and skip IgBLAST regions with stops — GenBank junk -(genomic/partial/pseudogene, e.g. mouse IGK AM0863xx with stops in FR3) otherwise -drags mouse to ~91%. On productive rearrangements arda ≈ IgBLAST everywhere. - -## Junction emitted only with both flanks -A junction is emitted ONLY when both conserved residues are present (Cys before -cdr3 AND [FW] opening FR4). Truncated queries past the Cys get empty junction -rather than a partial one — so `cdr3_aa == junction_aa[1:-1]` holds for every -emitted junction. (transfer.py `_set_junction`.) - -## 30M-read / 32-core estimate -~25k/s at 16 threads, 1% content → ~35-40k/s at 32 cores (mmseqs scales sublinearly) -→ **~10-20 min for 30M reads**. Same order of magnitude as a STAR genome-mapping -pass (STAR ~10^5 reads/s on many-core nodes, Dobin et al. Bioinformatics 2013 — -approximate, hardware/index dependent; STAR is faster per read but indexes the whole -genome). **Caveat**: current `annotate_file` loads all reads into memory (~10-15 GB -for 30M); needs streaming/sharding (Roadmap TODO) before running at that scale. diff --git a/memory/igblast-gotchas.md b/memory/igblast-gotchas.md deleted file mode 100644 index 8c6c8f0..0000000 --- a/memory/igblast-gotchas.md +++ /dev/null @@ -1,29 +0,0 @@ -# IgBLAST gotchas - -## VJ loci still require -germline_db_D -Even for D-less loci (TRA, TRG, IGK, IGL), `igblastn` errors without -`-germline_db_D` (it falls back to a built-in `_gl_D` not on our search -path). Fix: `airr_extract._dummy_d_db` builds a one-sequence placeholder D -database and passes it for VJ loci. See `refbuild/airr_extract.py`. - -## IGDATA -`igblast.igdata_env()` sets `IGDATA` to `bin/` so IgBLAST finds `internal_data/` -and `optional_file/`. The downloaded release is laid out by -`scripts/fetch_igblast.py` so that `bin/` contains the executables **and** those -two trees. - -## Germline DBs -Build with `makeblastdb -parse_seqids -dbtype nucl` from **ungapped** IMGT files -(via `edit_imgt_file.pl`). Gapped IMGT sequences (with `.`) must not go into -makeblastdb. - -## AIRR output (`-outfmt 19`) -- Coordinates (`*_start`/`*_end`, `v_sequence_start`…) are **1-based, closed**. -- Region nt seqs are `fwr1..cdr3`; AA versions are `fwr1_aa..cdr3_aa` (we reuse - these directly instead of re-translating regions). -- `productive` is `T`/`F`/empty; rev_comp `T`/`F`. - -## Binaries / versions -IgBLAST 1.22.0 (NCBI LATEST as of build). Organisms with internal_data: human, -mouse, rat, rabbit, rhesus_monkey. MMseqs2 = official static binary (no conda), -auto-fetched into `bin/`. diff --git a/memory/junction-markup.md b/memory/junction-markup.md deleted file mode 100644 index 0c61c27..0000000 --- a/memory/junction-markup.md +++ /dev/null @@ -1,111 +0,0 @@ -# Junction markup + repair (`arda.cdr3fix`) - -Marks up a record that has **no read behind it**: a CDR3 amino acid, a V call, a J call, a -species. The VDJdb case. - -## 1. The single biggest correctness trap: junction ≠ CDR3 - -VDJdb's `cdr3` column is the AIRR **junction** — Cys104 → Phe/Trp118, **both anchors -included**. arda's own `cdr3` field *excludes* both, and its `junction` includes them. -`junction_aa` is therefore two residues longer than `cdr3_aa`. - -Everything in `cdr3fix`, `dmap` and `dpost` works in **junction space**. Conflating the two -silently corrupts every coordinate, and downstream corrupts Pgen, clustering and matching. -Assert it in tests; the invariant is `cdr3_aa == junction_aa[1:-1]`. - -## 2. Anchors are derived per allele, and `markup.tsv` is not a valid source - -`markup.tsv` is deduped, comma-joined and scaffold-gated, and already contains a P allele. -VDJdb cites pseudogene and truncated alleles. Derive per allele instead -(`database/vdj//cdr3_anchors.tsv`), keep every functionality (F / ORF / P), and **flag, -never guess** (`status` = `ok` / `truncated` / `no_anchor`). - -Three rules that look right and are wrong: - -* **"IMGT position 104 == gapped nt 310..312"** is false for mouse and rhesus — V-QUEST - carries insertion positions. It silently produced 671 wrong rhesus anchors (first residue - `Y`, not `C`). Use IgBLAST's `internal_data//.ndm.imgt` col 11 (1-based FWR3 - stop) minus 3, with a conserved-FR3-aromatic motif fallback. 1183 agree, 0 disagree. -* **"take the last Cys"** has a 3.8 % error rate: the V CDR3 tail can hold a second Cys - (`YYC-AC-DT` in TRDV2*01, `YYCC…` in IGLV2-11*01). Take the **5′-most** Cys with an - aromatic at i−2 and a tail of 1..14 residues → 0 disagreements. -* **The aux `frame` column is unreliable.** For `TRAJ31*01` it contradicts its own - `cdr3_stop` and yields `QCQTHV`. Frame is `anchor % 3` by construction. Now `NNNARLMF`. - -V-anchor-not-Cys after all this: **0 / 3448**. - -## 2b. Repair always targets a canonical junction - -`cdr3_repaired` is accepted only if it opens with Cys104 and closes with Phe/Trp118 -(`_canonicalise`). A repair exists to restore the anchors, and everything downstream trusts -the result — so `good` implies canonical **by rule**. It used to hold by luck. - -Three constants, three different risks. Do not collapse them: - -* `_MAX_FIX = 2` — residues the repair may **invent** (substitute or restore from germline). - Three of those is a wrong allele, not three typos. -* `_MAX_TRIM = 3` — flanking framework residues it may **remove**. Much smaller risk: those - residues were never explained by the germline. VDJdb trims up to 3 (`GAG`). Bounding trims - by `_MAX_FIX` made the outcome depend on an accident — `YFCASSLGGNEQFF` trimmed but - `YFCASSSRGRGETQYF` did not, because the V's 5th templated residue had been exonuclease- - trimmed and the score fell to −1. -* `_TRIM = -1` — what one trimmed residue costs. **Never zero.** A free trim can *tie* the - untrimmed alignment, and the tie-break prefers consuming more query, so it takes the trim: - clean short IGK junctions came back as `CQQYYSYPF` → `CQQYYSY`, minus their Phe118. Costing - the trim forces it to pay for itself in matches it unlocks. Free-gap OLGA false repairs: - 0.85 %; costed: 0.35 % (the pre-change baseline). - -`v_end` / `j_start` index **`cdr3_repaired`**, not the submission. `v_end` used to be the count -of *consumed input* residues, so any V-side trim or insert shifted it — and `dpost` slices the -non-templated middle with exactly these two coordinates, so a repaired record got its D placed -3 nt off. - -`v_canonical` / `j_canonical` likewise describe the **repaired** junction, which is what -VDJdb's `vCanonical`/`jCanonical` mean. Reading them off the submission disagreed with VDJdb on -76 of 250 fixture rows — every record whose terminal Phe arda had just restored. - -## 3. Detection and repair are different decisions - -A mismatch inside the templated window is either a curation error or simply the V/N boundary -(exonuclease trims the germline, so `templated_aa` is an upper bound, not a promise). - -Repairing everything corrupted 84/3000 VDJdb records: a single mismatch needs only two -flanking matches to outscore stopping (~1/400 by chance), so `CASSPRRY-N-L-QFF` was rewritten -to `…NEQFF` against TRBJ2-1's `SYNEQFF`. Fix: `_MAX_REPLACE = 1` — only anchor-adjacent edits -are **applied**; everything deeper is **reported** with `Cdr3Error.applied = False`. - -Measured on 102,990 VDJdb records: `v_end` 98.66 %, `j_start` 96.52 %, idempotent 99.99 %. -~33k rec/s, pure Python. On the committed 250-row fixture, VDJdb's repair is reproduced on -**100/100** records it flags, with zero novel rewrites. - -## 4. The aligner - -Semi-global Needleman–Wunsch anchored at the conserved residue, free end gaps toward the N -region. `match=1, mismatch=-1, gap=-2`. - -* **Leading germline gaps are free** (`s[i][0] = 0`, `s[0][j] = j*gap`). Without this a - truncated submission is unrepairable — `CASSRGSVRLGTTDPQ` missing its `YF` scored 0. - 94.1 % → 96.4 %. -* **The tie-break must be `max((s[i][j], j, -i))`** — prefer consuming query, then less - germline. The obvious `max((s, i+j, i, j))` made the aligner *prepend* `CA` to - `CYVPGDRGGYTDKLIF`, which already began with the conserved Cys. False repairs 0.63 % → - 0.35 %, under-reports 1 → 0. - -**This bug was invisible to every VDJdb concordance check and only surfaced against OLGA -generative ground truth**, which knows the true `delV`/`delJ`. That is what the synthetic -tier is for. - -## 5. Ground truth harness - -`choose_random_recomb_events()` (not `gen_rnd_prod_CDR3()`, which discards the D and the -trims) returns `{V,D,J,delV,delDl,delDr,delJ,insVD,insDJ}`. Note `delV` can be **negative** — -palindromic P-nt *extend* the templated region. OLGA's germlines are not IMGT's (its -`TRBV3-1*01` anchor is 267 vs IMGT's 270), so filter on `germline_matches_imgt`. - -`$ARDA_VDJREARM` supplies human TRD (with D) and TRG, which OLGA lacks — together, all 7 -human chains plus mouse TRA/TRB. - -Invariants that hold on 3,180 generated junctions: **never under-reports**; boundaries within -1 residue 98.3 %; false repairs 0.35 %; an injected anchor-adjacent typo is detected *and* -repaired 93.7 % of the time; a deep typo is reported 99.7 % of the time and **rewritten 0 -times**. diff --git a/memory/markup-transfer.md b/memory/markup-transfer.md deleted file mode 100644 index f382ed6..0000000 --- a/memory/markup-transfer.md +++ /dev/null @@ -1,79 +0,0 @@ -# Runtime markup transfer (Phase 2) - -## Pipeline -`annotate_records` (annotate/mapper.py): read seqs → `mmseqs easy-search` query vs -`database/vdj//alleles{,.aa}.fasta` → best hit per query (max bits) → -`transfer_hit` projects reference region coords onto the query via the C++ -`_markup.transfer_regions` → AIRR TSV. - -- nt input → search-type 3, `markup.tsv` coords (nt space). -- aa input → search-type 1, `markup.aa.tsv` coords (aa space). Same projection code. - -## C++ `transfer_regions` (src/_markup/markup.cpp) -Single walk over `qaln`/`taln`. Inputs/outputs 1-based closed (AIRR). Verified: -- full identity → coords unchanged. -- insertion in query within a region → region span absorbs the inserted bases; - downstream regions shift by the insertion length (validated: +3nt in CDR3 → FR4 - shifts +3, FR4 stays `WGxG`). -- deletion → deleted ref positions contribute no query base. -- 5'-truncated query → uncovered regions return (-1,-1) → emitted as empty. -Self-hit (query == scaffold) reproduces reference coords (~209/210; rare 1-off at a -boundary, acceptable). - -## AIRR record assembly (annotate/transfer.py) -- Region nt seq = `query[qs-1:qe]`; aa via translating query from FR1 query-start. -- `productive` = no stop in the translated V..J span. -- `junction` = CDR3 ± one codon (nt) / ± one residue (aa). -- `v_call`/`j_call` carry the (possibly comma-joined) ambiguous allele set from - the deduped scaffold. - -## mmseqs binary discovery -`$ARDA_MMSEQS` → `bin/mmseqs` → PATH → auto-fetch a static build into `bin/`. No -conda: dev uses a uv `.venv`; `mmseqs` is the official static binary -(`scripts/fetch_mmseqs.py` / `arda._mmseqs_fetch`). Set `$ARDA_MMSEQS` to override. - -## Done since v1 -- **Reverse-complement** nt queries: `--strand 2`; reverse hits (qstart>qend) are - re-oriented (revcomp query, remap coords), `rev_comp=T`. See [[mmseqs-params]]. -- **Target-DB caching**: `_cached_target_db` createdb once under `data/mmseqs_db/`. -- **C++ seq primitives**: translate/detect_frame/reverse_complement/back_translate - in `src/_markup/markup.cpp`, re-exported by `refbuild/translate.py` (mirpy-compatible). - -## D-segment mapping (done) -After V/J transfer, `transfer._map_d` takes the V..J interior of the junction -(query coords `v_sequence_end+1 .. j_sequence_start-1`) and gapless-local-aligns -every locus D germline against it via the C++ `_markup.d_local_align` (match=+1/ -mismatch=-1; mmseqs is unreliable at ~8-31 nt). Best hit (score≥`_D_MIN_SCORE`=6) -→ `d_call` + `d_sequence_start`/`d_sequence_end` (AIRR, 1-based closed, query -space). For D-D loci (IGH/TRD) a second non-overlapping D above `_D2_MIN_SCORE`=7 -→ `d2_call`/`d2_sequence_*`; `np1`/`np2`/`np3` partition the interior between V, -the D(s), and J. D germlines ship in `database/vdj//d_germlines.fasta` -(`>locus|allele`, VDJ loci only), loaded into `Reference.d_germlines`; VJ loci get -no germlines so D mapping is skipped. `build._collect_d_germlines` writes the file -during refbuild (from `imgt.ungap_gene`). Concordance vs IgBLAST (committed -fixtures): TRB/TRD ~97% gene agreement among co-called; IGH ~46-69% (paralogous D -+ SHM). **Limitation**: a long junction can exceed what mmseqs aligns *through*, -so the projected interior collapses and D mapping silently no-ops (lowers recall, -never a wrong call). Tests: `tests/unit` (d_local_align), `tests/synthetic` -(single-D e2e + `_map_d` double-D logic + option toggle), `tests/realworld` -(per-org concordance), `tests/benchmark` (overhead). - -**Optional**: `map_d` (default True) threads through `annotate_records`/ -`annotate_file`/`adapter.annotate_sequences` and CLI `--map-d/--no-map-d`; gated in -`_annotate_chunk` (`dg = ... if seqtype=='nt' and map_d else None`). D mapping is a -short per-hit C++ local alignment, not a new mmseqs pass, so the cost is tiny: -benchmark (`test_d_mapping_overhead`, VDJ-only worst case, M3, 16 threads) shows -**+1.0% @ 2k, +1.3% @ 10k** wall time (~7-15 µs/seq). Negligible on real bulk -RNA-seq where most reads aren't VDJ hits. - -## TODO / known gaps -- `productive` is heuristic (stop-free V..J span), not full AIRR productivity. - -## Streaming I/O (done) -`annotate_file` streams the input in bounded chunks (`_CHUNK_SIZE=50k`) via a -background reader thread (prefetch queue, maxsize 2) that parses the next chunk -while mmseqs annotates the current (subprocess releases the GIL). Reference + -cached target DB load once and are reused across chunks. Memory is flat for -arbitrarily large FASTQ; output is written incrementally (`airr_header` + -`format_rows`). `annotate_records` (in-memory) shares the same `_annotate_chunk` -core. CLI `--chunk-size`. Multi-node sharding still TODO. diff --git a/memory/mmseqs-params.md b/memory/mmseqs-params.md deleted file mode 100644 index 62ba72c..0000000 --- a/memory/mmseqs-params.md +++ /dev/null @@ -1,43 +0,0 @@ -# MMseqs2 parameters & locus handling - -## One DB, all loci -The runtime target is a **single** mmseqs DB built from `database/vdj// -alleles{,.aa}.fasta`, which contains scaffolds for *all* loci (IGH/IGK/IGL/TRA/ -TRB/TRG/TRD). So one search annotates mixed **bulk RNA-seq** across every locus at -once; the per-query `locus` (AIRR field) comes from the best-hit scaffold. - -## Precompiled (shipped) indexes -The mmseqs target DBs are **committed** under `database/vdj//mmseqs//` -(createdb output + a `VERSION` marker), so annotation runs out of the box (~24 MB -total). `_cached_target_db` (annotate/mapper.py): prefer the committed DB **iff its -`VERSION` == local `mmseqs version`** (DBs are version-sensitive); else build once -into `data/mmseqs_db/_` (private cache, never dirties git). `arda -build-index [--force]` (mapper.build_index) regenerates the shipped DBs for the -local version — a maintainer tool, deliberately NOT in setup.sh so end users don't -dirty the committed blobs. CI uses a different mmseqs build → exercises the fallback. - -## Tuned defaults (annotate/mapper.py) -- nt: `--search-type 3`, `-s 7.0`, `--max-seqs 50`, `--strand 2` (both strands), - `-a`; no coverage filter (partial RNA-seq reads must still map). -- aa: `--search-type 1`, `-s 7.0`, `--max-seqs 50`, `-a`. -- Pipeline is `createdb (cached target) + createdb query + search + convertalis`, - NOT `easy-search`, to reuse the target DB across calls. - -Rationale: queries are germline-similar (85–95% id) and short; `-s 7.0` gives -reliable best-hit recall without the cost of 8.5. `--max-seqs 50` is plenty for a -single best germline. Default `-s 5.7` also worked (~97%); 7.0 buys margin on -divergent/short reads. Measured: 98.7% region concordance vs IgBLAST on 100 real -IGH mRNA; 98.9% on 10k synthetic. - -## convertalis gotcha -`convertalis` needs `--search-type 3` for nucleotide results too, or it errors -"unclear if translated or nucleotide search". Wired via `mmseqs.convertalis(search_type=)`. - -## Reverse strand -`--strand 2` finds reverse hits, reported with **qstart > qend** and qaln/taln -already on the coding strand. mapper detects this, reverse-complements the query, -remaps qstart/qend onto it, and sets `rev_comp=T`. CLI `--strand both|forward`. - -## Speedup vs IgBLAST (synthetic IGH, 16 threads) -~4.4× (10k) → 7.9× (100k); arda ~3.3k seq/s vs igblast ~0.42k seq/s. Throughput -rises with N as fixed DB-build cost amortizes. See `scripts/bench_vs_igblast.py`. diff --git a/memory/reading-frames.md b/memory/reading-frames.md deleted file mode 100644 index f2ed3cf..0000000 --- a/memory/reading-frames.md +++ /dev/null @@ -1,24 +0,0 @@ -# Reading frames: the subtle part of scaffold assembly - -## V frame is NOT always 0 -We assumed IMGT V-REGION starts at codon 1 (frame 0). **False for partial-5′ -alleles** — e.g. human `IGHV5-51*05` (245 nt, "partial in 5′ and in 3′") reads -cleanly only in frame 2. The IMGT header field 8 ("codon start") is **always 1** -and useless here; the real signal is the `partial in 5'` flag (field ~14). - -## Fix: auto-detect coding frame -Rather than parse fragile IMGT fields, `translate.detect_coding_frame` picks the -frame (0/1/2) with the fewest stop codons. For a real germline V exactly one -frame is stop-free and it is unique (verified across alleles). We trim the -leading `frame` nt so every V — and thus the whole scaffold — reads in frame 0. - -## J frame comes from the aux file -`bin/optional_file/_gl.aux` column 2 = 0-based "first coding frame start -position". We pad the V–J junction with N so `(len(V_trimmed) + n_pad + jframe) % -3 == 0`. This put the canonical FR4 `WGQGTLVTVSS` exactly at the J translation — -the sanity check that the frame logic is correct. - -## Validation signal -After the fix, ~18/20 human IGH smoke scaffolds were productive with correct -canonical FR4. A few J alleles still yield non-productive/incomplete markup; those -scaffolds are dropped from the reference and counted in `build.log`. diff --git a/memory/scaffold-enumeration.md b/memory/scaffold-enumeration.md deleted file mode 100644 index 8858e09..0000000 --- a/memory/scaffold-enumeration.md +++ /dev/null @@ -1,28 +0,0 @@ -# Scaffold enumeration: V×J only, dedup, N-spacer - -## Decision -The reference DB enumerates **V×J** scaffolds for every locus — NOT the full -V×D×J the original spec mentioned for IGH/TRB/TRD. - -## Why -The FR/CDR region *coordinates* we transfer at runtime are fully determined by: -- **V**: FR1, CDR1, FR2, CDR2, FR3, and the CDR3 start (conserved Cys104). -- **J**: the CDR3 end and FR4 (conserved `[FW]GXG`). - -The D segment lies **inside** CDR3, whose sequence is somatic/query-specific at -runtime. Enumerating D therefore adds no markup information but multiplies the DB -~50× (human IGH: 5,306 V×J vs 254,688 V×D×J). Confirmed with the user. - -For VDJ loci we still insert a short **frame-neutral N spacer** (`DEFAULT_D_SPACER_NT -= 9`, a multiple of 3) where D would sit, so IgBLAST annotates a plausible CDR3 + -FR4. See `refbuild/combinations.py`. - -## Dedup -Byte-identical assembled scaffolds collapse to one DB entry; `combinations.tsv` -records every (V,J) allele pair mapping to each scaffold. Dedup yield is low -(alleles differ by SNPs) — e.g. IGH 5,306 → ~5,124 — but it's free and lets the -AIRR output report ambiguous allele calls correctly (comma-joined `v_call`). - -## Scale (human, measured) -IGH ~5.1k, TRB ~2.6k, TRA ~8.6k unique scaffolds; ~25–30k total across 7 loci per -species. Watch committed-size if this grows; consider gzip if it becomes heavy.