diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..a9eea33 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,17 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- **DegNorm degradation normalization** (rustar-aligner extension, no STAR + counterpart): `--quantMode GeneCoverage` records per-gene, per-exonic-base + coverage during alignment and writes `GeneCoverage.out.bin`; `--runMode + degNorm` merges several samples' coverage files through a port of DegNorm's + rank-one NMF over-approximation and writes `DegNorm.out/` with + `DegradationIndex.tab`, `AdjustedCounts.tab`, `RawCounts.tab`, + `ScaleFactors.tab`, and `Summary.txt`. Parameters: `--degNormIter`, + `--degNormNmfIter`, `--degNormDownsampleRate`, `--degNormMinimaxCoverage`, + `--degNormSkipBaselineSelection`, `--degNormBins`, + `--degNormMinHighCoverage`, `--degNormSampleId`. Off by default; alignment + output is unchanged. See `DIVERGENCE.md` §4.3. - **CLI and output parity: SAM/SJ/read-input knobs and the STAR limit surface** — 30 further STAR 2.7.11b parameters. (`--outSAMorder` came from #145.) diff --git a/CLAUDE.md b/CLAUDE.md index 149b5a3..2568e41 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -32,7 +32,7 @@ Always run `cargo clippy --all-targets`, `cargo fmt --check`, and `cargo test` b ## Current Status -**396 tests passing, 0 clippy warnings.** SE: 8613/8926 compare_sam.py (96.5%; note: lower due to seeded-RNG tie-break PR diverging from STAR's mt19937), **99.815% faithfulness (tie-adjusted)** (8611/8627 non-tie reads exact), 299 tie-breaking diffs excluded. 1 CIGAR-only disagree (ERR12389696.13573895, insertion placement, seed-level tie). **0 STAR-only / 0 rustar-aligner-only SE reads**. PE: **8390 both-mapped** (STAR: 8390), **0 half-mapped**, 0 MAPQ inflations / 0 deflations, **99.883% PE exact faithfulness (tie-adjusted)** (16284/16306, 475 tie-breaking diffs excluded), **0 proper-pair diffs**, **0 NH diffs**. Phase 17.A: `scoreSeedBest` pre-extension. Phase 17.B: per-mate seeding. Phase 17.C: STAR-faithful SCORE-GATE + mappedFilter. Phase 17.D: combined-span penalty fix + dedup ordering. Phase 17.8: `--quantMode GeneCounts`. Phase E fix (2026-04-21): mate_id-aware diagonal dedup. Phase E2 (2026-04-22): STAR-faithful combined-read seeding. Phase E3 (2026-04-22): combined-threshold half-mapped fallback. Phase E4 (2026-04-22): PE-CHECK2 unconditional. Phase E5 (2026-04-23): split_combined_wt n_mismatch propagation. Phase E6 (2026-04-24): tie-adjusted faithfulness metric in assess_faithfulness.py. Phase F1: --runRNGseed + seeded primary tie-break (PR #5). Phase F2: --outSAMattrRGline (PR #6). Phase F3: --quantMode TranscriptomeSAM (PR #7). Phase F4: SJDB insertion into Genome+SA at genomeGenerate (PR #8). Phase G1 (2026-04-29): junction_shifts fix in split_combined_wt (rDNA cross-copy false-splice filter). Phase G2 (2026-04-29): MAX_RECURSION 10k→100k + sa_pos_to_forward overflow fix (ERR12389696.7118031 NH=3→9). Phase 17.2 (2026-04-29): coordinate-sorted BAM output (`--outSAMtype BAM SortedByCoordinate` → `Aligned.sortedByCoord.out.bam`). Phase 17.4 (2026-04-29): `--outReadsUnmapped Fastx` → `Unmapped.out.mate1` / `Unmapped.out.mate2`; writes unmapped + TooManyLoci reads; PE writes both mates for fully-unmapped and half-mapped pairs. Phase 17.6 (2026-05-01): `--outStd SAM/BAM_Unsorted/BAM_SortedByCoordinate` — routes primary alignment output to stdout via `Box` trait dispatch; `SamStdoutWriter`, `BamStdoutWriter`, `SortedBamStdoutWriter` in sam.rs/bam.rs; verified with samtools pipe (967 records). Phase G3 (2026-05-01): SA tie-breaking fix — `compare_suffixes` tie-breaker changed from `pos_b.cmp(&pos_a)` to `packed_a.cmp(&packed_b)` (ascending by packed SA value with strand bit); rustar-aligner SA is now **byte-for-byte identical** to STAR's SA for the yeast genome (10,862 → 0 entry diffs). diff AS: 6→4 cases (4 remaining are rustar-aligner improvements: .844151 VIII 0mm vs STAR VII 6mm, .4972950 spliced vs unspliced mate2). Phase 17.3 (2026-05-01): PE chimeric detection — `detect_inter_mate_chimeric` in `chimeric/detect.rs`; intra-mate multi-cluster chimeric via cluster splitting + mate2 read_pos adjustment; inter-mate chimeric for discordant pairs (diff chr, same strand, or >1Mb); `align_paired_read` returns 4-tuple including `Vec`; no benchmark regression (8390 both-mapped, 0 half-mapped). Phase 17.11 (2026-05-01): `--chimOutType WithinBAM` — chimeric alignments written as supplementary records (FLAG 0x800) in primary BAM; donor record has full SEQ + SA tag; acceptor has FLAG 0x800 + SA tag + empty SEQ; `build_within_bam_records` in `chimeric/output.rs`; `chim_out_junctions()` / `chim_out_within_bam()` helpers in params.rs; supports mixed `--chimOutType Junctions WithinBAM`. Phase 17.7 (2026-05-01): GTF tag parameters — `--sjdbGTFchrPrefix`, `--sjdbGTFfeatureExon`, `--sjdbGTFtagExonParentTranscript`, `--sjdbGTFtagExonParentGene`; `_configured` variants in `junction/gtf.rs`, `quant/mod.rs`, `quant/transcriptome.rs`, `junction/mod.rs`; all 4 production paths thread params; backward-compat wrappers preserve zero test disruption. Phase 17.9 (2026-05-01): `--outBAMcompression` (BGZF level -1–9, default 1; -1/0=NONE, 1-8=flate2 levels, ≥9=BEST) + `--limitBAMsortRAM` (bytes, 0=unlimited; aborts sort if ~400 bytes/record estimate exceeds limit); `bgzf_compression()` + `make_bgzf_writer()` helpers in `io/bam.rs`; threaded through all 4 BAM writers (unsorted file, sorted file, unsorted stdout, sorted stdout). PE chimericDetectionOld (2026-05-01): per-mate `detect_chimeric_old` called on `all_m1_transcripts` / `all_m2_transcripts` pools after `filter_paired_transcripts` in `read_align.rs`. Phase 17.12 (2026-05-01): BySJout disk buffering — `BySJReadMeta` struct + `NamedTempFile` SAM temp file replaces `Vec`; `create_bysj_writer` / `bysj_write_records` / `bysj_read_n_records` helpers in `io/sam.rs`; `tempfile` moved to `[dependencies]`. Phase 17.13 (2026-05-01): 8 integration tests in `tests/alignment_features.rs` — synthetic 20kb genome with planted GT-AG intron; tests cover BAM output, PE alignment, spliced reads, BySJout, GeneCounts, unmapped output, two-pass mode. Phase 12.2 (2026-05-04): SE chimeric Tier 1b soft-clip re-mapping — `detect_from_soft_clips` in `chimeric/detect.rs` re-seeds the primary alignment's soft-clipped bases when `detect_chimeric_old` finds no partner; `adjust_read_positions` helper shifts sub-seq coords into full-read space for right clips; called as Step 3c in `read_align.rs`. Phase 17.10 (2026-05-04): Chimeric Tier 3 — `detect_from_chimeric_residuals` in `chimeric/detect.rs` re-seeds outer uncovered read regions (before donor / after acceptor) of each found chimeric pair; enables 3-way gene-fusion detection; called as Step 3d in `read_align.rs`. See [ROADMAP.md](ROADMAP.md) for detailed phase tracking and [docs-old/](docs-old/) for per-phase development notes. The published Astro Starlight docs site is in [docs/](docs/). +**396 tests passing, 0 clippy warnings.** SE: 8613/8926 compare_sam.py (96.5%; note: lower due to seeded-RNG tie-break PR diverging from STAR's mt19937), **99.815% faithfulness (tie-adjusted)** (8611/8627 non-tie reads exact), 299 tie-breaking diffs excluded. 1 CIGAR-only disagree (ERR12389696.13573895, insertion placement, seed-level tie). **0 STAR-only / 0 rustar-aligner-only SE reads**. PE: **8390 both-mapped** (STAR: 8390), **0 half-mapped**, 0 MAPQ inflations / 0 deflations, **99.883% PE exact faithfulness (tie-adjusted)** (16284/16306, 475 tie-breaking diffs excluded), **0 proper-pair diffs**, **0 NH diffs**. Phase 17.A: `scoreSeedBest` pre-extension. Phase 17.B: per-mate seeding. Phase 17.C: STAR-faithful SCORE-GATE + mappedFilter. Phase 17.D: combined-span penalty fix + dedup ordering. Phase 17.8: `--quantMode GeneCounts`. Phase E fix (2026-04-21): mate_id-aware diagonal dedup. Phase E2 (2026-04-22): STAR-faithful combined-read seeding. Phase E3 (2026-04-22): combined-threshold half-mapped fallback. Phase E4 (2026-04-22): PE-CHECK2 unconditional. Phase E5 (2026-04-23): split_combined_wt n_mismatch propagation. Phase E6 (2026-04-24): tie-adjusted faithfulness metric in assess_faithfulness.py. Phase F1: --runRNGseed + seeded primary tie-break (PR #5). Phase F2: --outSAMattrRGline (PR #6). Phase F3: --quantMode TranscriptomeSAM (PR #7). Phase F4: SJDB insertion into Genome+SA at genomeGenerate (PR #8). Phase G1 (2026-04-29): junction_shifts fix in split_combined_wt (rDNA cross-copy false-splice filter). Phase G2 (2026-04-29): MAX_RECURSION 10k→100k + sa_pos_to_forward overflow fix (ERR12389696.7118031 NH=3→9). Phase 17.2 (2026-04-29): coordinate-sorted BAM output (`--outSAMtype BAM SortedByCoordinate` → `Aligned.sortedByCoord.out.bam`). Phase 17.4 (2026-04-29): `--outReadsUnmapped Fastx` → `Unmapped.out.mate1` / `Unmapped.out.mate2`; writes unmapped + TooManyLoci reads; PE writes both mates for fully-unmapped and half-mapped pairs. Phase 17.6 (2026-05-01): `--outStd SAM/BAM_Unsorted/BAM_SortedByCoordinate` — routes primary alignment output to stdout via `Box` trait dispatch; `SamStdoutWriter`, `BamStdoutWriter`, `SortedBamStdoutWriter` in sam.rs/bam.rs; verified with samtools pipe (967 records). Phase G3 (2026-05-01): SA tie-breaking fix — `compare_suffixes` tie-breaker changed from `pos_b.cmp(&pos_a)` to `packed_a.cmp(&packed_b)` (ascending by packed SA value with strand bit); rustar-aligner SA is now **byte-for-byte identical** to STAR's SA for the yeast genome (10,862 → 0 entry diffs). diff AS: 6→4 cases (4 remaining are rustar-aligner improvements: .844151 VIII 0mm vs STAR VII 6mm, .4972950 spliced vs unspliced mate2). Phase 17.3 (2026-05-01): PE chimeric detection — `detect_inter_mate_chimeric` in `chimeric/detect.rs`; intra-mate multi-cluster chimeric via cluster splitting + mate2 read_pos adjustment; inter-mate chimeric for discordant pairs (diff chr, same strand, or >1Mb); `align_paired_read` returns 4-tuple including `Vec`; no benchmark regression (8390 both-mapped, 0 half-mapped). Phase 17.11 (2026-05-01): `--chimOutType WithinBAM` — chimeric alignments written as supplementary records (FLAG 0x800) in primary BAM; donor record has full SEQ + SA tag; acceptor has FLAG 0x800 + SA tag + empty SEQ; `build_within_bam_records` in `chimeric/output.rs`; `chim_out_junctions()` / `chim_out_within_bam()` helpers in params.rs; supports mixed `--chimOutType Junctions WithinBAM`. Phase 17.7 (2026-05-01): GTF tag parameters — `--sjdbGTFchrPrefix`, `--sjdbGTFfeatureExon`, `--sjdbGTFtagExonParentTranscript`, `--sjdbGTFtagExonParentGene`; `_configured` variants in `junction/gtf.rs`, `quant/mod.rs`, `quant/transcriptome.rs`, `junction/mod.rs`; all 4 production paths thread params; backward-compat wrappers preserve zero test disruption. Phase 17.9 (2026-05-01): `--outBAMcompression` (BGZF level -1–9, default 1; -1/0=NONE, 1-8=flate2 levels, ≥9=BEST) + `--limitBAMsortRAM` (bytes, 0=unlimited; aborts sort if ~400 bytes/record estimate exceeds limit); `bgzf_compression()` + `make_bgzf_writer()` helpers in `io/bam.rs`; threaded through all 4 BAM writers (unsorted file, sorted file, unsorted stdout, sorted stdout). PE chimericDetectionOld (2026-05-01): per-mate `detect_chimeric_old` called on `all_m1_transcripts` / `all_m2_transcripts` pools after `filter_paired_transcripts` in `read_align.rs`. Phase 17.12 (2026-05-01): BySJout disk buffering — `BySJReadMeta` struct + `NamedTempFile` SAM temp file replaces `Vec`; `create_bysj_writer` / `bysj_write_records` / `bysj_read_n_records` helpers in `io/sam.rs`; `tempfile` moved to `[dependencies]`. Phase 17.13 (2026-05-01): 8 integration tests in `tests/alignment_features.rs` — synthetic 20kb genome with planted GT-AG intron; tests cover BAM output, PE alignment, spliced reads, BySJout, GeneCounts, unmapped output, two-pass mode. Phase 12.2 (2026-05-04): SE chimeric Tier 1b soft-clip re-mapping — `detect_from_soft_clips` in `chimeric/detect.rs` re-seeds the primary alignment's soft-clipped bases when `detect_chimeric_old` finds no partner; `adjust_read_positions` helper shifts sub-seq coords into full-read space for right clips; called as Step 3c in `read_align.rs`. Phase 17.10 (2026-05-04): Chimeric Tier 3 — `detect_from_chimeric_residuals` in `chimeric/detect.rs` re-seeds outer uncovered read regions (before donor / after acceptor) of each found chimeric pair; enables 3-way gene-fusion detection; called as Step 3d in `read_align.rs`. Phase 18.1 (2026-08-20): DegNorm degradation normalization (not a STAR feature) — `--quantMode GeneCoverage` accumulates per-gene, per-exonic-base coverage on the `GeneCounts` unique-hit path and writes `GeneCoverage.out.bin` (gzip; header, per-gene table, gene ids, flat u32 coverage); `--runMode degNorm` merges N sample files through a port of DegNorm's rank-one NMF over-approximation (`degnorm/nmf.rs`), baseline selection (`degnorm/baseline.rs`) and outer depth-factor loop (`degnorm/run.rs`), writing `DegNorm.out/{DegradationIndex,AdjustedCounts,RawCounts,ScaleFactors}.tab` + `Summary.txt`. DI needs >= 2 samples. See [DIVERGENCE.md](DIVERGENCE.md) §4.3. See [ROADMAP.md](ROADMAP.md) for detailed phase tracking and [docs-old/](docs-old/) for per-phase development notes. The published Astro Starlight docs site is in [docs/](docs/). ## Source Layout @@ -72,6 +72,12 @@ src/ gtf.rs -- GTF parser (internal) quant/ mod.rs -- Gene-level read counting (--quantMode GeneCounts, ReadsPerGene.out.tab) + coverage.rs -- Per-gene exonic coverage (--quantMode GeneCoverage, GeneCoverage.out.bin) + degnorm/ + mod.rs -- --runMode degNorm entry point (DegNorm degradation normalization, not STAR) + nmf.rs -- Rank-one NMF over-approximation (port of DegNorm nmf.py) + baseline.rs -- Baseline selection + per-gene DI scores + run.rs -- Multi-sample driver + DegNorm.out/ tables chimeric/ mod.rs -- Module exports detect.rs -- Chimeric detection (Tier 1: transcript-pair, Tier 1b: soft-clip re-seed, Tier 2: multi-cluster, Tier 3: residual re-seed) diff --git a/DIVERGENCE.md b/DIVERGENCE.md index 66f4088..de4ef61 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -127,6 +127,20 @@ For `--quantMode TranscriptomeSAM`, rustar-aligner builds the per-transcript exo rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` crate, avoiding the `getrandom`/`zerocopy`/`ppv-lite86` dependency chain. This is the generator underlying §1.1; it is called out separately because it is a dependency/implementation choice independent of the tie-break policy. It is not the only in-tree generator: `--soloCellFilter EmptyDrops_CR` samples with a bit-exact libc++ `mt19937` (`src/solo/libcxx_rng.rs`) so its Monte-Carlo null matches STAR's — a convergence with STAR rather than a divergence from it. +### 4.3 DegNorm degradation normalization (`--quantMode GeneCoverage`, `--runMode degNorm`) + +**What STAR does.** Nothing: STAR has no transcript-degradation model and no per-gene coverage output. + +**What rustar-aligner does.** Two optional, off-by-default additions. `--quantMode GeneCoverage` accumulates per-gene, per-exonic-base coverage during alignment and writes `GeneCoverage.out.bin`. `--runMode degNorm` merges several samples' coverage files and fits the rank-one NMF over-approximation of [DegNorm](https://nustatbioinfo.github.io/DegNorm/) (Xiong et al., *Genome Biology* 2019), writing Degradation Index scores and degradation-adjusted counts to `DegNorm.out/`. + +**Why.** Degradation bias is gene- and sample-specific, so a global size factor cannot remove it. DegNorm normally recovers coverage curves by re-reading sorted, indexed BAMs; the aligner already assigns each unique read to a gene for `ReadsPerGene.out.tab`, so the curves cost almost nothing to produce there. + +**Impact.** None on alignment output: coverage comes from the same uniquely mapped, unambiguously assigned reads that feed `ReadsPerGene.out.tab` column 1, and the integration suite asserts the SAM records are identical with and without the flag. A degradation index is defined across samples, so `--runMode degNorm` requires at least two coverage files and cannot run inside a single alignment. + +**Source.** `src/quant/coverage.rs`, `src/degnorm/` (ported from DegNorm's `degnorm/nmf.py`). No STAR source: this feature has no STAR counterpart. + +--- + ## 5. Known residual single-read differences These are **not** deliberate divergences — they are tracked residual diffs on the 10k yeast benchmark, kept here for completeness. Each is a single read; none is a systematic behaviour difference. diff --git a/README.md b/README.md index d6536e6..efc93c6 100644 --- a/README.md +++ b/README.md @@ -216,6 +216,7 @@ resident; the 16 GB sparse index is stable at ~54 s. - Chimeric alignment detection — SE and PE, 4-tier pipeline: transcript-pair search, multi-cluster, soft-clip re-seeding, residual outer re-seeding for multi-junction fusions (`--chimSegmentMin`) - Gene-level read counting (`--quantMode GeneCounts` → `ReadsPerGene.out.tab`) - Transcriptome-coordinate SAM output (`--quantMode TranscriptomeSAM`) +- Transcript-degradation normalization, DegNorm-style (`--quantMode GeneCoverage` during alignment, then `--runMode degNorm` across samples) — per-gene Degradation Index scores and degradation-adjusted counts; not a STAR feature - **Single-cell quantification (STARsolo)** — `--soloType CB_UMI_Simple`, `CB_UMI_Complex` (multi-segment barcodes), and `SmartSeq` (plate-based, SE + PE); features `Gene`, `GeneFull` (pre-mRNA), `SJ`, and `Velocyto` (spliced/unspliced/ambiguous); barcode correction (`--soloCBmatchWLtype` Exact/1MM/1MM_multi/…), UMI dedup (`--soloUMIdedup` 1MM_All/1MM_CR/1MM_Directional/…), multi-gene UMI filtering, multi-mapper resolution (`--soloMultiMappers` Uniform/PropUnique/EM/Rescue), cell calling (`--soloCellFilter` CellRanger2.2/TopCells/EmptyDrops_CR), gzip output, and `Summary.csv` — writes STARsolo-compatible `Solo.out//{raw,filtered}/{matrix.mtx, barcodes.tsv, features.tsv}` - WASP allele-specific-mapping filter (`--waspOutputMode SAMtag`, `--varVCFfile`) — vW/vA/vG tags, SE + PE - Paired-end mate-overlap merging (`--peOverlapNbasesMin`, `--peOverlapMMp`) diff --git a/ROADMAP.md b/ROADMAP.md index d419f61..a993871 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -348,3 +348,22 @@ Single-cell quantification layered around the existing aligner: the cDNA read al **Phase 14.5–14.11 + performance** (2026-07): completed the feature-parity set — `Summary.csv` (STARsolo-faithful, CellRanger funnel split to its own file), `--soloCellFilter` CellRanger2.2/TopCells/**EmptyDrops_CR** (Monte-Carlo ambient rescue in the `filtered/` writer), `--soloFeatures` **GeneFull/SJ/Velocyto** (spliced/unspliced/ambiguous per Sullivan 2025), `--soloMultiMappers` Uniform/PropUnique/EM/Rescue, chemistries **CB_UMI_Complex** (multi-segment) and **SmartSeq** (plate-based, SE + PE fragment counts), and a rustar-vs-STARsolo SJ + multi-mapper diff harness. Performance: pipelined solo FASTQ decode, parallelized matrix build + EmptyDrops MC, libdeflate/zlib-rs for matrix gzip + BGZF, and an **O(log n + k) segment-tree gene-overlap query** (replacing STAR's linear scan — the #1 solo hotspot, ~14% wall reduction). Sparse suffix array (`--genomeSAsparseD`, byte-identical to STAR's D=2) for a 31% smaller index. 516 tests, 0 clippy warnings. **Native three-way benchmark** (2026-07, `test/aws/`): fresh single-instance EC2 comparison on a real 10x dataset (`5k_Mouse_PBMCs_5p_gem-x_GEX`, 5′ GEM-X, GRCm39-2024-A), all native x86_64, 10 threads, NVMe, page cache dropped, no BAM. Wall / peak RSS / cells: **STARsolo 2.7.11b** 87 s / 28.3 GB / 4,061; **rustar-aligner** 121 s / 25.7 GB / 3,689 (→ ~105 s with the segment-tree query merged after this run); **rustar `--genomeSAsparseD 2`** 119 s / **17.7 GB** / 3,692; **CellRanger 10.0.0** 347 s / 13.1 GB / 3,858. `Gene/raw` matrix byte-identical to STARsolo's. Supersedes the earlier Docker-emulation numbers above (those were penalized by Rosetta/virtiofs). Remaining gap to STARsolo is small and output-identical; rustar owns the memory frontier via sparse SA. + +--- + +## Phase 18: DegNorm degradation normalization (not a STAR feature) + +| Sub-phase | Item | Status | +|-----------|------|--------| +| 18.1 | `--quantMode GeneCoverage` — per-gene, per-exonic-base coverage during alignment | ✅ Complete | +| 18.2 | `--runMode degNorm` — multi-sample NMF-OA fit, DI scores, adjusted counts | ✅ Complete | +| 18.3 | Validation against the Python DegNorm on a real multi-sample dataset | ⬜ Planned | +| 18.4 | Coverage-curve plots / estimated coverage matrices | ⬜ Not planned | + +**Phase 18.1–18.2** (2026-08-20): DegNorm (Xiong et al., *Genome Biology* 2019) corrects RNA-seq counts for sample- and gene-specific transcript degradation. Since the model is a rank-one over-approximation *across samples*, a Degradation Index cannot be produced inside a single alignment run; the work splits in two. + +Phase 1, during alignment: `--quantMode GeneCoverage` allocates a flat `AtomicU32` array over the concatenated merged exons of every gene (transcript space) and accumulates coverage on the same unique-hit path `GeneCounts` uses (unique alignment, exactly one overlapping gene). Paired-end mate blocks are merged first, so an overlapping pair contributes one per base. Output `GeneCoverage.out.bin`: a gzip stream with magic `RSDGNCOV`, a per-gene table (length + raw count), the gene id block, and the flat coverage array. `src/quant/coverage.rs`. + +Phase 2, `--runMode degNorm`: loads N coverage files, validates a shared gene model, and runs the port in `src/degnorm/` — `nmf.rs` (leading singular triplet by power iteration on the `p x p` Gram matrix, dual-ascent NMF over-approximation, `ratio_svd`), `baseline.rs` (high-coverage filter, bin-dropping baseline search, DI with DegNorm's `+1` denominator and `[0, 0.9]` clamp), `run.rs` (initial depth factors from low-DI genes, `--degNormIter` outer rounds, sample-average DI for genes that skip baseline selection, `rayon` fan-out over genes). Writes `DegNorm.out/{DegradationIndex,AdjustedCounts,RawCounts,ScaleFactors}.tab` and `Summary.txt`. + +Alignment output is unchanged when `GeneCoverage` is on (asserted in `tests/degnorm.rs`). Recorded in [DIVERGENCE.md](DIVERGENCE.md) §4.3; user guide at `docs/src/content/docs/guides/degnorm.md`; design spec at `docs/superpowers/specs/2026-08-20-degnorm-design.md`. diff --git a/docs/astro.config.mjs b/docs/astro.config.mjs index 548198d..cb38882 100644 --- a/docs/astro.config.mjs +++ b/docs/astro.config.mjs @@ -55,6 +55,7 @@ export default defineConfig({ { label: 'Two-pass mode', slug: 'guides/two-pass' }, { label: 'Chimeric detection', slug: 'guides/chimeric' }, { label: 'Gene quantification', slug: 'guides/quantification' }, + { label: 'Degradation normalization (DegNorm)', slug: 'guides/degnorm' }, { label: 'Migrating from STAR', slug: 'guides/migrating-from-star' }, ], }, diff --git a/docs/src/content/docs/guides/degnorm.md b/docs/src/content/docs/guides/degnorm.md new file mode 100644 index 0000000..5198231 --- /dev/null +++ b/docs/src/content/docs/guides/degnorm.md @@ -0,0 +1,136 @@ +--- +title: Degradation normalization (DegNorm) +description: Capture per-gene coverage during alignment, then correct read counts for sample- and gene-specific transcript degradation. +--- + +RNA degrades unevenly. Long transcripts lose more material than short ones, and +the loss differs from sample to sample, so a single global size factor cannot +undo it. [DegNorm](https://nustatbioinfo.github.io/DegNorm/) (Xiong et al., +*Genome Biology* 2019) models this directly: it fits each gene's coverage matrix +across samples with a rank-one over-approximation, splitting a shared coverage +envelope from per-sample abundances, and reports a Degradation Index (DI) score +per gene per sample plus degradation-adjusted counts. + +rustar-aligner implements this in two phases: + +1. **During alignment**, `--quantMode GeneCoverage` records per-gene, per-exonic + base coverage and writes `GeneCoverage.out.bin`. +2. **After alignment**, `--runMode degNorm` merges several samples' coverage + files and fits the model. + +:::note[This is not a STAR feature] +Both flags are rustar-aligner extensions with no STAR counterpart, off by +default, and documented in `DIVERGENCE.md`. Alignment output is unchanged when +they are enabled. +::: + +## Why two phases + +The DI score is defined *across* samples. With one sample the rank-one fit is +exact and every DI is zero, so the model needs at least two coverage files and +cannot run inside a single alignment. What a single run can do, essentially for +free, is the expensive part DegNorm normally pays for by writing a sorted BAM, +indexing it, and re-reading it: computing the coverage curves themselves. + +## Phase 1: capture coverage + +```bash +for sample in ctrl1 ctrl2 treat1 treat2; do + rustar-aligner \ + --genomeDir /path/to/index \ + --readFilesIn ${sample}.fastq.gz \ + --readFilesCommand zcat \ + --sjdbGTFfile annotations.gtf \ + --quantMode GeneCounts GeneCoverage \ + --outFileNamePrefix out/${sample}_ +done +``` + +This writes `out/${sample}_GeneCoverage.out.bin` alongside the usual output. + +Coverage is accumulated from the same reads that feed `ReadsPerGene.out.tab` +column 1: uniquely mapped, and assigned to exactly one gene. Multimappers and +gene-ambiguous reads are skipped, matching DegNorm's default. For paired-end +data the fragment is one observation: mate blocks are merged first, so an +overlapping pair contributes one to each base it covers, not two. + +| Flag | Default | Meaning | +|---|---|---| +| `--quantMode GeneCoverage` | off | Enable coverage capture (needs `--sjdbGTFfile`) | +| `--degNormSampleId` | basename of `--outFileNamePrefix` | Sample name stored in the file and used as a column header later | + +Memory: about 4 bytes per exonic base of the annotation, roughly 280 MB for +human GENCODE. Yeast or a chromosome-scale index costs a few megabytes. + +## Phase 2: fit the model + +```bash +rustar-aligner \ + --runMode degNorm \ + --degNormCoverageFiles out/ctrl1_GeneCoverage.out.bin \ + out/ctrl2_GeneCoverage.out.bin \ + out/treat1_GeneCoverage.out.bin \ + out/treat2_GeneCoverage.out.bin \ + --outFileNamePrefix out/ +``` + +This phase touches no genome index and no read file: it is pure CPU work over +the coverage matrices, parallelised over genes. All samples must have been +aligned against the same GTF; a mismatching gene model is a fatal error naming +the first gene that differs. + +| Flag | Default | DegNorm CLI equivalent | +|---|---|---| +| `--degNormCoverageFiles` | (required, at least 2) | `--bam-files` | +| `--degNormIter` | 5 | `--iter` | +| `--degNormNmfIter` | 100 | `--nmf-iter` | +| `--degNormDownsampleRate` | 1 (no downsampling) | `--downsample-rate` | +| `--degNormMinimaxCoverage` | 0 | `--minimax-coverage` | +| `--degNormSkipBaselineSelection` | off | `--skip-baseline-selection` | +| `--degNormBins` | 20 | (baseline-selection bins) | +| `--degNormMinHighCoverage` | 50 | (minimum high-coverage positions) | + +`--runRNGseed` sets the offset used by systematic downsampling, so a run is +reproducible regardless of thread count. + +## Output + +Everything lands in `DegNorm.out/`: + +| File | Contents | +|---|---| +| `DegradationIndex.tab` | DI per gene (rows) per sample (columns). 0 = no detected degradation, capped at 0.9 | +| `AdjustedCounts.tab` | Depth-normalised counts divided by `1 - DI` | +| `RawCounts.tab` | The unique counts that entered the model, for provenance | +| `ScaleFactors.tab` | Final per-sample sequencing-depth factors | +| `Summary.txt` | Parameters, gene counts, genes sent through baseline selection per iteration, median DI per sample | + +`AdjustedCounts.tab` is the matrix to hand to a downstream differential +expression tool. + +## How the fit works + +For each gene, the coverage matrix `F` (samples x exonic positions) is scaled by +the current depth factors, then fitted by a rank-one over-approximation: dual +ascent lifts the fit above the observed coverage, so the envelope represents the +undegraded curve rather than the average one. The DI is the area between the +scaled envelope and the observed coverage, divided by the area under the +envelope. + +Baseline selection refines this. When a gene shows degradation, the transcript +is split into bins and the bin with the worst relative residual is dropped +repeatedly, until what remains is a region where every sample agrees. The +envelope re-estimated from that baseline region gives a cleaner DI for the whole +transcript. Genes that never enter baseline selection inherit the sample-average +DI, as upstream does. + +The outer loop then folds the degradation correction back into the depth +factors, and repeats `--degNormIter` times. + +## Limitations + +- At least two samples are required. +- Coverage must come from phase 1; third-party BAMs are not read. +- Coverage-curve plots and the estimated coverage matrices that upstream stores + as `.pkl` files are not produced. +- No MPI or warm-start directories. diff --git a/docs/superpowers/plans/2026-08-20-degnorm.md b/docs/superpowers/plans/2026-08-20-degnorm.md new file mode 100644 index 0000000..8302312 --- /dev/null +++ b/docs/superpowers/plans/2026-08-20-degnorm.md @@ -0,0 +1,2159 @@ +# DegNorm Implementation Plan + +> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. + +**Goal:** Add DegNorm-style transcript-degradation normalization to rustar-aligner: per-gene exonic coverage captured during alignment, then a merge run that computes Degradation Index (DI) scores and degradation-adjusted read counts across samples. + +**Architecture:** Phase A adds `--quantMode GeneCoverage`, an accumulator hanging off the existing `QuantContext` that records per-gene, per-exonic-base coverage on the same unique-hit path `GeneCounts` uses, and writes `GeneCoverage.out.bin`. Phase B adds `--runMode degNorm`, a pure-CPU pass that loads N such files and runs a faithful port of DegNorm's rank-one NMF over-approximation with baseline selection. + +**Tech Stack:** Rust 2024, clap derive params, `flate2` (already a dependency) for the coverage file, `rayon` for per-gene parallelism, `std::sync::atomic` for lock-free accumulation. + +**Spec:** `docs/superpowers/specs/2026-08-20-degnorm-design.md` + +## Global Constraints + +- Crate is `rustar-aligner`, lib name `rustar_aligner`; STAR-style camelCase identifiers are allowed (`#![allow(non_snake_case)]`). +- New CLI flags use clap `#[arg(long = "camelCase")]`; multi-value flags need explicit `num_args`. +- No new dependencies. `flate2`, `rayon`, `tempfile` (dev) are already in `Cargo.toml`. +- `cargo clippy --all-targets` must stay at **0 warnings**; `cargo fmt --check` must pass; all tests must pass before a task is complete. +- Alignment output must remain byte-identical whether or not `GeneCoverage` is enabled. +- Everything here diverges from STAR and must be recorded in `DIVERGENCE.md`; never present it as STAR-faithful, and never invent a STAR flag name. +- Numeric behaviour is ported from DegNorm `degnorm/nmf.py` (v0.1.4, GitHub NUStatBioinfo/DegNorm) including its quirks (`+1` denominators, `[0, 0.9]` clamping). Do not "fix" them. +- Do not commit anything the user did not ask for; commit messages must not use em-dashes. + +--- + +### Task 1: Transcript-space coordinate mapping + +**Files:** +- Create: `src/quant/coverage.rs` +- Modify: `src/quant/mod.rs` (add `pub mod coverage;` near the top, after the existing `use` block) + +**Interfaces:** +- Consumes: `GeneAnnotation.gene_exons: Vec>` (merged, sorted, absolute coords), `GeneAnnotation::n_genes()`. +- Produces: `GeneCoverage::new(&GeneAnnotation) -> GeneCoverage`, `GeneCoverage::gene_len(&self, g: usize) -> u32`, `GeneCoverage::add_block(&self, g: usize, start: u64, end: u64)`, `GeneCoverage::gene_slice(&self, g: usize) -> Vec`. + +- [ ] **Step 1: Write the failing test** + +Append to `src/quant/coverage.rs`: + +```rust +#[cfg(test)] +mod tests { + use super::*; + use crate::quant::GeneAnnotation; + + /// Two-exon gene: [100,110) and [200,205). Transcript length 15. + fn two_exon_ann() -> GeneAnnotation { + GeneAnnotation { + gene_ids: vec!["G1".into()], + gene_names: vec!["G1".into()], + gene_is_reverse: vec![false], + gene_exons: vec![vec![(100, 110), (200, 205)]], + ..Default::default() + } + } + + #[test] + fn block_inside_first_exon_maps_to_transcript_prefix() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + assert_eq!(cov.gene_len(0), 15); + cov.add_block(0, 102, 105); + assert_eq!(cov.gene_slice(0), vec![0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]); + } + + #[test] + fn block_spanning_intron_covers_both_exons_only() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + // A block covering 105..202 in genome space: 5 bases of exon 1, 2 of exon 2. + cov.add_block(0, 105, 202); + assert_eq!(cov.gene_slice(0), vec![0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]); + } + + #[test] + fn block_outside_all_exons_is_ignored() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.add_block(0, 120, 150); + assert_eq!(cov.gene_slice(0).iter().sum::(), 0); + } +} +``` + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --lib quant::coverage` +Expected: FAIL, module `coverage` not found / `GeneCoverage` undefined. + +- [ ] **Step 3: Write minimal implementation** + +Create `src/quant/coverage.rs` above the test module: + +```rust +//! Per-gene, per-exonic-base coverage accumulation for `--quantMode GeneCoverage`. +//! +//! Not a STAR feature. This is the input half of the DegNorm degradation +//! normalization pipeline (see `src/degnorm/`), captured during alignment so +//! that no sorted BAM has to be written, indexed, and re-read. + +use std::sync::atomic::{AtomicU32, AtomicU64, Ordering}; + +use crate::quant::GeneAnnotation; + +/// Lock-free per-gene coverage in transcript space (merged exons, concatenated). +pub struct GeneCoverage { + /// Prefix sums of per-gene transcript lengths; len = n_genes + 1. + offsets: Vec, + /// Flat coverage array, len = offsets[n_genes]. + cov: Vec, + /// Per-gene raw unique read/fragment counts (GeneCounts column 1 rule). + counts: Vec, + /// Total reads/fragments counted into any gene (library size). + n_counted: AtomicU64, +} + +impl GeneCoverage { + pub fn new(ann: &GeneAnnotation) -> Self { + let n = ann.n_genes(); + let mut offsets = Vec::with_capacity(n + 1); + let mut acc: u64 = 0; + offsets.push(0); + for g in 0..n { + let len: u64 = ann.gene_exons[g].iter().map(|&(s, e)| e - s).sum(); + acc += len; + offsets.push(acc); + } + GeneCoverage { + offsets, + cov: (0..acc).map(|_| AtomicU32::new(0)).collect(), + counts: (0..n).map(|_| AtomicU32::new(0)).collect(), + n_counted: AtomicU64::new(0), + } + } + + pub fn n_genes(&self) -> usize { + self.counts.len() + } + + pub fn gene_len(&self, g: usize) -> u32 { + (self.offsets[g + 1] - self.offsets[g]) as u32 + } + + pub fn total_len(&self) -> u64 { + self.offsets[self.offsets.len() - 1] + } + + pub fn n_counted(&self) -> u64 { + self.n_counted.load(Ordering::Relaxed) + } + + pub fn gene_count(&self, g: usize) -> u32 { + self.counts[g].load(Ordering::Relaxed) + } + + /// Increment the coverage of gene `g` over the genomic interval + /// `[start, end)`, restricted to the gene's merged exons. + pub fn add_block(&self, g: usize, start: u64, end: u64, ann: &GeneAnnotation) { + let exons = &ann.gene_exons[g]; + // Transcript offset of the start of each exon. + let mut tx_off: u64 = 0; + for &(es, ee) in exons { + let exon_len = ee - es; + if end > es && start < ee { + let s = start.max(es); + let e = end.min(ee); + let base = self.offsets[g] + tx_off + (s - es); + for i in 0..(e - s) { + self.cov[(base + i) as usize].fetch_add(1, Ordering::Relaxed); + } + } + tx_off += exon_len; + } + } + + /// Read back one gene's coverage vector (test / output helper). + pub fn gene_slice(&self, g: usize) -> Vec { + (self.offsets[g]..self.offsets[g + 1]) + .map(|i| self.cov[i as usize].load(Ordering::Relaxed)) + .collect() + } +} +``` + +Note the test calls `add_block(g, start, end)` with three arguments while the implementation takes `ann` as well; update the tests to pass `&ann` (the annotation is needed to walk the exons). Adjust the three test bodies to `cov.add_block(0, 102, 105, &ann);` and so on. + +Add to `src/quant/mod.rs`, immediately after the module doc comment / first `use` block: + +```rust +pub mod coverage; +``` + +Also add `#[derive(Default)]`-compatible construction: the test builds `GeneAnnotation` with `..Default::default()`, which already works because `GeneAnnotation` derives `Default`. + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --lib quant::coverage` +Expected: 3 passed. + +- [ ] **Step 5: Lint and commit** + +```bash +cargo fmt +cargo clippy --all-targets 2>&1 | tail -5 +git add src/quant/coverage.rs src/quant/mod.rs +git commit -m "feat(degnorm): per-gene transcript-space coverage accumulator" +``` + +--- + +### Task 2: Read-to-coverage rule (SE and PE) + +**Files:** +- Modify: `src/quant/coverage.rs` + +**Interfaces:** +- Consumes: `GeneCoverage::add_block`, `Transcript { exons: Vec, is_reverse }`, `Exon { genome_start, genome_end }`, `PairedAlignment { mate1_transcript, mate2_transcript }`, `GeneAnnotation::overlapping_genes_into`. +- Produces: `GeneCoverage::count_se_read(&self, transcripts: &[Transcript], ann: &GeneAnnotation)`, `GeneCoverage::count_pe_read(&self, both_mapped: &[&PairedAlignment], ann: &GeneAnnotation)`. + +- [ ] **Step 1: Write the failing test** + +Add to the `tests` module in `src/quant/coverage.rs`: + +```rust + use crate::align::transcript::{Exon, Transcript}; + + fn tr(blocks: &[(u64, u64)]) -> Transcript { + let mut t = Transcript::default(); + t.chr_idx = 0; + t.genome_start = blocks[0].0; + t.genome_end = blocks[blocks.len() - 1].1; + t.exons = blocks + .iter() + .map(|&(s, e)| Exon { + genome_start: s, + genome_end: e, + read_start: 0, + read_end: (e - s) as usize, + i_frag: 0, + }) + .collect(); + t + } + + #[test] + fn multimapper_is_not_counted() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.count_se_read(&[tr(&[(100, 105)]), tr(&[(100, 105)])], &ann); + assert_eq!(cov.gene_slice(0).iter().sum::(), 0); + assert_eq!(cov.gene_count(0), 0); + assert_eq!(cov.n_counted(), 0); + } + + #[test] + fn unique_read_increments_coverage_and_count() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.count_se_read(&[tr(&[(100, 105)])], &ann); + assert_eq!(cov.gene_slice(0)[..5], [1, 1, 1, 1, 1]); + assert_eq!(cov.gene_count(0), 1); + assert_eq!(cov.n_counted(), 1); + } + + #[test] + fn overlapping_mates_are_counted_once_per_base() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + let pair = crate::align::read_align::PairedAlignment { + mate1_transcript: tr(&[(100, 106)]), + mate2_transcript: tr(&[(104, 110)]), + mate1_region: (0, 6), + mate2_region: (0, 6), + is_proper_pair: true, + insert_size: 10, + combined_wt_score: 0, + combined_n_match: 12, + }; + cov.count_pe_read(&[&pair], &ann); + // Bases 100..110 covered exactly once each, including the 104..106 overlap. + assert_eq!(cov.gene_slice(0)[..10], [1; 10]); + assert_eq!(cov.gene_count(0), 1); + } +``` + +If `Transcript` does not implement `Default`, build it with the struct literal used by the existing tests in `src/quant/mod.rs` (`make_transcript`) instead; copy that helper rather than adding a `Default` impl. + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --lib quant::coverage` +Expected: FAIL, `count_se_read` / `count_pe_read` undefined. + +- [ ] **Step 3: Write minimal implementation** + +Add to `impl GeneCoverage` in `src/quant/coverage.rs`: + +```rust + /// Single-end: accumulate coverage for a uniquely mapped read assigned to + /// exactly one gene. Mirrors the `GeneCounts` column-1 rule so counts and + /// coverage always agree; multimappers and gene-ambiguous reads are skipped + /// (DegNorm's default: unique alignments only). + pub fn count_se_read(&self, transcripts: &[Transcript], ann: &GeneAnnotation) { + if transcripts.len() != 1 { + return; + } + let t = &transcripts[0]; + let mut genes = Vec::new(); + ann.overlapping_genes_into(t, &mut genes); + if genes.len() != 1 { + return; + } + let g = genes[0]; + for ex in &t.exons { + self.add_block(g, ex.genome_start, ex.genome_end, ann); + } + self.counts[g].fetch_add(1, Ordering::Relaxed); + self.n_counted.fetch_add(1, Ordering::Relaxed); + } + + /// Paired-end: the fragment is one observation. Mate blocks are merged + /// before accumulation so an overlapping pair contributes 1 per base, which + /// is what DegNorm's paired-read coverage does. + pub fn count_pe_read(&self, both_mapped: &[&PairedAlignment], ann: &GeneAnnotation) { + if both_mapped.len() != 1 { + return; + } + let pair = both_mapped[0]; + let mut genes = Vec::new(); + let mut genes2 = Vec::new(); + ann.overlapping_genes_into(&pair.mate1_transcript, &mut genes); + ann.overlapping_genes_into(&pair.mate2_transcript, &mut genes2); + genes.extend_from_slice(&genes2); + genes.sort_unstable(); + genes.dedup(); + if genes.len() != 1 { + return; + } + let g = genes[0]; + + // Merge both mates' blocks so overlapping bases are counted once. + let mut blocks: Vec<(u64, u64)> = pair + .mate1_transcript + .exons + .iter() + .chain(pair.mate2_transcript.exons.iter()) + .map(|e| (e.genome_start, e.genome_end)) + .collect(); + blocks.sort_unstable(); + let mut merged: Vec<(u64, u64)> = Vec::with_capacity(blocks.len()); + for (s, e) in blocks { + match merged.last_mut() { + Some(last) if s <= last.1 => last.1 = last.1.max(e), + _ => merged.push((s, e)), + } + } + for (s, e) in merged { + self.add_block(g, s, e, ann); + } + self.counts[g].fetch_add(1, Ordering::Relaxed); + self.n_counted.fetch_add(1, Ordering::Relaxed); + } +``` + +Add the imports at the top of the file: + +```rust +use crate::align::read_align::PairedAlignment; +use crate::align::transcript::Transcript; +``` + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --lib quant::coverage` +Expected: 6 passed. + +- [ ] **Step 5: Lint and commit** + +```bash +cargo fmt && cargo clippy --all-targets 2>&1 | tail -5 +git add src/quant/coverage.rs +git commit -m "feat(degnorm): SE and PE read-to-coverage rules with mate merging" +``` + +--- + +### Task 3: `GeneCoverage.out.bin` writer and reader + +**Files:** +- Modify: `src/quant/coverage.rs` + +**Interfaces:** +- Consumes: `GeneCoverage`, `GeneAnnotation.gene_ids`, `Genome` (unused; chr index comes from the annotation build). +- Produces: + - `GeneCoverage::write_file(&self, path: &Path, ann: &GeneAnnotation, sample_id: &str, paired: bool) -> Result<(), Error>` + - `pub struct CoverageFile { pub sample_id: String, pub paired: bool, pub n_counted: u64, pub gene_ids: Vec, pub gene_lens: Vec, pub counts: Vec, pub cov: Vec, pub offsets: Vec }` + - `CoverageFile::read(path: &Path) -> Result` + - `CoverageFile::gene(&self, g: usize) -> &[u32]` + +- [ ] **Step 1: Write the failing test** + +Add to the `tests` module: + +```rust + #[test] + fn coverage_file_round_trips() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.count_se_read(&[tr(&[(100, 105)])], &ann); + + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("GeneCoverage.out.bin"); + cov.write_file(&path, &ann, "sampleA", false).unwrap(); + + let f = CoverageFile::read(&path).unwrap(); + assert_eq!(f.sample_id, "sampleA"); + assert!(!f.paired); + assert_eq!(f.n_counted, 1); + assert_eq!(f.gene_ids, vec!["G1".to_string()]); + assert_eq!(f.gene_lens, vec![15]); + assert_eq!(f.counts, vec![1]); + assert_eq!(f.gene(0), &cov.gene_slice(0)[..]); + } + + #[test] + fn bad_magic_is_an_error_not_a_panic() { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("junk.bin"); + std::fs::write(&path, b"not a coverage file").unwrap(); + assert!(CoverageFile::read(&path).is_err()); + } +``` + +`tempfile` is a dev-dependency, so use `#[cfg(test)] use tempfile;` implicitly (it is available in unit tests without a `use`). + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --lib quant::coverage` +Expected: FAIL, `write_file` / `CoverageFile` undefined. + +- [ ] **Step 3: Write minimal implementation** + +Add to `src/quant/coverage.rs`: + +```rust +use std::io::{BufReader, BufWriter, Read, Write}; +use std::path::Path; + +use flate2::Compression; +use flate2::read::GzDecoder; +use flate2::write::GzEncoder; + +use crate::error::Error; + +const COV_MAGIC: &[u8; 8] = b"RSDGNCOV"; +const COV_VERSION: u32 = 1; + +impl GeneCoverage { + /// Write `GeneCoverage.out.bin`: a gzip stream holding a fixed header, a + /// per-gene table, the gene id block, and the flat coverage array. + pub fn write_file( + &self, + path: &Path, + ann: &GeneAnnotation, + sample_id: &str, + paired: bool, + ) -> Result<(), Error> { + let file = std::fs::File::create(path).map_err(|e| Error::io(e, path))?; + let mut w = GzEncoder::new(BufWriter::new(file), Compression::new(6)); + let io = |e: std::io::Error| Error::io(e, path); + + w.write_all(COV_MAGIC).map_err(io)?; + w.write_all(&COV_VERSION.to_le_bytes()).map_err(io)?; + w.write_all(&(paired as u32).to_le_bytes()).map_err(io)?; + let n = self.n_genes() as u64; + w.write_all(&n.to_le_bytes()).map_err(io)?; + w.write_all(&self.total_len().to_le_bytes()).map_err(io)?; + w.write_all(&self.n_counted().to_le_bytes()).map_err(io)?; + + let sid = sample_id.as_bytes(); + w.write_all(&(sid.len() as u16).to_le_bytes()).map_err(io)?; + w.write_all(sid).map_err(io)?; + + for g in 0..self.n_genes() { + w.write_all(&self.gene_len(g).to_le_bytes()).map_err(io)?; + w.write_all(&self.gene_count(g).to_le_bytes()).map_err(io)?; + } + for id in &ann.gene_ids { + let b = id.as_bytes(); + w.write_all(&(b.len() as u16).to_le_bytes()).map_err(io)?; + w.write_all(b).map_err(io)?; + } + for c in &self.cov { + w.write_all(&c.load(Ordering::Relaxed).to_le_bytes()) + .map_err(io)?; + } + w.finish().map_err(io)?; + Ok(()) + } +} + +/// A loaded `GeneCoverage.out.bin`. +pub struct CoverageFile { + pub sample_id: String, + pub paired: bool, + pub n_counted: u64, + pub gene_ids: Vec, + pub gene_lens: Vec, + pub counts: Vec, + pub offsets: Vec, + pub cov: Vec, +} + +impl CoverageFile { + pub fn read(path: &Path) -> Result { + let file = std::fs::File::open(path).map_err(|e| Error::io(e, path))?; + let mut r = GzDecoder::new(BufReader::new(file)); + let io = |e: std::io::Error| Error::io(e, path); + + let mut magic = [0u8; 8]; + r.read_exact(&mut magic).map_err(io)?; + if &magic != COV_MAGIC { + return Err(Error::Parameter(format!( + "{} is not a rustar-aligner GeneCoverage file", + path.display() + ))); + } + let mut u32buf = [0u8; 4]; + let mut u64buf = [0u8; 8]; + let mut u16buf = [0u8; 2]; + + r.read_exact(&mut u32buf).map_err(io)?; + let version = u32::from_le_bytes(u32buf); + if version != COV_VERSION { + return Err(Error::Parameter(format!( + "{} has GeneCoverage version {version}, expected {COV_VERSION}", + path.display() + ))); + } + r.read_exact(&mut u32buf).map_err(io)?; + let paired = u32::from_le_bytes(u32buf) != 0; + r.read_exact(&mut u64buf).map_err(io)?; + let n_genes = u64::from_le_bytes(u64buf) as usize; + r.read_exact(&mut u64buf).map_err(io)?; + let total_len = u64::from_le_bytes(u64buf) as usize; + r.read_exact(&mut u64buf).map_err(io)?; + let n_counted = u64::from_le_bytes(u64buf); + + r.read_exact(&mut u16buf).map_err(io)?; + let mut sid = vec![0u8; u16::from_le_bytes(u16buf) as usize]; + r.read_exact(&mut sid).map_err(io)?; + let sample_id = String::from_utf8_lossy(&sid).into_owned(); + + let mut gene_lens = Vec::with_capacity(n_genes); + let mut counts = Vec::with_capacity(n_genes); + for _ in 0..n_genes { + r.read_exact(&mut u32buf).map_err(io)?; + gene_lens.push(u32::from_le_bytes(u32buf)); + r.read_exact(&mut u32buf).map_err(io)?; + counts.push(u32::from_le_bytes(u32buf)); + } + let mut gene_ids = Vec::with_capacity(n_genes); + for _ in 0..n_genes { + r.read_exact(&mut u16buf).map_err(io)?; + let mut b = vec![0u8; u16::from_le_bytes(u16buf) as usize]; + r.read_exact(&mut b).map_err(io)?; + gene_ids.push(String::from_utf8_lossy(&b).into_owned()); + } + let mut raw = vec![0u8; total_len * 4]; + r.read_exact(&mut raw).map_err(io)?; + let cov: Vec = raw + .chunks_exact(4) + .map(|c| u32::from_le_bytes([c[0], c[1], c[2], c[3]])) + .collect(); + + let mut offsets = Vec::with_capacity(n_genes + 1); + let mut acc = 0u64; + offsets.push(0); + for &l in &gene_lens { + acc += l as u64; + offsets.push(acc); + } + + Ok(CoverageFile { + sample_id, + paired, + n_counted, + gene_ids, + gene_lens, + counts, + offsets, + cov, + }) + } + + pub fn gene(&self, g: usize) -> &[u32] { + &self.cov[self.offsets[g] as usize..self.offsets[g + 1] as usize] + } + + pub fn n_genes(&self) -> usize { + self.gene_ids.len() + } +} +``` + +Check the exact constructor name for parameter errors in `src/error.rs` (the enum has a `Parameter` variant per CLAUDE.md); if it is `Error::parameter(...)` or takes a different shape, match the existing call sites in `src/params/mod.rs`. + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --lib quant::coverage` +Expected: 8 passed. + +- [ ] **Step 5: Lint and commit** + +```bash +cargo fmt && cargo clippy --all-targets 2>&1 | tail -5 +git add src/quant/coverage.rs +git commit -m "feat(degnorm): GeneCoverage.out.bin writer and reader" +``` + +--- + +### Task 4: Wire `--quantMode GeneCoverage` into the aligner + +**Files:** +- Modify: `src/params/mod.rs` (accessor next to `quant_gene_counts` around line 1768; validation next to the GeneCounts GTF check around line 1484; new `--degNormSampleId` flag near the other quant flags) +- Modify: `src/quant/mod.rs` (`QuantContext` gains an optional coverage field) +- Modify: `src/lib.rs` (context build ~line 274, SE call site ~line 1835, PE call site ~line 3207, output write ~line 399) + +**Interfaces:** +- Consumes: `GeneCoverage`, `CoverageFile` from Tasks 1-3. +- Produces: `Parameters::quant_gene_coverage(&self) -> bool`, `QuantContext.coverage: Option`. + +- [ ] **Step 1: Write the failing test** + +Add to the params test module in `src/params/mod.rs`: + +```rust + #[test] + fn quant_mode_gene_coverage_is_recognised() { + let p = try_parse(&[ + "--readFilesIn", + "r.fq", + "--quantMode", + "GeneCoverage", + "--sjdbGTFfile", + "a.gtf", + ]) + .unwrap(); + assert!(p.quant_gene_coverage()); + assert!(!p.quant_gene_counts()); + } + + #[test] + fn quant_mode_gene_coverage_requires_gtf() { + let err = try_parse(&["--readFilesIn", "r.fq", "--quantMode", "GeneCoverage"]); + assert!(err.is_err()); + } +``` + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --lib params::tests::quant_mode_gene_coverage` +Expected: FAIL, `quant_gene_coverage` undefined. + +- [ ] **Step 3: Write minimal implementation** + +In `src/params/mod.rs`, next to `quant_gene_counts`: + +```rust + /// Returns true if `--quantMode GeneCoverage` was requested. + /// + /// Not a STAR mode: this is the rustar-aligner extension that captures + /// per-gene exonic coverage for the DegNorm pipeline (`--runMode degNorm`). + pub fn quant_gene_coverage(&self) -> bool { + self.quant_mode.iter().any(|m| m == "GeneCoverage") + } +``` + +Extend the existing GeneCounts GTF validation so it also fires for `GeneCoverage`: + +```rust + // quantMode GeneCounts / GeneCoverage require a GTF file + if (params.quant_gene_counts() || params.quant_gene_coverage()) + && params.sjdb_gtf_file.is_none() + { + return Err(Error::Parameter( + "--quantMode GeneCounts/GeneCoverage requires --sjdbGTFfile".to_string(), + )); + } +``` + +(Keep the exact wording of the existing message for GeneCounts-only runs if other tests assert on it; check with `grep -rn "requires --sjdbGTFfile" src tests` before editing.) + +Add the sample-id flag near the other quant params: + +```rust + /// Sample name recorded in `GeneCoverage.out.bin`. Defaults to the + /// basename of `--outFileNamePrefix`, or `sample` if that is empty. + #[arg(long = "degNormSampleId", default_value = "")] + pub deg_norm_sample_id: String, +``` + +In `src/quant/mod.rs`, extend `QuantContext`: + +```rust +pub struct QuantContext { + pub gene_ann: GeneAnnotation, + pub counts: GeneCounts, + /// Present only when `--quantMode GeneCoverage` was requested. + pub coverage: Option, +} +``` + +and give `build` a `with_coverage: bool` argument that sets it: + +```rust + let coverage = if with_coverage { + Some(coverage::GeneCoverage::new(&gene_ann)) + } else { + None + }; + Ok(QuantContext { gene_ann, counts, coverage }) +``` + +Update the existing `QuantContext::build` call in `src/lib.rs` (~line 281) to pass `params.quant_gene_coverage()`, and change the guard at line 274 from `params.quant_gene_counts()` to `params.quant_gene_counts() || params.quant_gene_coverage()`. Any other `QuantContext::build` call sites (grep for them, including solo) pass `false`. + +At the SE call site in `src/lib.rs` (~line 1835), just below the `count_se_read` call: + +```rust + if let Some(ref cov) = q.coverage { + cov.count_se_read(&transcripts, &q.gene_ann); + } +``` + +At the PE call site (~line 3207), below `count_pe_read`: + +```rust + if let Some(ref cov) = q.coverage { + cov.count_pe_read(&bm_deref, &q.gene_ann); + } +``` + +Guard the `ReadsPerGene.out.tab` write so it only runs for GeneCounts, and add the coverage write (~line 399): + +```rust + if let Some(ref ctx) = quant_ctx { + if params.quant_gene_counts() { + let quant_path = params.output_path("ReadsPerGene.out.tab"); + ctx.counts.write_output(&quant_path, &ctx.gene_ann)?; + info!("Wrote {}", quant_path.display()); + } + if let Some(ref cov) = ctx.coverage { + let cov_path = params.output_path("GeneCoverage.out.bin"); + let sample_id = params.deg_norm_sample_id_or_default(); + cov.write_file(&cov_path, &ctx.gene_ann, &sample_id, params.is_paired_end())?; + info!("Wrote {}", cov_path.display()); + } + } +``` + +Add the helper to `src/params/mod.rs` (check the actual name of the paired-end predicate with `grep -n "fn is_paired\|read_files_in.len() > 1" src/params/mod.rs` and use it): + +```rust + /// Sample id for the DegNorm coverage file: `--degNormSampleId` when set, + /// otherwise the trimmed `--outFileNamePrefix`, otherwise `sample`. + pub fn deg_norm_sample_id_or_default(&self) -> String { + if !self.deg_norm_sample_id.is_empty() { + return self.deg_norm_sample_id.clone(); + } + let p = self.out_file_name_prefix.trim_end_matches('/'); + let base = p.rsplit('/').next().unwrap_or(""); + if base.is_empty() { + "sample".to_string() + } else { + base.trim_end_matches('.').to_string() + } + } +``` + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --lib params && cargo test --lib quant` +Expected: PASS. + +- [ ] **Step 5: Verify alignment output is unchanged and commit** + +```bash +cargo build --release +cargo test +cargo fmt && cargo clippy --all-targets 2>&1 | tail -5 +git add src/params/mod.rs src/quant/mod.rs src/lib.rs +git commit -m "feat(degnorm): --quantMode GeneCoverage writes GeneCoverage.out.bin" +``` + +--- + +### Task 5: Rank-one approximation and NMF-OA + +**Files:** +- Create: `src/degnorm/mod.rs` (module declarations only for now: `pub mod nmf;`) +- Create: `src/degnorm/nmf.rs` +- Modify: `src/lib.rs` (add `pub mod degnorm;` next to the other module declarations) + +**Interfaces:** +- Produces: + - `pub struct Mat { pub p: usize, pub l: usize, pub data: Vec }` with `Mat::new(p, l)`, `Mat::get(i, j)`, `Mat::set(i, j, v)`, `Mat::row_sums() -> Vec`, `Mat::col_max() -> Vec`, `Mat::max() -> f64` + - `pub fn rank_one(x: &Mat) -> (Vec, Vec)` returning `(k, e)` with `k.len() == p`, `e.len() == l`, such that `k[i] * e[j]` approximates `x` + - `pub fn nmf_oa(x: &Mat, iters: usize) -> (Vec, Vec)` + - `pub fn ratio_svd(x: &Mat) -> Mat` (rank-one estimate, raised elementwise to at least `x`) + +- [ ] **Step 1: Write the failing test** + +Create `src/degnorm/nmf.rs` with the test module: + +```rust +#[cfg(test)] +mod tests { + use super::*; + + fn mat_from(rows: &[&[f64]]) -> Mat { + let p = rows.len(); + let l = rows[0].len(); + let mut m = Mat::new(p, l); + for (i, r) in rows.iter().enumerate() { + for (j, &v) in r.iter().enumerate() { + m.set(i, j, v); + } + } + m + } + + #[test] + fn rank_one_recovers_a_planted_rank_one_matrix() { + // x = k e^T with k = [1, 2, 3], e = [4, 5] + let x = mat_from(&[&[4.0, 5.0], &[8.0, 10.0], &[12.0, 15.0]]); + let (k, e) = rank_one(&x); + for i in 0..3 { + for j in 0..2 { + assert!((k[i] * e[j] - x.get(i, j)).abs() < 1e-8); + } + } + assert!(k.iter().all(|&v| v > 0.0)); + assert!(e.iter().all(|&v| v > 0.0)); + } + + #[test] + fn nmf_oa_over_approximates() { + // Sample 2 is degraded at the 3' end: rank-one cannot fit it exactly. + let x = mat_from(&[ + &[10.0, 10.0, 10.0, 10.0], + &[10.0, 10.0, 2.0, 1.0], + ]); + let (k, e) = nmf_oa(&x, 100); + let mut n_below = 0; + for i in 0..x.p { + for j in 0..x.l { + if k[i] * e[j] < x.get(i, j) - 1e-6 { + n_below += 1; + } + } + } + assert_eq!(n_below, 0, "NMF-OA estimate must dominate the input"); + } + + #[test] + fn ratio_svd_dominates_input() { + let x = mat_from(&[&[5.0, 1.0], &[1.0, 5.0]]); + let est = ratio_svd(&x); + for i in 0..2 { + for j in 0..2 { + assert!(est.get(i, j) >= x.get(i, j) - 1e-9); + } + } + } +} +``` + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --lib degnorm::nmf` +Expected: FAIL, module not declared / `Mat` undefined. + +- [ ] **Step 3: Write minimal implementation** + +Above the tests in `src/degnorm/nmf.rs`: + +```rust +//! Rank-one NMF over-approximation, ported from DegNorm's `degnorm/nmf.py` +//! (NUStatBioinfo/DegNorm). Pure math: no I/O, no globals. +//! +//! DegNorm calls `scipy.sparse.linalg.svds(x, k=1)`. Sample counts are tiny +//! (p is the number of RNA-seq libraries), so the leading singular triplet is +//! obtained here by power iteration on the p x p Gram matrix `x x^T`, which is +//! both cheaper and deterministic. + +/// A dense p x l matrix in row-major order. +#[derive(Clone, Debug)] +pub struct Mat { + pub p: usize, + pub l: usize, + pub data: Vec, +} + +impl Mat { + pub fn new(p: usize, l: usize) -> Self { + Mat { p, l, data: vec![0.0; p * l] } + } + + #[inline] + pub fn get(&self, i: usize, j: usize) -> f64 { + self.data[i * self.l + j] + } + + #[inline] + pub fn set(&mut self, i: usize, j: usize, v: f64) { + self.data[i * self.l + j] = v; + } + + pub fn row_sums(&self) -> Vec { + (0..self.p) + .map(|i| self.data[i * self.l..(i + 1) * self.l].iter().sum()) + .collect() + } + + pub fn col_max(&self) -> Vec { + (0..self.l) + .map(|j| (0..self.p).map(|i| self.get(i, j)).fold(f64::MIN, f64::max)) + .collect() + } + + pub fn max(&self) -> f64 { + self.data.iter().copied().fold(f64::MIN, f64::max) + } + + /// Keep only the columns listed in `idx`, in order. + pub fn select_cols(&self, idx: &[usize]) -> Mat { + let mut out = Mat::new(self.p, idx.len()); + for i in 0..self.p { + for (jj, &j) in idx.iter().enumerate() { + out.set(i, jj, self.get(i, j)); + } + } + out + } +} + +/// Leading singular triplet: returns `(k, e)` with `k[i]*e[j] ~= x[i][j]`, +/// signs fixed so both vectors are non-negative for a non-negative input. +pub fn rank_one(x: &Mat) -> (Vec, Vec) { + let p = x.p; + // Gram matrix G = x x^T (p x p). + let mut g = vec![0.0f64; p * p]; + for a in 0..p { + for b in a..p { + let mut s = 0.0; + for j in 0..x.l { + s += x.get(a, j) * x.get(b, j); + } + g[a * p + b] = s; + g[b * p + a] = s; + } + } + // Power iteration for the dominant eigenvector. + let mut u = vec![1.0 / (p as f64).sqrt(); p]; + for _ in 0..1000 { + let mut next = vec![0.0f64; p]; + for a in 0..p { + let mut s = 0.0; + for b in 0..p { + s += g[a * p + b] * u[b]; + } + next[a] = s; + } + let norm = next.iter().map(|v| v * v).sum::().sqrt(); + if norm <= 0.0 { + return (vec![0.0; p], vec![0.0; x.l]); + } + for v in next.iter_mut() { + *v /= norm; + } + let delta: f64 = next.iter().zip(&u).map(|(a, b)| (a - b).abs()).sum(); + u = next; + if delta < 1e-12 { + break; + } + } + // Flip so the dominant direction is positive. + if u.iter().sum::() < 0.0 { + for v in u.iter_mut() { + *v = -*v; + } + } + // e_raw = x^T u; s = ||e_raw||; k = u * s; e = e_raw / s. + let mut e: Vec = (0..x.l) + .map(|j| (0..p).map(|i| x.get(i, j) * u[i]).sum()) + .collect(); + let s = e.iter().map(|v| v * v).sum::().sqrt(); + if s <= 0.0 { + return (vec![0.0; p], vec![0.0; x.l]); + } + for v in e.iter_mut() { + *v /= s; + } + let k: Vec = u.iter().map(|v| v * s).collect(); + (k, e) +} + +/// Elementwise product of the rank-one factors into a full matrix. +pub fn outer(k: &[f64], e: &[f64]) -> Mat { + let mut m = Mat::new(k.len(), e.len()); + for i in 0..k.len() { + for j in 0..e.len() { + m.set(i, j, k[i] * e[j]); + } + } + m +} + +/// NMF over-approximation. Faithful port of `GeneNMFOA.nmf`: dual ascent on the +/// non-negative multiplier `lambda` with step `1/sqrt(iters)`. +pub fn nmf_oa(x: &Mat, iters: usize) -> (Vec, Vec) { + let (mut k, mut e) = rank_one(x); + let mut est = outer(&k, &e); + let mut lambda = Mat::new(x.p, x.l); + let c = 1.0 / (iters as f64).sqrt(); + + for _ in 0..iters { + for i in 0..x.p { + for j in 0..x.l { + let res = est.get(i, j) - x.get(i, j); + let v = lambda.get(i, j) - c * res; + lambda.set(i, j, if v < 0.0 { 0.0 } else { v }); + } + } + let mut shifted = Mat::new(x.p, x.l); + for i in 0..x.p { + for j in 0..x.l { + shifted.set(i, j, x.get(i, j) + lambda.get(i, j)); + } + } + let (k2, e2) = rank_one(&shifted); + k = k2; + e = e2; + est = outer(&k, &e); + } + (k, e) +} + +/// One-shot rank-one estimate raised to at least the input (DegNorm's +/// `ratio_svd`, used only for initialisation). +pub fn ratio_svd(x: &Mat) -> Mat { + let (k, e) = rank_one(x); + let mut est = outer(&k, &e); + for i in 0..x.p { + for j in 0..x.l { + if est.get(i, j) < x.get(i, j) { + est.set(i, j, x.get(i, j)); + } + } + } + est +} +``` + +Note: `nmf_oa`'s test asserts the estimate dominates the input. DegNorm's `nmf` does not itself clamp (the clamping line is commented out upstream), so if the test fails, clamp in the *caller* and change the test to assert on `ratio_svd`-style clamped output instead of weakening the port. Record which choice was made in the module doc comment. + +Create `src/degnorm/mod.rs`: + +```rust +//! DegNorm-style transcript degradation normalization (`--runMode degNorm`). +//! +//! Not a STAR feature; see `DIVERGENCE.md`. + +pub mod nmf; +``` + +Add `pub mod degnorm;` to `src/lib.rs` beside the other module declarations. + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --lib degnorm::nmf` +Expected: 3 passed. + +- [ ] **Step 5: Lint and commit** + +```bash +cargo fmt && cargo clippy --all-targets 2>&1 | tail -5 +git add src/degnorm/ src/lib.rs +git commit -m "feat(degnorm): rank-one NMF over-approximation core" +``` + +--- + +### Task 6: Baseline selection and per-gene DI + +**Files:** +- Create: `src/degnorm/baseline.rs` +- Modify: `src/degnorm/mod.rs` (add `pub mod baseline;`) + +**Interfaces:** +- Consumes: `Mat`, `rank_one`, `nmf_oa`, `outer` from Task 5. +- Produces: + - `pub struct BaselineParams { pub nmf_iter: usize, pub bins: usize, pub min_high_coverage: usize, pub downsample_rate: usize, pub skip: bool }` + - `pub struct BaselineResult { pub rho: Vec, pub ran_baseline: bool }` + - `pub fn baseline_selection(f: &Mat, params: &BaselineParams, seed_offset: u64) -> BaselineResult` + +- [ ] **Step 1: Write the failing test** + +Create `src/degnorm/baseline.rs` with: + +```rust +#[cfg(test)] +mod tests { + use super::*; + use crate::degnorm::nmf::Mat; + + fn params() -> BaselineParams { + BaselineParams { + nmf_iter: 100, + bins: 20, + min_high_coverage: 50, + downsample_rate: 1, + skip: false, + } + } + + /// Flat coverage in both samples: no degradation, DI ~ 0. + #[test] + fn undegraded_gene_has_near_zero_di() { + let l = 400; + let mut f = Mat::new(2, l); + for j in 0..l { + f.set(0, j, 20.0); + f.set(1, j, 20.0); + } + let out = baseline_selection(&f, ¶ms(), 0); + assert!(out.rho.iter().all(|&r| r.abs() < 0.05), "rho = {:?}", out.rho); + } + + /// Sample 1 loses coverage over the last 60% of the transcript. + #[test] + fn degraded_sample_gets_higher_di_than_control() { + let l = 400; + let mut f = Mat::new(2, l); + for j in 0..l { + f.set(0, j, 20.0); + f.set(1, j, if j < 160 { 20.0 } else { 2.0 }); + } + let out = baseline_selection(&f, ¶ms(), 0); + assert!( + out.rho[1] > out.rho[0] + 0.2, + "degraded sample DI {} should exceed control DI {}", + out.rho[1], + out.rho[0] + ); + assert!(out.rho.iter().all(|&r| (0.0..=0.9).contains(&r))); + } + + /// Too few high-coverage positions: DegNorm returns all-zero DI. + #[test] + fn short_gene_returns_zero_di() { + let mut f = Mat::new(2, 10); + for j in 0..10 { + f.set(0, j, 5.0); + f.set(1, j, 5.0); + } + let out = baseline_selection(&f, ¶ms(), 0); + assert_eq!(out.rho, vec![0.0, 0.0]); + assert!(!out.ran_baseline); + } +} +``` + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --lib degnorm::baseline` +Expected: FAIL, `baseline_selection` undefined. + +- [ ] **Step 3: Write minimal implementation** + +Above the tests in `src/degnorm/baseline.rs`: + +```rust +//! Baseline selection: DegNorm's per-gene search for a transcript region where +//! degradation is minimal, so the coverage envelope (and therefore the DI +//! scores) can be estimated from undegraded positions. Port of +//! `GeneNMFOA.baseline_selection` in `degnorm/nmf.py`, quirks included: the +//! `+1` in the DI denominator, the 0.1 / 0.2 / 0.9 thresholds, and the +//! bin-dropping loop. + +use crate::degnorm::nmf::{Mat, nmf_oa, outer, rank_one}; + +pub struct BaselineParams { + pub nmf_iter: usize, + pub bins: usize, + pub min_high_coverage: usize, + pub downsample_rate: usize, + pub skip: bool, +} + +pub struct BaselineResult { + pub rho: Vec, + pub ran_baseline: bool, +} + +/// Columns whose sample-wise maximum exceeds 10% of the matrix maximum. +fn high_coverage_idx(f: &Mat) -> Vec { + let thresh = 0.1 * f.max(); + let cm = f.col_max(); + (0..f.l).filter(|&j| cm[j] > thresh).collect() +} + +/// Deterministic systematic sample: start offset derived from `seed_offset` +/// rather than a shared RNG, so results do not depend on thread scheduling. +fn systematic_sample(n: usize, take_every: usize, seed_offset: u64) -> Vec { + if take_every <= 1 || take_every >= n { + return (0..n).collect(); + } + let start = (seed_offset % take_every as u64) as usize; + (start..n).step_by(take_every).collect() +} + +/// Split `0..n` into `bins` roughly equal consecutive chunks (DegNorm's +/// `split_into_chunks`: chunk size is `ceil(n / bins)`). +fn split_into_chunks(n: usize, bins: usize) -> Vec> { + let csize = n.div_ceil(bins.max(1)); + let mut out = Vec::new(); + let mut i = 0; + while i < n { + out.push((i..(i + csize).min(n)).collect()); + i += csize; + } + out +} + +/// DI vector: `1 - rowsum(F) / (rowsum(estimate) + 1)`, DegNorm's formula. +fn di_scores(f: &Mat, est: &Mat) -> Vec { + let fs = f.row_sums(); + let es = est.row_sums(); + (0..f.p).map(|i| 1.0 - fs[i] / (es[i] + 1.0)).collect() +} + +fn clamp_rho(rho: &mut [f64]) { + for r in rho.iter_mut() { + if *r > 0.9 { + *r = 0.9; + } + if *r < 0.0 { + *r = 0.0; + } + } +} + +fn over_approximate(est: &mut Mat, f: &Mat) { + for i in 0..f.p { + for j in 0..f.l { + if est.get(i, j) < f.get(i, j) { + est.set(i, j, f.get(i, j)); + } + } + } +} + +pub fn baseline_selection(f: &Mat, params: &BaselineParams, seed_offset: u64) -> BaselineResult { + let zero = BaselineResult { rho: vec![0.0; f.p], ran_baseline: false }; + + let mut hi = high_coverage_idx(f); + if params.downsample_rate > 1 { + let sampled = systematic_sample(f.l, params.downsample_rate, seed_offset); + let set: std::collections::HashSet = sampled.into_iter().collect(); + hi.retain(|j| set.contains(j)); + } + if hi.len() < params.min_high_coverage.max(2) { + return zero; + } + + let f_start = f.select_cols(&hi); + if f_start.row_sums().iter().any(|&s| s <= 0.0) { + return zero; + } + + let (k_start, e_start) = nmf_oa(&f_start, params.nmf_iter); + let mut f_bin = f_start.clone(); + let (mut k, mut e) = (k_start.clone(), e_start.clone()); + let mut ke = outer(&k, &e); + let mut rho = di_scores(&f_bin, &ke); + + // Exclude extreme non-convergence, as upstream does. + let mut one_minus: Vec = rho.iter().map(|r| 1.0 - r).collect(); + one_minus.sort_by(|a, b| a.partial_cmp(b).unwrap()); + if one_minus[one_minus.len() / 2] > 1.0 { + return zero; + } + + let min_gene_len = (200.0 / params.downsample_rate as f64).ceil().max(2.0) as usize; + let min_bins = (params.bins as f64 * 0.2).ceil() as usize; + let mut ran_baseline = false; + + let can_run = hi.len() >= min_gene_len + && rho.iter().copied().fold(f64::MAX, f64::min) <= 0.2 + && !params.skip; + + if can_run { + let mut bin_segs = split_into_chunks(f_bin.l, params.bins); + while rho.iter().copied().fold(f64::MIN, f64::max) > 0.1 { + ran_baseline = true; + + // Per-column worst squared relative residual, averaged per bin. + let res: Vec = (0..f_bin.l) + .map(|j| { + (0..f_bin.p) + .map(|i| { + let d = (ke.get(i, j) - f_bin.get(i, j)) / (f_bin.get(i, j) + 1.0); + d * d + }) + .fold(f64::MIN, f64::max) + }) + .collect(); + let ss_r: Vec = bin_segs + .iter() + .map(|b| b.iter().map(|&j| res[j]).sum::() / b.len() as f64) + .collect(); + if ss_r.iter().copied().fold(f64::MIN, f64::max) == 0.0 { + break; + } + + let drop = ss_r + .iter() + .enumerate() + .max_by(|a, b| a.1.partial_cmp(b.1).unwrap()) + .map(|(i, _)| i) + .unwrap(); + + // Drop that bin's columns and rebuild consecutive bin indices. + let dropped: std::collections::HashSet = + bin_segs[drop].iter().copied().collect(); + let keep: Vec = (0..f_bin.l).filter(|j| !dropped.contains(j)).collect(); + f_bin = f_bin.select_cols(&keep); + bin_segs.remove(drop); + let sizes: Vec = bin_segs.iter().map(|b| b.len()).collect(); + bin_segs = { + let mut out = Vec::with_capacity(sizes.len()); + let mut i = 0; + for s in sizes { + out.push((i..i + s).collect::>()); + i += s; + } + out + }; + + if f_bin.l == 0 || bin_segs.is_empty() { + break; + } + + let (k2, e2) = nmf_oa(&f_bin, params.nmf_iter); + k = k2; + e = e2; + ke = outer(&k, &e); + if ke.row_sums().iter().any(|&s| s == 0.0) { + break; + } + over_approximate(&mut ke, &f_bin); + rho = di_scores(&f_bin, &ke); + + if bin_segs.len() <= min_bins || f_bin.l < min_gene_len { + break; + } + } + + if rho.iter().copied().fold(f64::MIN, f64::max) < 0.2 { + // Refine the envelope on the full high-coverage matrix using the + // baseline-region abundances. + let mut kk: Vec = k.iter().map(|v| v.abs()).collect(); + let min_pos = kk.iter().copied().filter(|&v| v >= 1e-5).fold(f64::MAX, f64::min); + let floor = if min_pos.is_finite() { min_pos } else { 1e-5 }; + for v in kk.iter_mut() { + if *v < 1e-5 { + *v = floor; + } + } + let ee: Vec = (0..f_start.l) + .map(|j| { + (0..f_start.p) + .map(|i| f_start.get(i, j) / kk[i]) + .fold(f64::MIN, f64::max) + }) + .collect(); + let est = outer(&kk, &ee); + rho = di_scores(&f_start, &est); + + if rho.iter().copied().fold(f64::MIN, f64::max) > 0.9 { + let mut est = outer(&k_start, &e_start); + over_approximate(&mut est, &f_start); + rho = di_scores(&f_start, &est); + } + } else { + let mut est = outer(&k_start, &e_start); + over_approximate(&mut est, &f_start); + rho = di_scores(&f_start, &est); + } + } + + clamp_rho(&mut rho); + BaselineResult { rho, ran_baseline } +} +``` + +Add `pub mod baseline;` to `src/degnorm/mod.rs`. + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --lib degnorm::baseline` +Expected: 3 passed. If `undegraded_gene_has_near_zero_di` fails by a small margin, check the `+1` denominators before touching thresholds: with coverage 20 across 400 bases the DI floor is about `1/(8000+1)`, far below 0.05. + +- [ ] **Step 5: Lint and commit** + +```bash +cargo fmt && cargo clippy --all-targets 2>&1 | tail -5 +git add src/degnorm/ +git commit -m "feat(degnorm): baseline selection and per-gene DI scores" +``` + +--- + +### Task 7: DegNorm driver and output tables + +**Files:** +- Create: `src/degnorm/run.rs` +- Modify: `src/degnorm/mod.rs` (add `pub mod run;`) + +**Interfaces:** +- Consumes: `CoverageFile` (Task 3), `Mat`, `ratio_svd`, `baseline_selection`. +- Produces: + - `pub struct DegNormConfig { pub iter: usize, pub nmf_iter: usize, pub downsample_rate: usize, pub minimax_coverage: u32, pub skip_baseline: bool, pub bins: usize, pub min_high_coverage: usize, pub seed: u64 }` + - `pub struct DegNormOutput { pub gene_ids: Vec, pub sample_ids: Vec, pub rho: Vec>, pub raw_counts: Vec>, pub adjusted_counts: Vec>, pub scale_factors: Vec, pub n_baseline_selected: Vec }` + - `pub fn run_degnorm(files: &[CoverageFile], cfg: &DegNormConfig) -> Result` + - `pub fn write_outputs(out: &DegNormOutput, dir: &Path, cfg: &DegNormConfig) -> Result<(), Error>` + +- [ ] **Step 1: Write the failing test** + +Create `src/degnorm/run.rs` with: + +```rust +#[cfg(test)] +mod tests { + use super::*; + use crate::quant::coverage::CoverageFile; + + /// Build an in-memory CoverageFile for `n_genes` genes of length `len`, + /// with per-gene coverage supplied by `f(gene, position)`. + fn synth(sample: &str, n_genes: usize, len: usize, f: impl Fn(usize, usize) -> u32) -> CoverageFile { + let mut cov = Vec::new(); + let mut counts = Vec::new(); + let mut offsets = vec![0u64]; + for g in 0..n_genes { + let mut total = 0u64; + for j in 0..len { + let v = f(g, j); + total += v as u64; + cov.push(v); + } + counts.push((total / 100).max(1) as u32); + offsets.push(offsets[g] + len as u64); + } + CoverageFile { + sample_id: sample.to_string(), + paired: false, + n_counted: counts.iter().map(|&c| c as u64).sum(), + gene_ids: (0..n_genes).map(|g| format!("G{g}")).collect(), + gene_lens: vec![len as u32; n_genes], + counts, + offsets, + cov, + } + } + + fn cfg() -> DegNormConfig { + DegNormConfig { + iter: 2, + nmf_iter: 50, + downsample_rate: 1, + minimax_coverage: 0, + skip_baseline: false, + bins: 20, + min_high_coverage: 50, + seed: 777, + } + } + + #[test] + fn degraded_sample_gets_higher_di_and_upweighted_counts() { + let len = 400; + let a = synth("ctrl", 2, len, |_, _| 20); + let b = synth("degraded", 2, len, |g, j| { + if g == 0 && j >= 160 { 2 } else { 20 } + }); + let out = run_degnorm(&[a, b], &cfg()).unwrap(); + + assert_eq!(out.sample_ids, vec!["ctrl".to_string(), "degraded".to_string()]); + assert_eq!(out.gene_ids.len(), 2); + // Gene 0 is degraded in sample 1 only. + assert!(out.rho[0][1] > out.rho[0][0] + 0.2); + // Adjusted counts scale raw counts up by 1/(1-DI). + assert!(out.adjusted_counts[0][1] > out.raw_counts[0][1]); + } + + #[test] + fn mismatched_gene_sets_are_rejected() { + let a = synth("a", 2, 400, |_, _| 10); + let mut b = synth("b", 2, 400, |_, _| 10); + b.gene_ids[1] = "OTHER".to_string(); + assert!(run_degnorm(&[a, b], &cfg()).is_err()); + } + + #[test] + fn single_sample_is_rejected() { + let a = synth("a", 2, 400, |_, _| 10); + assert!(run_degnorm(&[a], &cfg()).is_err()); + } +} +``` + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --lib degnorm::run` +Expected: FAIL, `run_degnorm` undefined. + +- [ ] **Step 3: Write minimal implementation** + +Above the tests in `src/degnorm/run.rs`: + +```rust +//! DegNorm driver: load per-sample coverage, run the NMF-OA iterations, and +//! write DI scores plus degradation-adjusted counts. Port of +//! `GeneNMFOA.run` in DegNorm's `degnorm/nmf.py`. + +use std::io::Write; +use std::path::Path; + +use rayon::prelude::*; + +use crate::degnorm::baseline::{BaselineParams, baseline_selection}; +use crate::degnorm::nmf::{Mat, ratio_svd}; +use crate::error::Error; +use crate::quant::coverage::CoverageFile; + +pub struct DegNormConfig { + pub iter: usize, + pub nmf_iter: usize, + pub downsample_rate: usize, + pub minimax_coverage: u32, + pub skip_baseline: bool, + pub bins: usize, + pub min_high_coverage: usize, + pub seed: u64, +} + +pub struct DegNormOutput { + pub gene_ids: Vec, + pub sample_ids: Vec, + /// DI score per included gene, per sample. + pub rho: Vec>, + pub raw_counts: Vec>, + pub adjusted_counts: Vec>, + pub scale_factors: Vec, + /// Genes sent through baseline selection, per outer iteration. + pub n_baseline_selected: Vec, +} + +fn median(v: &[f64]) -> f64 { + let mut s = v.to_vec(); + s.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let n = s.len(); + if n == 0 { + return 0.0; + } + if n % 2 == 1 { s[n / 2] } else { (s[n / 2 - 1] + s[n / 2]) / 2.0 } +} + +pub fn run_degnorm(files: &[CoverageFile], cfg: &DegNormConfig) -> Result { + let p = files.len(); + if p < 2 { + return Err(Error::Parameter( + "--runMode degNorm requires at least 2 --degNormCoverageFiles: a degradation \ + index is defined across samples, not within one" + .to_string(), + )); + } + // Validate identical gene models. + for f in &files[1..] { + if f.gene_ids != files[0].gene_ids || f.gene_lens != files[0].gene_lens { + let bad = f + .gene_ids + .iter() + .zip(&files[0].gene_ids) + .position(|(a, b)| a != b) + .map(|i| f.gene_ids[i].clone()) + .unwrap_or_else(|| "".to_string()); + return Err(Error::Parameter(format!( + "coverage file '{}' has a different gene model than '{}' (first difference: {}); \ + all samples must be aligned against the same GTF", + f.sample_id, files[0].sample_id, bad + ))); + } + } + + let sample_ids: Vec = files.iter().map(|f| f.sample_id.clone()).collect(); + + // Select genes: every sample must reach minimax_coverage, and the gene must + // be long enough for a rank-one fit to mean anything. + let mut kept: Vec = Vec::new(); + for g in 0..files[0].n_genes() { + if files[0].gene_lens[g] < 2 { + continue; + } + let ok = files + .iter() + .all(|f| f.gene(g).iter().copied().max().unwrap_or(0) >= cfg.minimax_coverage.max(1)); + if ok { + kept.push(g); + } + } + if kept.is_empty() { + return Err(Error::Parameter( + "no gene passed --degNormMinimaxCoverage in every sample".to_string(), + )); + } + + // Coverage matrices, one per kept gene, shape p x L. + let mats: Vec = kept + .iter() + .map(|&g| { + let l = files[0].gene_lens[g] as usize; + let mut m = Mat::new(p, l); + for (i, f) in files.iter().enumerate() { + for (j, &v) in f.gene(g).iter().enumerate() { + m.set(i, j, v as f64); + } + } + m + }) + .collect(); + + // Raw count matrix (genes x samples). + let x: Vec> = kept + .iter() + .map(|&g| files.iter().map(|f| f.counts[g] as f64).collect()) + .collect(); + let n_genes = kept.len(); + + // Initialisation: single-shot over-approximation to get starting DI scores. + let mut rho: Vec> = mats + .par_iter() + .map(|m| { + let est = ratio_svd(m); + let fs = m.row_sums(); + let es = est.row_sums(); + (0..p).map(|i| 1.0 - fs[i] / (es[i] + 1.0)).collect() + }) + .collect(); + for r in rho.iter_mut() { + for v in r.iter_mut() { + *v = v.clamp(0.0, 0.9); + } + } + + let low_di: Vec = (0..n_genes) + .filter(|&g| rho[g].iter().copied().fold(f64::MIN, f64::max) < 0.1) + .collect(); + let count_sums: Vec = (0..p) + .map(|i| { + if low_di.is_empty() { + (0..n_genes).map(|g| x[g][i]).sum() + } else { + low_di.iter().map(|&g| x[g][i]).sum() + } + }) + .collect(); + let med = median(&count_sums); + let mut norm_factors: Vec = count_sums + .iter() + .map(|&c| if med > 0.0 { c / med } else { 1.0 }) + .collect(); + let mut x_weighted: Vec> = (0..n_genes) + .map(|g| (0..p).map(|i| x[g][i] / norm_factors[i].max(1e-12)).collect()) + .collect(); + let mut scale_factors = norm_factors.clone(); + let mut x_adj = x_weighted.clone(); + let mut n_baseline_selected = Vec::with_capacity(cfg.iter); + + let bp = BaselineParams { + nmf_iter: cfg.nmf_iter, + bins: cfg.bins, + min_high_coverage: cfg.min_high_coverage, + downsample_rate: cfg.downsample_rate, + skip: cfg.skip_baseline, + }; + + for _it in 0..cfg.iter { + // Scale coverage curves by the per-sample depth factors. + let results: Vec<(Vec, bool)> = mats + .par_iter() + .enumerate() + .map(|(gi, m)| { + let mut adj = Mat::new(m.p, m.l); + for i in 0..m.p { + let s = scale_factors[i].max(1e-12); + for j in 0..m.l { + adj.set(i, j, m.get(i, j) / s); + } + } + let r = baseline_selection(&adj, &bp, cfg.seed.wrapping_add(gi as u64)); + (r.rho, r.ran_baseline) + }) + .collect(); + + n_baseline_selected.push(results.iter().filter(|(_, b)| *b).count()); + for (g, (r, _)) in results.into_iter().enumerate() { + rho[g] = r; + } + + // Genes that never went through baseline selection inherit the + // sample-average DI (DegNorm's correct_di_scores). + x_adj = (0..n_genes) + .map(|g| { + (0..p) + .map(|i| x_weighted[g][i] / (1.0 - rho[g][i]).max(1e-12)) + .collect() + }) + .collect(); + let sample_avg_di: Vec = (0..p) + .map(|i| { + let w: f64 = (0..n_genes).map(|g| x_weighted[g][i]).sum(); + let a: f64 = (0..n_genes).map(|g| x_adj[g][i]).sum(); + if a > 0.0 { 1.0 - w / a } else { 0.0 } + }) + .collect(); + for g in 0..n_genes { + if rho[g].iter().copied().fold(f64::MIN, f64::max) == 0.0 { + rho[g] = sample_avg_di.clone(); + } + } + + x_adj = (0..n_genes) + .map(|g| { + (0..p) + .map(|i| x_weighted[g][i] / (1.0 - rho[g][i]).max(1e-12)) + .collect() + }) + .collect(); + + let col_sums: Vec = (0..p).map(|i| (0..n_genes).map(|g| x_adj[g][i]).sum()).collect(); + let med = median(&col_sums); + norm_factors = col_sums + .iter() + .map(|&c| if med > 0.0 { c / med } else { 1.0 }) + .collect(); + for g in 0..n_genes { + for i in 0..p { + x_weighted[g][i] /= norm_factors[i].max(1e-12); + } + } + for i in 0..p { + scale_factors[i] *= norm_factors[i]; + } + } + + Ok(DegNormOutput { + gene_ids: kept.iter().map(|&g| files[0].gene_ids[g].clone()).collect(), + sample_ids, + rho, + raw_counts: x, + adjusted_counts: x_adj, + scale_factors, + n_baseline_selected, + }) +} + +fn write_matrix( + path: &Path, + row_ids: &[String], + header: &[String], + rows: &[Vec], +) -> Result<(), Error> { + let mut f = std::fs::File::create(path).map_err(|e| Error::io(e, path))?; + let io = |e: std::io::Error| Error::io(e, path); + writeln!(f, "gene\t{}", header.join("\t")).map_err(io)?; + for (g, id) in row_ids.iter().enumerate() { + let cells: Vec = rows[g].iter().map(|v| format!("{v:.6}")).collect(); + writeln!(f, "{id}\t{}", cells.join("\t")).map_err(io)?; + } + Ok(()) +} + +pub fn write_outputs(out: &DegNormOutput, dir: &Path, cfg: &DegNormConfig) -> Result<(), Error> { + std::fs::create_dir_all(dir).map_err(|e| Error::io(e, dir))?; + write_matrix(&dir.join("DegradationIndex.tab"), &out.gene_ids, &out.sample_ids, &out.rho)?; + write_matrix(&dir.join("AdjustedCounts.tab"), &out.gene_ids, &out.sample_ids, &out.adjusted_counts)?; + write_matrix(&dir.join("RawCounts.tab"), &out.gene_ids, &out.sample_ids, &out.raw_counts)?; + + let sf_path = dir.join("ScaleFactors.tab"); + let mut f = std::fs::File::create(&sf_path).map_err(|e| Error::io(e, &sf_path))?; + let io = |e: std::io::Error| Error::io(e, sf_path.as_path()); + writeln!(f, "sample\tscaleFactor").map_err(io)?; + for (i, s) in out.sample_ids.iter().enumerate() { + writeln!(f, "{s}\t{:.6}", out.scale_factors[i]).map_err(io)?; + } + + let sum_path = dir.join("Summary.txt"); + let mut f = std::fs::File::create(&sum_path).map_err(|e| Error::io(e, &sum_path))?; + let io = |e: std::io::Error| Error::io(e, sum_path.as_path()); + writeln!(f, "samples\t{}", out.sample_ids.len()).map_err(io)?; + writeln!(f, "genesIncluded\t{}", out.gene_ids.len()).map_err(io)?; + writeln!(f, "degNormIter\t{}", cfg.iter).map_err(io)?; + writeln!(f, "degNormNmfIter\t{}", cfg.nmf_iter).map_err(io)?; + writeln!(f, "downsampleRate\t{}", cfg.downsample_rate).map_err(io)?; + writeln!(f, "skipBaselineSelection\t{}", cfg.skip_baseline).map_err(io)?; + for (i, n) in out.n_baseline_selected.iter().enumerate() { + writeln!(f, "baselineSelectedIter{}\t{}", i + 1, n).map_err(io)?; + } + for (i, s) in out.sample_ids.iter().enumerate() { + let col: Vec = out.rho.iter().map(|r| r[i]).collect(); + writeln!(f, "medianDI_{}\t{:.6}", s, median(&col)).map_err(io)?; + } + Ok(()) +} +``` + +Add `pub mod run;` to `src/degnorm/mod.rs`. `CoverageFile` fields must be public for the test's struct literal; they already are per Task 3. + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --lib degnorm` +Expected: all degnorm tests pass. + +- [ ] **Step 5: Lint and commit** + +```bash +cargo fmt && cargo clippy --all-targets 2>&1 | tail -5 +git add src/degnorm/ +git commit -m "feat(degnorm): NMF-OA driver, DI scores and adjusted count output" +``` + +--- + +### Task 8: `--runMode degNorm` CLI + +**Files:** +- Modify: `src/params/mod.rs` (`RunMode` enum around line 35, `FromStr`, `Display`, new flags, `validate`) +- Modify: `src/lib.rs` (dispatch arm next to `RunMode::LiftOver => liftover::run(params)` at line 81) +- Modify: `src/degnorm/mod.rs` (add `pub fn run(params: &Parameters) -> anyhow::Result<()>`) + +**Interfaces:** +- Consumes: `run_degnorm`, `write_outputs`, `CoverageFile::read`. +- Produces: `RunMode::DegNorm`, `Parameters::deg_norm_*` fields, `degnorm::run`. + +- [ ] **Step 1: Write the failing test** + +Add to the params tests: + +```rust + #[test] + fn run_mode_degnorm_parses_and_validates() { + let p = try_parse(&[ + "--runMode", + "degNorm", + "--degNormCoverageFiles", + "a.bin", + "b.bin", + ]) + .unwrap(); + assert_eq!(p.run_mode(), RunMode::DegNorm); + assert_eq!(p.deg_norm_coverage_files.len(), 2); + assert_eq!(p.deg_norm_iter, 5); + assert_eq!(p.deg_norm_nmf_iter, 100); + } + + #[test] + fn run_mode_degnorm_requires_two_coverage_files() { + assert!(try_parse(&["--runMode", "degNorm", "--degNormCoverageFiles", "a.bin"]).is_err()); + assert!(try_parse(&["--runMode", "degNorm"]).is_err()); + } +``` + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --lib params::tests::run_mode_degnorm` +Expected: FAIL, `RunMode::DegNorm` undefined. + +- [ ] **Step 3: Write minimal implementation** + +In `src/params/mod.rs`, add `DegNorm` to `RunMode`, to `FromStr` (`"degNorm" => Ok(Self::DegNorm)`, and extend the error message to list it), and to `Display` (`Self::DegNorm => write!(f, "degNorm")`). + +Add the flags (near the quant flags): + +```rust + /// `GeneCoverage.out.bin` files, one per sample, for `--runMode degNorm`. + /// At least two are required: a degradation index is defined across samples. + #[arg(long = "degNormCoverageFiles", num_args = 0.., value_name = "FILE")] + pub deg_norm_coverage_files: Vec, + + /// Number of outer DegNorm iterations (DegNorm's `--iter`). + #[arg(long = "degNormIter", default_value_t = 5)] + pub deg_norm_iter: usize, + + /// NMF-OA iterations per gene per outer iteration (DegNorm's `--nmf-iter`). + #[arg(long = "degNormNmfIter", default_value_t = 100)] + pub deg_norm_nmf_iter: usize, + + /// Systematic "take every" downsample rate over transcript positions. + #[arg(long = "degNormDownsampleRate", default_value_t = 1)] + pub deg_norm_downsample_rate: usize, + + /// A gene is included only if its maximum coverage reaches this value in + /// every sample (DegNorm's `--minimax-coverage`). + #[arg(long = "degNormMinimaxCoverage", default_value_t = 0)] + pub deg_norm_minimax_coverage: u32, + + /// Skip the baseline-selection search (faster, less accurate DI scores). + #[arg(long = "degNormSkipBaselineSelection", default_value_t = false)] + pub deg_norm_skip_baseline_selection: bool, + + /// Bins used by baseline selection. + #[arg(long = "degNormBins", default_value_t = 20)] + pub deg_norm_bins: usize, + + /// Minimum number of high-coverage positions for baseline selection. + #[arg(long = "degNormMinHighCoverage", default_value_t = 50)] + pub deg_norm_min_high_coverage: usize, +``` + +In `validate`: + +```rust + if params.run_mode() == RunMode::DegNorm && params.deg_norm_coverage_files.len() < 2 { + return Err(Error::Parameter( + "--runMode degNorm requires at least 2 --degNormCoverageFiles (one per sample); \ + a degradation index is defined across samples, not within one" + .to_string(), + )); + } +``` + +Make sure the `alignReads`-only validations (reads files, genome dir) do not fire for `RunMode::DegNorm`: grep for `run_mode() == RunMode::AlignReads` and for any `else` branch that assumes alignReads, and extend the guards. + +Add to `src/degnorm/mod.rs`: + +```rust +pub mod baseline; +pub mod nmf; +pub mod run; + +use crate::params::Parameters; +use crate::quant::coverage::CoverageFile; + +/// `--runMode degNorm` entry point. +pub fn run_mode(params: &Parameters) -> anyhow::Result<()> { + let mut files = Vec::with_capacity(params.deg_norm_coverage_files.len()); + for p in ¶ms.deg_norm_coverage_files { + log::info!("degNorm: loading {}", p.display()); + files.push(CoverageFile::read(p)?); + } + let cfg = run::DegNormConfig { + iter: params.deg_norm_iter, + nmf_iter: params.deg_norm_nmf_iter, + downsample_rate: params.deg_norm_downsample_rate.max(1), + minimax_coverage: params.deg_norm_minimax_coverage, + skip_baseline: params.deg_norm_skip_baseline_selection, + bins: params.deg_norm_bins, + min_high_coverage: params.deg_norm_min_high_coverage, + seed: params.run_rng_seed, + }; + let out = run::run_degnorm(&files, &cfg)?; + let dir = std::path::PathBuf::from(params.output_path("DegNorm.out")); + run::write_outputs(&out, &dir, &cfg)?; + log::info!( + "degNorm: {} genes x {} samples written to {}", + out.gene_ids.len(), + out.sample_ids.len(), + dir.display() + ); + Ok(()) +} +``` + +Check the actual type returned by `params.output_path(...)` (it is used as `&Path` at line 400 of `src/lib.rs`) and the exact field name of the RNG seed (`grep -n "run_rng_seed" src/params/mod.rs`). + +In `src/lib.rs`, add the dispatch arm: + +```rust + RunMode::DegNorm => degnorm::run_mode(params), +``` + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --lib params && cargo build --release` +Expected: PASS. + +- [ ] **Step 5: Lint and commit** + +```bash +cargo fmt && cargo clippy --all-targets 2>&1 | tail -5 +git add src/params/mod.rs src/lib.rs src/degnorm/mod.rs +git commit -m "feat(degnorm): --runMode degNorm entry point and parameters" +``` + +--- + +### Task 9: End-to-end integration test + +**Files:** +- Create: `tests/degnorm.rs` + +**Interfaces:** +- Consumes: the `rustar-aligner` binary via `assert_cmd`, following the fixture-building pattern in `tests/alignment_features.rs`. + +- [ ] **Step 1: Write the failing test** + +Create `tests/degnorm.rs`. Read `tests/alignment_features.rs` first and reuse its genome/GTF/FASTQ builders verbatim rather than inventing new ones; the sketch below names what the test must assert: + +```rust +//! End-to-end DegNorm: align two synthetic samples with --quantMode +//! GeneCoverage, then merge them with --runMode degNorm. + +mod common; // if tests/alignment_features.rs has no shared module, inline the helpers instead. + +use std::process::Command; + +use assert_cmd::prelude::*; +use tempfile::TempDir; + +#[test] +fn degnorm_end_to_end_flags_the_degraded_sample() { + let dir = TempDir::new().unwrap(); + // 1. Build a 20 kb synthetic genome with two genes, write genome.fa + genes.gtf. + // 2. genomeGenerate into dir/index. + // 3. Sample A: reads spread evenly over both genes. + // Sample B: reads for gene1 drawn only from its 5' 40%, gene2 even. + // 4. Align each with: + // --quantMode GeneCounts GeneCoverage --sjdbGTFfile genes.gtf + // --outFileNamePrefix dir/A_ (and dir/B_) + // 5. Run: + // --runMode degNorm --degNormCoverageFiles dir/A_GeneCoverage.out.bin + // dir/B_GeneCoverage.out.bin --degNormMinHighCoverage 20 + // --outFileNamePrefix dir/ + // 6. Assert dir/DegNorm.out/DegradationIndex.tab exists and parse it: + // gene1's DI in sample B exceeds gene1's DI in sample A by > 0.1, + // and gene2's DI stays below 0.1 in both. + // 7. Assert AdjustedCounts.tab, RawCounts.tab, ScaleFactors.tab, Summary.txt exist. +} + +#[test] +fn gene_coverage_does_not_change_alignment_output() { + // Align the same sample twice, once with --quantMode GeneCounts and once + // with --quantMode GeneCounts GeneCoverage, and assert the two + // Aligned.out.sam files are byte-identical. +} +``` + +Fill in every step with real code before running; a comment-only test body is a plan failure carried into the codebase. + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test --test degnorm` +Expected: FAIL until Tasks 1-8 are wired and the fixture builders are in place. + +- [ ] **Step 3: Implement the fixtures** + +Copy the synthetic-genome, GTF, and FASTQ writers from `tests/alignment_features.rs`. Generate reads with a fixed, non-random layout (fixed stride over the gene body) so the test is deterministic. Use `--degNormMinHighCoverage 20` and `--degNormIter 2` to keep runtime under a few seconds. + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test --test degnorm` +Expected: 2 passed. + +- [ ] **Step 5: Run the whole suite and commit** + +```bash +cargo test +cargo fmt && cargo clippy --all-targets 2>&1 | tail -5 +git add tests/degnorm.rs +git commit -m "test(degnorm): end-to-end coverage capture and DI scoring" +``` + +--- + +### Task 10: Documentation + +**Files:** +- Modify: `DIVERGENCE.md` (section 4) +- Modify: `README.md` (feature list) +- Modify: `CHANGELOG.md` (Unreleased section) +- Modify: `ROADMAP.md` (new phase entry) +- Modify: `CLAUDE.md` (Current Status paragraph, Source Layout tree) +- Create: `docs/src/content/docs/guides/degnorm.md` (match the existing Starlight docs layout; check `ls docs/src/content/docs` first and follow whatever structure is there) + +- [ ] **Step 1: Write the DIVERGENCE.md entry** + +Under section 4, add: + +```markdown +### 4.N DegNorm degradation normalization (`--quantMode GeneCoverage`, `--runMode degNorm`) + +STAR has no equivalent. rustar-aligner can capture per-gene, per-exonic-base +coverage during alignment (`--quantMode GeneCoverage` writes +`GeneCoverage.out.bin`) and merge several samples' coverage into Degradation +Index scores and degradation-adjusted read counts (`--runMode degNorm`), porting +the NMF over-approximation of +[DegNorm](https://nustatbioinfo.github.io/DegNorm/). + +Both are off by default and neither changes alignment output: coverage is +accumulated from the same uniquely mapped, unambiguously assigned reads that +feed `ReadsPerGene.out.tab`. A degradation index is defined across samples, so +`--runMode degNorm` requires at least two coverage files and cannot run inside a +single alignment. +``` + +- [ ] **Step 2: Write the user guide** + +The guide covers: what degradation bias is, the two-phase workflow with a copy-pasteable example (align N samples, then merge), every `--degNorm*` flag with its default and its DegNorm CLI counterpart, the output files with a column description, the memory note (about 4 bytes per exonic base, roughly 280 MB for human GENCODE), and the explicit statement that DI needs at least two samples. + +- [ ] **Step 3: Update README, CHANGELOG, ROADMAP, CLAUDE.md** + +README: one bullet in the feature list linking the guide. CHANGELOG: an `Added` entry naming both flags. ROADMAP: a phase entry describing the two phases and pointing at the spec. CLAUDE.md: append to the Current Status paragraph and add `degnorm/` plus `quant/coverage.rs` to the Source Layout tree. + +- [ ] **Step 4: Verify** + +Run: `cargo test && cargo clippy --all-targets 2>&1 | tail -5` +Expected: green. Check the docs site builds if it has a build step (`ls docs/package.json`). + +- [ ] **Step 5: Commit** + +```bash +git add DIVERGENCE.md README.md CHANGELOG.md ROADMAP.md CLAUDE.md docs/ +git commit -m "docs(degnorm): document GeneCoverage capture and the degNorm run mode" +``` + +--- + +### Task 11: Open the pull request + +- [ ] **Step 1: Full verification** + +```bash +cargo fmt --check +cargo clippy --all-targets +cargo test +``` + +All three must be clean. + +- [ ] **Step 2: Push and open the PR** + +```bash +git push -u origin feat/degnorm +gh pr create --base main --title "feat(degnorm): DegNorm-style degradation normalization" --body "..." +``` + +The body states: what DegNorm is, the two-phase split and why DI cannot be computed inside a single alignment run, the new flags, the outputs, that alignment output is unchanged, and that this is a documented divergence from STAR. End the body with the Claude Code attribution line. diff --git a/docs/superpowers/specs/2026-08-20-degnorm-design.md b/docs/superpowers/specs/2026-08-20-degnorm-design.md new file mode 100644 index 0000000..af7a7ea --- /dev/null +++ b/docs/superpowers/specs/2026-08-20-degnorm-design.md @@ -0,0 +1,294 @@ +# DegNorm in rustar-aligner: design + +Date: 2026-08-20 +Status: approved, awaiting implementation plan + +## 1. Goal + +Reproduce, natively in rustar-aligner, what the +[DegNorm](https://nustatbioinfo.github.io/DegNorm/) pipeline does: correct RNA-seq +read counts for sample- and gene-specific transcript degradation, producing a +Degradation Index (DI) score per gene per sample plus a degradation-adjusted read +count matrix. + +The user request was "do during processing what DegNorm does". The half that can +genuinely move into alignment time is coverage computation; the DI math cannot, +because it is defined across samples (see section 2). + +## 2. What can and cannot be done online + +DegNorm's estimator is a rank-one over-approximation of a gene's coverage matrix +`F` of shape `p x L_i` (p = samples, L_i = exonic length of gene i). With `p = 1` +the rank-one fit is exact and every DI score is 0. A DI score therefore requires +at least two samples and cannot be produced by a single alignment run. + +What a single run *can* produce, essentially for free, is the input DegNorm spends +most of its wall clock computing: the per-gene, per-exonic-base coverage vector, +plus the per-gene raw read count. rustar-aligner already assigns each uniquely +mapped read to a gene for `--quantMode GeneCounts`; adding coverage accumulation +on that same path removes the need to write a sorted BAM, index it, and re-read it +with pysam. + +Split, therefore: + +- **Phase A (alignment time)**: `--quantMode GeneCoverage` emits + `GeneCoverage.out.bin` per sample. +- **Phase B (merge time)**: `--runMode degNorm` reads N such files and runs the + NMF-OA pipeline, emitting DI scores and adjusted counts. + +Phase B is a pure-CPU pass over coverage matrices; it never touches the genome +index, the SA, or any read file. + +## 3. Phase A: coverage capture during alignment + +### 3.1 Parameter + +`--quantMode GeneCoverage`, composable with the existing values exactly like STAR +composes `GeneCounts TranscriptomeSAM` (`quant_mode_in: Vec`). Requires a +GTF (`--sjdbGTFfile`), same as `GeneCounts`. Enabling it implies building the +`QuantContext`. When absent, zero cost: no allocation, no branch in the hot path +beyond the existing `Option` check. + +### 3.2 Data structure + +New module `src/quant/coverage.rs`: + +```rust +pub struct GeneCoverage { + /// Prefix sums of merged-exon lengths, len = n_genes + 1 (transcript space). + offsets: Vec, + /// Flat per-base coverage, len = offsets[n_genes]. + cov: Vec, + /// Per-gene raw unique read/fragment count (same rule as GeneCounts col 1). + counts: Vec, + /// Total reads/fragments counted into any gene (library size). + n_counted: AtomicU64, +} +``` + +`offsets` is built from `GeneAnnotation.gene_exons`, which already holds merged, +sorted, absolute-coordinate exon intervals per gene. Transcript-space position of +an absolute coordinate `x` inside gene `g` is +`sum(len(e) for e in exons[..k]) + (x - exons[k].0)` where `k` is the exon +containing `x`, found by `partition_point` (same technique as the existing +`block_is_exonic`). + +Memory: human GENCODE, ~20k genes, ~70 Mb merged exonic bases → ~280 MB for `cov`. +Documented in the flag's help text and in the docs page. + +### 3.3 Read-to-coverage rule + +Reuse the `GeneCounts` unique-hit rule verbatim so counts and coverage always +agree: + +- alignment count != 1 → skipped (multimapper; DegNorm's default + `--non-unique-alignments` off); +- overlapping genes != 1 → skipped (ambiguous); +- otherwise, for each aligned block of the transcript, intersect with the gene's + merged exons and increment `cov` over the intersection. + +Paired-end: union of both mates' blocks, so an overlapping mate pair contributes 1 +to each base it covers, not 2. This matches DegNorm's paired-read handling. It is +implemented by collecting both mates' blocks, sorting, merging overlaps, then +incrementing. + +Increments are `fetch_add(1, Relaxed)` on the flat array. Contention is negligible +(distinct genes across threads dominate), and correctness does not depend on +ordering. + +### 3.4 Exons shared by several genes + +DegNorm discards exons that overlap multiple genes when it builds its gene model. +rustar-aligner's `GeneAnnotation` keeps them. The unique-gene rule in 3.3 already +drops any read that hits two genes, which is the same effect at read level. +Recorded as a documented behavioural note, not a divergence in alignment output. + +### 3.5 Output file `GeneCoverage.out.bin` + +A single gzip stream (`flate2`, already a dependency), containing: + +| field | type | note | +|---|---|---| +| magic | 8 bytes | `RSDGNCOV` | +| version | u32 | 1 | +| flags | u32 | bit 0 = paired-end run | +| n_genes | u64 | | +| total_len | u64 | sum of exonic lengths | +| n_counted | u64 | library size, for depth normalisation | +| sample_id | u16 len + bytes | derived from `--outFileNamePrefix`, overridable with `--degNormSampleId` | +| gene table | n_genes x { offset u64, len u32, chr_idx u32, strand u8, count u32 } | | +| gene id block | n_genes x { u16 len, bytes } | `gene_ids[i]` | +| coverage block | total_len x u32 LE | | + +All integers little-endian. Reader and writer live in the same module and are +round-trip tested. + +## 4. Phase B: `--runMode degNorm` + +### 4.1 Parameters + +| flag | default | meaning | +|---|---|---| +| `--degNormCoverageFiles` | (required, >= 2) | Phase A files, one per sample | +| `--degNormIter` | 5 | outer DegNorm iterations (DegNorm `--iter`) | +| `--degNormNmfIter` | 100 | inner NMF-OA iterations (`--nmf-iter`) | +| `--degNormDownsampleRate` | 1 | take-every systematic sampling (`--downsample-rate`) | +| `--degNormMinimaxCoverage` | 0 | gene included only if min-over-samples of max coverage >= this (`--minimax-coverage`) | +| `--degNormSkipBaselineSelection` | off | `--skip-baseline-selection` | +| `--degNormBins` | 20 | baseline-selection bins | +| `--degNormMinHighCoverage` | 50 | minimum high-coverage positions for baseline selection | +| `--runRNGseed` | existing | reused for downsampling's random offset | +| `--outFileNamePrefix` | existing | output location | + +Names mirror DegNorm's own CLI so the mapping is obvious; the `degNorm` prefix +keeps them clearly non-STAR. + +Validation: at least two files; identical gene id vectors and identical exonic +lengths across files (otherwise error naming the first mismatching gene); files +readable and version-compatible. + +### 4.2 Algorithm (ported from DegNorm `degnorm/nmf.py`) + +Notation follows the source: `x` = raw count matrix (n_genes x p), `F_i` = gene +i's coverage matrix (p x L_i), `rho` = DI matrix (n_genes x p). + +**Rank-one approximation.** DegNorm calls `scipy.sparse.linalg.svds(x, k=1)` and +uses `K = u*s` (p x 1), `E = v` (1 x L). Because `p` is tiny (samples), we compute +it by power iteration on the Gram matrix `G = F F^T` (p x p, built in O(p^2 L)): +iterate `u <- G u / ||G u||` to convergence (tol 1e-10, cap 1000 iters), then +`s = ||F^T u||`, `E = (F^T u / s)^T`, `K = u * s`. Sign is fixed so the dominant +component is non-negative. This is numerically the same leading singular triplet +`svds` returns; unit tests assert agreement on planted rank-one matrices. + +**NMF-OA (`nmf`).** With `c = 1/sqrt(nmf_iter)` and `lambda = 0`: + +``` +K, E = rank1(F); est = K E +repeat nmf_iter times: + res = est - F + lambda = max(lambda - c*res, 0) + K, E = rank1(F + lambda) + est = K E +``` + +**Initialisation (`ratio_svd` path).** Per gene, `est = max(rank1(F), F)` +elementwise; `rho = 1 - rowsum(F)/(rowsum(est) + 1)`. Then +`low_di = rowmax(rho) < 0.1`; `count_sums = colsum(x[low_di])` (all genes if none +qualify); `norm_factors = count_sums / median(count_sums)`; +`x_weighted = x / norm_factors`; `scale_factors = norm_factors`. + +**Baseline selection (per gene, per outer iteration).** Faithful port, including +its quirks: + +1. high-coverage columns: `colmax(F) > 0.1 * max(F)`; intersect with the + systematic downsample when `downsample_rate > 1`. +2. bail out with `rho = 0` if fewer than `min_high_coverage` such columns, or if + any sample has zero coverage over them. +3. `K, E = nmf(F_bin)`; `rho = 1 - rowsum(F_bin)/(rowsum(K E) + 1)` (the `+1` is + DegNorm's, kept). +4. bail out if `median(1 - rho) > 1`. +5. run the bin-dropping loop only if + `n_hi_cov >= max(2, ceil(200/downsample_rate))` and `min(rho) <= 0.2` and + baseline selection is not skipped. +6. loop while `max(rho) > 0.1`: per column, `res = max_over_samples(((KE - F)/(F+1))^2)`; + per bin, mean of `res`; drop the argmax bin; shift remaining bins; refit NMF; + over-approximate `KE = max(KE, F)`; recompute `rho`; stop when bins <= + `ceil(0.2*bins)` or the surviving length drops below the minimum. +7. on convergence (`max(rho) < 0.2`): clamp tiny `K`, re-derive the envelope on + the *full* high-coverage matrix as `E = colwise-max(F_start^T / K)`, recompute + `rho` from it; if `max(rho) > 0.9`, revert to the pre-baseline fit. +8. otherwise revert to the pre-baseline fit with over-approximation. + +**Outer loop** (`degnorm_iter` times): + +``` +F_adj_i = F_i / scale_factors (row-wise) +rho = baseline_selection(F_adj) (clamped to [0, 0.9]) +genes that never ran baseline selection get the sample-average DI + 1 - colsum(x_weighted)/colsum(x_adj) +x_adj = x_weighted / (1 - rho) +norm_factors = colsum(x_adj) / median(colsum(x_adj)) +x_weighted = x_weighted / norm_factors +scale_factors *= norm_factors +``` + +Parallelised over genes with `rayon` (DegNorm uses joblib threading over gene +chunks); the per-gene fit is independent, so the result is deterministic +regardless of thread count, except for the downsampling draw, which is derived +from `--runRNGseed` and the gene index rather than a shared RNG. + +### 4.3 Outputs + +Into `DegNorm.out/`: + +- `DegradationIndex.tab` — genes x samples, tab-separated, header row of sample ids. +- `AdjustedCounts.tab` — `x_adj`, same shape and header. +- `RawCounts.tab` — `x`, for provenance. +- `ScaleFactors.tab` — final per-sample sequencing-depth scale factors. +- `Summary.txt` — parameters, gene counts (total, included, baseline-selected per + iteration), per-sample median DI. + +Estimated coverage matrices (DegNorm's `.pkl` files, used only for plotting) are +out of scope. + +## 5. Module layout + +``` +src/quant/coverage.rs GeneCoverage accumulator + binary file writer/reader +src/degnorm/mod.rs run_degnorm(): load, validate, drive, write outputs +src/degnorm/nmf.rs rank1 / nmf_oa / ratio_svd / baseline_selection (pure math) +src/degnorm/io.rs output tables +``` + +`src/degnorm/nmf.rs` takes and returns plain `&[f64]` matrices with explicit +shapes and has no I/O, so it is unit-testable in isolation. `lib.rs` gains one +dispatch arm for `RunMode::DegNorm` and one call site for coverage output next to +the existing `ReadsPerGene.out.tab` write. + +## 6. Testing + +Unit: + +- absolute-coordinate to transcript-coordinate mapping across multi-exon genes, + including block boundaries and blocks spanning an intron; +- paired-end mate overlap counted once; +- `GeneCoverage.out.bin` round trip (write, read, compare all fields); +- `rank1` recovers a planted rank-one matrix (agreement with a reference singular + triplet computed by hand for a small case); +- `nmf_oa` produces an over-approximation and converges; +- DI is ~0 for an undegraded synthetic gene, and recovers the expected ordering + for a planted 3'-biased truncation in one of three samples; +- gene-set mismatch across coverage files is a clean error, not a panic. + +Integration (`tests/degnorm.rs`, following `tests/alignment_features.rs`): + +- synthetic genome plus GTF, two simulated samples (one with reads truncated + towards the 3' end of one gene), align both with `--quantMode GeneCounts + GeneCoverage`, run `--runMode degNorm`, assert the output files exist, the DI of + the degraded gene in the degraded sample exceeds that of the control, and + adjusted counts move in the expected direction. + +Optional validation script `scripts/compare_degnorm.py`: runs the Python DegNorm +on the same BAMs and correlates DI matrices. Not part of `cargo test`. + +## 7. Divergence and documentation + +This is a rustar-aligner extension, not STAR behaviour. It gets: + +- an entry in `DIVERGENCE.md` under section 4 (implementation divergences with no + intended alignment-output difference), stating that `GeneCoverage` and + `degNorm` are additions with no STAR counterpart, are off by default, and do not + affect alignment output; +- a docs page under `docs/` describing the two-phase workflow; +- `README.md` feature list and `CHANGELOG.md` entries; +- `ROADMAP.md` phase entry. + +Alignment output is bit-identical whether or not `GeneCoverage` is enabled; the +integration suite asserts this on the synthetic fixture. + +## 8. Out of scope + +- Coverage-curve plots and `.pkl` estimated coverage matrices. +- MPI / multi-node execution (DegNorm's `degnorm_mpi`). +- Warm-start directories. +- Reading third-party BAMs as Phase B input; coverage must come from Phase A. diff --git a/src/degnorm/baseline.rs b/src/degnorm/baseline.rs new file mode 100644 index 0000000..0c69809 --- /dev/null +++ b/src/degnorm/baseline.rs @@ -0,0 +1,318 @@ +//! Baseline selection: DegNorm's per-gene search for a transcript region where +//! degradation is minimal, so the coverage envelope (and therefore the DI +//! scores) can be estimated from undegraded positions. +//! +//! Port of `GeneNMFOA.baseline_selection` in DegNorm's `degnorm/nmf.py`, quirks +//! included: the `+ 1` in the DI denominator, the 0.1 / 0.2 / 0.9 thresholds, +//! and the bin-dropping loop. + +use std::collections::HashSet; + +use crate::degnorm::nmf::{Mat, nmf_oa, outer, over_approximate}; + +pub struct BaselineParams { + pub nmf_iter: usize, + pub bins: usize, + pub min_high_coverage: usize, + pub downsample_rate: usize, + pub skip: bool, +} + +pub struct BaselineResult { + /// DI score per sample for this gene. + pub rho: Vec, + /// Whether the bin-dropping search actually ran. + pub ran_baseline: bool, +} + +/// Columns whose sample-wise maximum exceeds 10% of the matrix maximum. +fn high_coverage_idx(f: &Mat) -> Vec { + let thresh = 0.1 * f.max(); + let cm = f.col_max(); + (0..f.l).filter(|&j| cm[j] > thresh).collect() +} + +/// Deterministic systematic sample. DegNorm draws the start offset at random; +/// here it is derived from the caller's seed and the gene index so results do +/// not depend on thread scheduling. +fn systematic_sample(n: usize, take_every: usize, seed_offset: u64) -> Vec { + if take_every <= 1 || take_every >= n { + return (0..n).collect(); + } + let start = (seed_offset % take_every as u64) as usize; + (start..n).step_by(take_every).collect() +} + +/// Split `0..n` into consecutive chunks of `ceil(n / bins)` (DegNorm's +/// `split_into_chunks`, which can yield fewer than `bins` chunks). +fn split_into_chunks(n: usize, bins: usize) -> Vec> { + let csize = n.div_ceil(bins.max(1)).max(1); + let mut out = Vec::new(); + let mut i = 0; + while i < n { + out.push((i..(i + csize).min(n)).collect()); + i += csize; + } + out +} + +/// DegNorm's DI formula: `1 - rowsum(F) / (rowsum(estimate) + 1)`. +fn di_scores(f: &Mat, est: &Mat) -> Vec { + let fs = f.row_sums(); + let es = est.row_sums(); + (0..f.p).map(|i| 1.0 - fs[i] / (es[i] + 1.0)).collect() +} + +fn vmax(v: &[f64]) -> f64 { + v.iter().copied().fold(f64::NEG_INFINITY, f64::max) +} + +fn vmin(v: &[f64]) -> f64 { + v.iter().copied().fold(f64::INFINITY, f64::min) +} + +fn median(v: &[f64]) -> f64 { + if v.is_empty() { + return 0.0; + } + let mut s = v.to_vec(); + s.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let n = s.len(); + if n % 2 == 1 { + s[n / 2] + } else { + f64::midpoint(s[n / 2 - 1], s[n / 2]) + } +} + +fn clamp_rho(rho: &mut [f64]) { + for r in rho.iter_mut() { + *r = r.clamp(0.0, 0.9); + } +} + +/// Fit gene `f` (shape `p x L`, already depth-scaled) and return its DI scores. +pub fn baseline_selection(f: &Mat, params: &BaselineParams, seed_offset: u64) -> BaselineResult { + let zero = || BaselineResult { + rho: vec![0.0; f.p], + ran_baseline: false, + }; + + let mut hi = high_coverage_idx(f); + if params.downsample_rate > 1 { + let sampled: HashSet = systematic_sample(f.l, params.downsample_rate, seed_offset) + .into_iter() + .collect(); + hi.retain(|j| sampled.contains(j)); + } + if hi.len() < params.min_high_coverage.max(2) { + return zero(); + } + + let f_start = f.select_cols(&hi); + if f_start.row_sums().iter().any(|&s| s <= 0.0) { + return zero(); + } + + let (k_start, e_start) = nmf_oa(&f_start, params.nmf_iter); + let mut f_bin = f_start.clone(); + let (mut k, _e) = (k_start.clone(), e_start.clone()); + let mut ke = outer(&k_start, &e_start); + let mut rho = di_scores(&f_bin, &ke); + + // Exclude extreme cases where the fit did not converge (upstream check). + let one_minus: Vec = rho.iter().map(|r| 1.0 - r).collect(); + if median(&one_minus) > 1.0 { + return zero(); + } + + let min_gene_len = (200.0 / params.downsample_rate as f64).ceil().max(2.0) as usize; + let min_bins = (params.bins as f64 * 0.2).ceil() as usize; + let mut ran_baseline = false; + + let can_run = hi.len() >= min_gene_len && vmin(&rho) <= 0.2 && !params.skip; + + if can_run { + let mut bin_segs = split_into_chunks(f_bin.l, params.bins); + + while vmax(&rho) > 0.1 { + ran_baseline = true; + + // Per-column worst squared relative residual, averaged per bin. + let res: Vec = (0..f_bin.l) + .map(|j| { + (0..f_bin.p) + .map(|i| { + let d = (ke.get(i, j) - f_bin.get(i, j)) / (f_bin.get(i, j) + 1.0); + d * d + }) + .fold(f64::NEG_INFINITY, f64::max) + }) + .collect(); + let ss_r: Vec = bin_segs + .iter() + .map(|b| b.iter().map(|&j| res[j]).sum::() / b.len() as f64) + .collect(); + if vmax(&ss_r) <= 0.0 { + break; + } + + let drop = ss_r + .iter() + .enumerate() + .max_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal)) + .map_or(0, |(i, _)| i); + + // Drop that bin's columns, then renumber the surviving bins so they + // index the shrunken matrix (DegNorm's `shift_bins`). + let dropped: HashSet = bin_segs[drop].iter().copied().collect(); + let keep: Vec = (0..f_bin.l).filter(|j| !dropped.contains(j)).collect(); + f_bin = f_bin.select_cols(&keep); + bin_segs.remove(drop); + let sizes: Vec = bin_segs.iter().map(Vec::len).collect(); + bin_segs = { + let mut out = Vec::with_capacity(sizes.len()); + let mut i = 0; + for s in sizes { + out.push((i..i + s).collect::>()); + i += s; + } + out + }; + + if f_bin.l == 0 || bin_segs.is_empty() { + break; + } + + let (k2, e2) = nmf_oa(&f_bin, params.nmf_iter); + k = k2; + ke = outer(&k, &e2); + if ke.row_sums().iter().any(|&s| s <= 0.0) { + break; + } + over_approximate(&mut ke, &f_bin); + rho = di_scores(&f_bin, &ke); + + if bin_segs.len() <= min_bins || f_bin.l < min_gene_len { + break; + } + } + + if vmax(&rho) < 0.2 { + // A baseline region was found: reuse its per-sample abundances to + // re-derive the envelope over the whole high-coverage transcript. + let mut kk: Vec = k.iter().map(|v| v.abs()).collect(); + let min_pos = kk + .iter() + .copied() + .filter(|&v| v >= 1e-5) + .fold(f64::INFINITY, f64::min); + let floor = if min_pos.is_finite() { min_pos } else { 1e-5 }; + for v in &mut kk { + if *v < 1e-5 { + *v = floor; + } + } + let ee: Vec = (0..f_start.l) + .map(|j| { + (0..f_start.p) + .map(|i| f_start.get(i, j) / kk[i]) + .fold(f64::NEG_INFINITY, f64::max) + }) + .collect(); + rho = di_scores(&f_start, &outer(&kk, &ee)); + + // Long, shallow genes can produce implausibly high DI this way; + // upstream falls back to the plain fit. + if vmax(&rho) > 0.9 { + let mut est = outer(&k_start, &e_start); + over_approximate(&mut est, &f_start); + rho = di_scores(&f_start, &est); + } + } else { + let mut est = outer(&k_start, &e_start); + over_approximate(&mut est, &f_start); + rho = di_scores(&f_start, &est); + } + } else { + // No baseline search: the DI comes from the plain over-approximation. + let mut est = outer(&k_start, &e_start); + over_approximate(&mut est, &f_start); + rho = di_scores(&f_start, &est); + } + + clamp_rho(&mut rho); + BaselineResult { rho, ran_baseline } +} + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +#[cfg(test)] +mod tests { + use super::*; + + fn params() -> BaselineParams { + BaselineParams { + nmf_iter: 100, + bins: 20, + min_high_coverage: 50, + downsample_rate: 1, + skip: false, + } + } + + #[test] + fn undegraded_gene_has_near_zero_di() { + let l = 400; + let mut f = Mat::new(2, l); + for j in 0..l { + f.set(0, j, 20.0); + f.set(1, j, 20.0); + } + let out = baseline_selection(&f, ¶ms(), 0); + assert!( + out.rho.iter().all(|&r| r.abs() < 0.05), + "rho = {:?}", + out.rho + ); + } + + #[test] + fn degraded_sample_gets_higher_di_than_control() { + let l = 400; + let mut f = Mat::new(2, l); + for j in 0..l { + f.set(0, j, 20.0); + f.set(1, j, if j < 160 { 20.0 } else { 2.0 }); + } + let out = baseline_selection(&f, ¶ms(), 0); + assert!( + out.rho[1] > out.rho[0] + 0.2, + "degraded DI {} should exceed control DI {}", + out.rho[1], + out.rho[0] + ); + assert!(out.rho.iter().all(|&r| (0.0..=0.9).contains(&r))); + } + + #[test] + fn short_gene_returns_zero_di() { + let mut f = Mat::new(2, 10); + for j in 0..10 { + f.set(0, j, 5.0); + f.set(1, j, 5.0); + } + let out = baseline_selection(&f, ¶ms(), 0); + assert_eq!(out.rho, vec![0.0, 0.0]); + assert!(!out.ran_baseline); + } + + #[test] + fn split_into_chunks_covers_every_index() { + let chunks = split_into_chunks(97, 20); + let flat: Vec = chunks.iter().flatten().copied().collect(); + assert_eq!(flat, (0..97).collect::>()); + } +} diff --git a/src/degnorm/mod.rs b/src/degnorm/mod.rs new file mode 100644 index 0000000..8782b1d --- /dev/null +++ b/src/degnorm/mod.rs @@ -0,0 +1,56 @@ +//! DegNorm-style transcript degradation normalization. +//! +//! **Not a STAR feature** (see `DIVERGENCE.md`). Two phases: +//! +//! 1. `--quantMode GeneCoverage` during alignment writes per-gene, per-exonic +//! base coverage plus raw unique counts to `GeneCoverage.out.bin` +//! ([`crate::quant::coverage`]). +//! 2. `--runMode degNorm` merges several samples' coverage files and fits +//! DegNorm's rank-one NMF over-approximation, producing Degradation Index +//! (DI) scores and degradation-adjusted read counts. +//! +//! The DI is defined *across* samples: the model separates a shared coverage +//! envelope from per-sample abundances, so at least two samples are required +//! and phase 2 cannot run inside a single alignment. +//! +//! Reference: Xiong et al., "Normalization of generalized transcript +//! degradation improves accuracy in RNA-seq analysis", *Genome Biology* 2019; +//! implementation ported from . + +pub mod baseline; +pub mod nmf; +pub mod run; + +use crate::params::Parameters; +use crate::quant::coverage::CoverageFile; + +/// `--runMode degNorm` entry point. +pub fn run_mode(params: &Parameters) -> anyhow::Result<()> { + let mut files = Vec::with_capacity(params.deg_norm_coverage_files.len()); + for path in ¶ms.deg_norm_coverage_files { + log::info!("degNorm: loading {}", path.display()); + files.push(CoverageFile::read(path)?); + } + + let cfg = run::DegNormConfig { + iter: params.deg_norm_iter, + nmf_iter: params.deg_norm_nmf_iter, + downsample_rate: params.deg_norm_downsample_rate.max(1), + minimax_coverage: params.deg_norm_minimax_coverage, + skip_baseline: params.deg_norm_skip_baseline_selection, + bins: params.deg_norm_bins, + min_high_coverage: params.deg_norm_min_high_coverage, + seed: params.run_rng_seed, + }; + + let out = run::run_degnorm(&files, &cfg)?; + let dir = params.output_path("DegNorm.out"); + run::write_outputs(&out, &dir, &cfg)?; + log::info!( + "degNorm: {} genes x {} samples written to {}", + out.gene_ids.len(), + out.sample_ids.len(), + dir.display() + ); + Ok(()) +} diff --git a/src/degnorm/nmf.rs b/src/degnorm/nmf.rs new file mode 100644 index 0000000..6608feb --- /dev/null +++ b/src/degnorm/nmf.rs @@ -0,0 +1,287 @@ +//! Rank-one NMF over-approximation, ported from DegNorm's `degnorm/nmf.py` +//! (NUStatBioinfo/DegNorm, GPL-free MIT-compatible research code; the algorithm +//! is described in Xiong et al., *Genome Biology* 2019). +//! +//! Pure math: no I/O, no globals, no randomness. +//! +//! DegNorm calls `scipy.sparse.linalg.svds(x, k=1)` for the leading singular +//! triplet. The number of rows here is the number of RNA-seq libraries (small), +//! so the triplet is obtained by power iteration on the `p x p` Gram matrix +//! `x x^T`, which is cheaper and deterministic. + +/// A dense `p x l` matrix in row-major order. +#[derive(Clone, Debug)] +pub struct Mat { + pub p: usize, + pub l: usize, + pub data: Vec, +} + +impl Mat { + pub fn new(p: usize, l: usize) -> Self { + Mat { + p, + l, + data: vec![0.0; p * l], + } + } + + #[inline] + pub fn get(&self, i: usize, j: usize) -> f64 { + self.data[i * self.l + j] + } + + #[inline] + pub fn set(&mut self, i: usize, j: usize, v: f64) { + self.data[i * self.l + j] = v; + } + + pub fn row_sums(&self) -> Vec { + (0..self.p) + .map(|i| self.data[i * self.l..(i + 1) * self.l].iter().sum()) + .collect() + } + + /// Sample-wise maximum per transcript position. + pub fn col_max(&self) -> Vec { + (0..self.l) + .map(|j| { + (0..self.p) + .map(|i| self.get(i, j)) + .fold(f64::NEG_INFINITY, f64::max) + }) + .collect() + } + + pub fn max(&self) -> f64 { + self.data.iter().copied().fold(f64::NEG_INFINITY, f64::max) + } + + /// Keep only the columns listed in `idx`, in the order given. + #[must_use] + pub fn select_cols(&self, idx: &[usize]) -> Mat { + let mut out = Mat::new(self.p, idx.len()); + for i in 0..self.p { + for (jj, &j) in idx.iter().enumerate() { + out.set(i, jj, self.get(i, j)); + } + } + out + } +} + +/// Leading singular triplet as `(k, e)` with `k[i] * e[j] ~= x[i][j]`, signs +/// fixed so both factors are non-negative for a non-negative input. +pub fn rank_one(mat: &Mat) -> (Vec, Vec) { + let n_samples = mat.p; + if n_samples == 0 || mat.l == 0 { + return (vec![0.0; n_samples], vec![0.0; mat.l]); + } + + // Gram matrix `gram = mat mat^T` (n_samples x n_samples). + let mut gram = vec![0.0f64; n_samples * n_samples]; + for row_a in 0..n_samples { + for row_b in row_a..n_samples { + let mut acc = 0.0; + for j in 0..mat.l { + acc += mat.get(row_a, j) * mat.get(row_b, j); + } + gram[row_a * n_samples + row_b] = acc; + gram[row_b * n_samples + row_a] = acc; + } + } + + // Power iteration for the dominant eigenvector of the Gram matrix. + let mut u_vec = vec![1.0 / (n_samples as f64).sqrt(); n_samples]; + for _ in 0..1000 { + let mut next = vec![0.0f64; n_samples]; + for (row, slot) in next.iter_mut().enumerate() { + *slot = (0..n_samples) + .map(|col| gram[row * n_samples + col] * u_vec[col]) + .sum(); + } + let norm = next.iter().map(|v| v * v).sum::().sqrt(); + if norm <= 0.0 { + return (vec![0.0; n_samples], vec![0.0; mat.l]); + } + for v in &mut next { + *v /= norm; + } + let delta: f64 = next.iter().zip(&u_vec).map(|(a, b)| (a - b).abs()).sum(); + u_vec = next; + if delta < 1e-12 { + break; + } + } + + // Point the dominant direction the positive way. + if u_vec.iter().sum::() < 0.0 { + for v in &mut u_vec { + *v = -*v; + } + } + + // envelope_raw = mat^T u; sigma = ||envelope_raw||; abundance = u * sigma. + let mut envelope: Vec = (0..mat.l) + .map(|j| (0..n_samples).map(|i| mat.get(i, j) * u_vec[i]).sum()) + .collect(); + let sigma = envelope.iter().map(|v| v * v).sum::().sqrt(); + if sigma <= 0.0 { + return (vec![0.0; n_samples], vec![0.0; mat.l]); + } + for v in &mut envelope { + *v /= sigma; + } + let abundance: Vec = u_vec.iter().map(|v| v * sigma).collect(); + (abundance, envelope) +} + +/// Outer product of the rank-one factors. +pub fn outer(abundance: &[f64], envelope: &[f64]) -> Mat { + let mut m = Mat::new(abundance.len(), envelope.len()); + for (i, &a) in abundance.iter().enumerate() { + for (j, &e) in envelope.iter().enumerate() { + m.set(i, j, a * e); + } + } + m +} + +/// NMF over-approximation: faithful port of `GeneNMFOA.nmf` — dual ascent on a +/// non-negative multiplier `lambda` with step `1 / sqrt(iters)`, re-fitting the +/// rank-one factors of `x + lambda` each round. +/// +/// As upstream, the returned factors are *not* clamped to dominate `x`; callers +/// that need a strict over-approximation apply [`over_approximate`] (DegNorm +/// does the same, its in-function clamp is commented out). +pub fn nmf_oa(mat: &Mat, iters: usize) -> (Vec, Vec) { + let (mut abundance, mut envelope) = rank_one(mat); + let mut est = outer(&abundance, &envelope); + let mut lambda = Mat::new(mat.p, mat.l); + let step = 1.0 / (iters.max(1) as f64).sqrt(); + + let mut shifted = Mat::new(mat.p, mat.l); + for _ in 0..iters { + for i in 0..mat.p { + for j in 0..mat.l { + let residual = est.get(i, j) - mat.get(i, j); + let dual = (lambda.get(i, j) - step * residual).max(0.0); + lambda.set(i, j, dual); + shifted.set(i, j, mat.get(i, j) + dual); + } + } + let (a2, e2) = rank_one(&shifted); + abundance = a2; + envelope = e2; + est = outer(&abundance, &envelope); + } + (abundance, envelope) +} + +/// Raise `est` elementwise to at least `f` (DegNorm's over-approximation +/// quality control). +pub fn over_approximate(est: &mut Mat, f: &Mat) { + for i in 0..f.p { + for j in 0..f.l { + if est.get(i, j) < f.get(i, j) { + est.set(i, j, f.get(i, j)); + } + } + } +} + +/// One-shot rank-one estimate raised to at least the input: DegNorm's +/// `ratio_svd`, used only to initialise the depth scale factors. +pub fn ratio_svd(mat: &Mat) -> Mat { + let (abundance, envelope) = rank_one(mat); + let mut est = outer(&abundance, &envelope); + over_approximate(&mut est, mat); + est +} + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +#[cfg(test)] +mod tests { + use super::*; + + fn mat_from(rows: &[&[f64]]) -> Mat { + let p = rows.len(); + let l = rows[0].len(); + let mut m = Mat::new(p, l); + for (i, r) in rows.iter().enumerate() { + for (j, &v) in r.iter().enumerate() { + m.set(i, j, v); + } + } + m + } + + #[test] + fn rank_one_recovers_a_planted_rank_one_matrix() { + // x = k e^T with k = [1, 2, 3], e = [4, 5]. + let x = mat_from(&[&[4.0, 5.0], &[8.0, 10.0], &[12.0, 15.0]]); + let (k, e) = rank_one(&x); + for (i, &ki) in k.iter().enumerate() { + for (j, &ej) in e.iter().enumerate() { + assert!((ki * ej - x.get(i, j)).abs() < 1e-8); + } + } + assert!(k.iter().all(|&v| v > 0.0)); + assert!(e.iter().all(|&v| v > 0.0)); + } + + #[test] + fn rank_one_of_a_zero_matrix_is_zero() { + let x = Mat::new(2, 3); + let (k, e) = rank_one(&x); + assert!(k.iter().all(|&v| v == 0.0)); + assert!(e.iter().all(|&v| v == 0.0)); + } + + #[test] + fn nmf_oa_moves_the_estimate_above_the_rank_one_fit() { + // Sample 2 is degraded at the 3' end: a rank-one fit cannot cover it. + let x = mat_from(&[&[10.0, 10.0, 10.0, 10.0], &[10.0, 10.0, 2.0, 1.0]]); + let (k0, e0) = rank_one(&x); + let (k, e) = nmf_oa(&x, 100); + let plain = outer(&k0, &e0); + let lifted = outer(&k, &e); + let deficit = |m: &Mat| -> f64 { + let mut d = 0.0; + for i in 0..x.p { + for j in 0..x.l { + d += (x.get(i, j) - m.get(i, j)).max(0.0); + } + } + d + }; + assert!( + deficit(&lifted) < deficit(&plain), + "NMF-OA should reduce the under-approximation deficit: {} vs {}", + deficit(&lifted), + deficit(&plain) + ); + } + + #[test] + fn ratio_svd_dominates_input() { + let x = mat_from(&[&[5.0, 1.0], &[1.0, 5.0]]); + let est = ratio_svd(&x); + for i in 0..2 { + for j in 0..2 { + assert!(est.get(i, j) >= x.get(i, j) - 1e-9); + } + } + } + + #[test] + fn select_cols_keeps_requested_positions() { + let x = mat_from(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]]); + let sub = x.select_cols(&[0, 2]); + assert_eq!(sub.l, 2); + assert!((sub.get(1, 1) - 6.0).abs() < 1e-12); + } +} diff --git a/src/degnorm/run.rs b/src/degnorm/run.rs new file mode 100644 index 0000000..004ab7c --- /dev/null +++ b/src/degnorm/run.rs @@ -0,0 +1,447 @@ +//! DegNorm driver: load per-sample coverage, run the NMF-OA iterations, and +//! write DI scores plus degradation-adjusted counts. +//! +//! Port of `GeneNMFOA.run` in DegNorm's `degnorm/nmf.py`. + +use std::io::Write; +use std::path::Path; + +use rayon::prelude::*; + +use crate::degnorm::baseline::{BaselineParams, baseline_selection}; +use crate::degnorm::nmf::{Mat, ratio_svd}; +use crate::error::Error; +use crate::quant::coverage::CoverageFile; + +pub struct DegNormConfig { + pub iter: usize, + pub nmf_iter: usize, + pub downsample_rate: usize, + pub minimax_coverage: u32, + pub skip_baseline: bool, + pub bins: usize, + pub min_high_coverage: usize, + pub seed: u64, +} + +pub struct DegNormOutput { + pub gene_ids: Vec, + pub sample_ids: Vec, + /// DI score per included gene, per sample. + pub rho: Vec>, + pub raw_counts: Vec>, + pub adjusted_counts: Vec>, + pub scale_factors: Vec, + /// Genes sent through baseline selection, per outer iteration. + pub n_baseline_selected: Vec, +} + +fn median(v: &[f64]) -> f64 { + if v.is_empty() { + return 0.0; + } + let mut s = v.to_vec(); + s.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let n = s.len(); + if n % 2 == 1 { + s[n / 2] + } else { + f64::midpoint(s[n / 2 - 1], s[n / 2]) + } +} + +/// Run the full DegNorm pipeline over `files` (one per sample). +pub fn run_degnorm(files: &[CoverageFile], cfg: &DegNormConfig) -> Result { + let p = files.len(); + if p < 2 { + return Err(Error::Parameter( + "--runMode degNorm requires at least 2 coverage files: a degradation index is \ + defined across samples, not within one" + .to_string(), + )); + } + + // All samples must share one gene model. + for f in &files[1..] { + if f.gene_ids != files[0].gene_ids || f.gene_lens != files[0].gene_lens { + let first_diff = f + .gene_ids + .iter() + .zip(&files[0].gene_ids) + .position(|(a, b)| a != b) + .map_or_else( + || "".to_string(), + |i| f.gene_ids[i].clone(), + ); + return Err(Error::Parameter(format!( + "coverage file '{}' has a different gene model than '{}' (first difference: {}); \ + all samples must be aligned against the same GTF", + f.sample_id, files[0].sample_id, first_diff + ))); + } + } + + let sample_ids: Vec = files.iter().map(|f| f.sample_id.clone()).collect(); + + // Gene selection: every sample must reach the minimax coverage threshold. + let min_cov = cfg.minimax_coverage.max(1); + let kept: Vec = (0..files[0].n_genes()) + .filter(|&g| { + files[0].gene_lens[g] >= 2 + && files + .iter() + .all(|f| f.gene(g).iter().copied().max().unwrap_or(0) >= min_cov) + }) + .collect(); + if kept.is_empty() { + return Err(Error::Parameter( + "no gene reached --degNormMinimaxCoverage in every sample; nothing to normalize" + .to_string(), + )); + } + let n_genes = kept.len(); + + // Coverage matrices, one per kept gene, shape p x L. + let mats: Vec = kept + .iter() + .map(|&g| { + let l = files[0].gene_lens[g] as usize; + let mut m = Mat::new(p, l); + for (i, f) in files.iter().enumerate() { + for (j, &v) in f.gene(g).iter().enumerate() { + m.set(i, j, f64::from(v)); + } + } + m + }) + .collect(); + + // Raw count matrix (genes x samples). + let x: Vec> = kept + .iter() + .map(|&g| files.iter().map(|f| f64::from(f.counts[g])).collect()) + .collect(); + + // Initialisation: one-shot over-approximation gives the starting DI scores. + let mut rho: Vec> = mats + .par_iter() + .map(|m| { + let est = ratio_svd(m); + let fs = m.row_sums(); + let es = est.row_sums(); + (0..m.p) + .map(|i| (1.0 - fs[i] / (es[i] + 1.0)).clamp(0.0, 0.9)) + .collect() + }) + .collect(); + + let low_di: Vec = (0..n_genes) + .filter(|&g| rho[g].iter().copied().fold(f64::NEG_INFINITY, f64::max) < 0.1) + .collect(); + let count_sums: Vec = (0..p) + .map(|i| { + if low_di.is_empty() { + (0..n_genes).map(|g| x[g][i]).sum() + } else { + low_di.iter().map(|&g| x[g][i]).sum() + } + }) + .collect(); + let med = median(&count_sums); + let mut norm_factors: Vec = count_sums + .iter() + .map(|&c| if med > 0.0 { (c / med).max(1e-12) } else { 1.0 }) + .collect(); + let mut x_weighted: Vec> = (0..n_genes) + .map(|g| (0..p).map(|i| x[g][i] / norm_factors[i]).collect()) + .collect(); + let mut scale_factors = norm_factors.clone(); + let mut x_adj = x_weighted.clone(); + let mut n_baseline_selected = Vec::with_capacity(cfg.iter); + + let bp = BaselineParams { + nmf_iter: cfg.nmf_iter, + bins: cfg.bins, + min_high_coverage: cfg.min_high_coverage, + downsample_rate: cfg.downsample_rate, + skip: cfg.skip_baseline, + }; + + for it in 0..cfg.iter { + // Scale each sample's coverage curve by its depth factor, then fit. + let results: Vec<(Vec, bool)> = mats + .par_iter() + .enumerate() + .map(|(gi, m)| { + let mut adj = Mat::new(m.p, m.l); + for (i, factor) in scale_factors.iter().enumerate().take(m.p) { + let s = factor.max(1e-12); + for j in 0..m.l { + adj.set(i, j, m.get(i, j) / s); + } + } + let r = baseline_selection(&adj, &bp, cfg.seed.wrapping_add(gi as u64)); + (r.rho, r.ran_baseline) + }) + .collect(); + + let n_sel = results.iter().filter(|(_, b)| *b).count(); + n_baseline_selected.push(n_sel); + log::info!( + "degNorm iteration {}: {n_sel} genes through baseline selection", + it + 1 + ); + for (g, (r, _)) in results.into_iter().enumerate() { + rho[g] = r; + } + + // Genes that never went through baseline selection inherit the + // sample-average DI (DegNorm's `correct_di_scores`). + x_adj = adjust_counts(&x_weighted, &rho, n_genes, p); + let sample_avg_di: Vec = (0..p) + .map(|i| { + let w: f64 = (0..n_genes).map(|g| x_weighted[g][i]).sum(); + let a: f64 = (0..n_genes).map(|g| x_adj[g][i]).sum(); + if a > 0.0 { 1.0 - w / a } else { 0.0 } + }) + .collect(); + for r in rho.iter_mut().take(n_genes) { + if r.iter().copied().fold(f64::NEG_INFINITY, f64::max) <= 0.0 { + r.clone_from(&sample_avg_di); + } + } + + x_adj = adjust_counts(&x_weighted, &rho, n_genes, p); + + // Fold the degradation correction back into the depth factors. + let col_sums: Vec = (0..p) + .map(|i| (0..n_genes).map(|g| x_adj[g][i]).sum()) + .collect(); + let med = median(&col_sums); + norm_factors = col_sums + .iter() + .map(|&c| if med > 0.0 { (c / med).max(1e-12) } else { 1.0 }) + .collect(); + for row in x_weighted.iter_mut().take(n_genes) { + for (v, nf) in row.iter_mut().zip(&norm_factors) { + *v /= *nf; + } + } + for (s, nf) in scale_factors.iter_mut().zip(&norm_factors) { + *s *= *nf; + } + } + + Ok(DegNormOutput { + gene_ids: kept.iter().map(|&g| files[0].gene_ids[g].clone()).collect(), + sample_ids, + rho, + raw_counts: x, + adjusted_counts: x_adj, + scale_factors, + n_baseline_selected, + }) +} + +fn adjust_counts( + x_weighted: &[Vec], + rho: &[Vec], + n_genes: usize, + p: usize, +) -> Vec> { + (0..n_genes) + .map(|g| { + (0..p) + .map(|i| x_weighted[g][i] / (1.0 - rho[g][i]).max(1e-12)) + .collect() + }) + .collect() +} + +fn write_matrix( + path: &Path, + row_ids: &[String], + header: &[String], + rows: &[Vec], +) -> Result<(), Error> { + let mut f = std::fs::File::create(path).map_err(|e| Error::io(e, path))?; + let io = |e: std::io::Error| Error::io(e, path); + writeln!(f, "gene\t{}", header.join("\t")).map_err(io)?; + for (g, id) in row_ids.iter().enumerate() { + let cells: Vec = rows[g].iter().map(|v| format!("{v:.6}")).collect(); + writeln!(f, "{id}\t{}", cells.join("\t")).map_err(io)?; + } + Ok(()) +} + +/// Write `DegradationIndex.tab`, `AdjustedCounts.tab`, `RawCounts.tab`, +/// `ScaleFactors.tab`, and `Summary.txt` into `dir`. +pub fn write_outputs(out: &DegNormOutput, dir: &Path, cfg: &DegNormConfig) -> Result<(), Error> { + std::fs::create_dir_all(dir).map_err(|e| Error::io(e, dir))?; + write_matrix( + &dir.join("DegradationIndex.tab"), + &out.gene_ids, + &out.sample_ids, + &out.rho, + )?; + write_matrix( + &dir.join("AdjustedCounts.tab"), + &out.gene_ids, + &out.sample_ids, + &out.adjusted_counts, + )?; + write_matrix( + &dir.join("RawCounts.tab"), + &out.gene_ids, + &out.sample_ids, + &out.raw_counts, + )?; + + let sf_path = dir.join("ScaleFactors.tab"); + let mut f = std::fs::File::create(&sf_path).map_err(|e| Error::io(e, &sf_path))?; + { + let io = |e: std::io::Error| Error::io(e, sf_path.as_path()); + writeln!(f, "sample\tscaleFactor").map_err(io)?; + for (i, s) in out.sample_ids.iter().enumerate() { + writeln!(f, "{s}\t{:.6}", out.scale_factors[i]).map_err(io)?; + } + } + + let sum_path = dir.join("Summary.txt"); + let mut f = std::fs::File::create(&sum_path).map_err(|e| Error::io(e, &sum_path))?; + let io = |e: std::io::Error| Error::io(e, sum_path.as_path()); + writeln!(f, "samples\t{}", out.sample_ids.len()).map_err(io)?; + writeln!(f, "genesIncluded\t{}", out.gene_ids.len()).map_err(io)?; + writeln!(f, "degNormIter\t{}", cfg.iter).map_err(io)?; + writeln!(f, "degNormNmfIter\t{}", cfg.nmf_iter).map_err(io)?; + writeln!(f, "downsampleRate\t{}", cfg.downsample_rate).map_err(io)?; + writeln!(f, "minimaxCoverage\t{}", cfg.minimax_coverage).map_err(io)?; + writeln!(f, "skipBaselineSelection\t{}", cfg.skip_baseline).map_err(io)?; + for (i, n) in out.n_baseline_selected.iter().enumerate() { + writeln!(f, "baselineSelectedIter{}\t{}", i + 1, n).map_err(io)?; + } + for (i, s) in out.sample_ids.iter().enumerate() { + let col: Vec = out.rho.iter().map(|r| r[i]).collect(); + writeln!(f, "medianDI_{}\t{:.6}", s, median(&col)).map_err(io)?; + } + Ok(()) +} + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +#[cfg(test)] +mod tests { + use super::*; + + /// Build an in-memory `CoverageFile` with `n_genes` genes of length `len`, + /// coverage supplied by `f(gene, position)`. + fn synth( + sample: &str, + n_genes: usize, + len: usize, + f: impl Fn(usize, usize) -> u32, + ) -> CoverageFile { + let mut cov = Vec::new(); + let mut counts = Vec::new(); + let mut offsets = vec![0u64]; + for g in 0..n_genes { + let mut total = 0u64; + for j in 0..len { + let v = f(g, j); + total += u64::from(v); + cov.push(v); + } + counts.push((total / 100).max(1) as u32); + offsets.push(offsets[g] + len as u64); + } + CoverageFile { + sample_id: sample.to_string(), + paired: false, + n_counted: counts.iter().map(|&c| u64::from(c)).sum(), + gene_ids: (0..n_genes).map(|g| format!("G{g}")).collect(), + gene_lens: vec![len as u32; n_genes], + counts, + offsets, + cov, + } + } + + fn cfg() -> DegNormConfig { + DegNormConfig { + iter: 2, + nmf_iter: 50, + downsample_rate: 1, + minimax_coverage: 0, + skip_baseline: false, + bins: 20, + min_high_coverage: 50, + seed: 777, + } + } + + #[test] + fn degraded_sample_gets_higher_di_and_upweighted_counts() { + let len = 400; + let a = synth("ctrl", 2, len, |_, _| 20); + let b = synth( + "degraded", + 2, + len, + |g, j| { + if g == 0 && j >= 160 { 2 } else { 20 } + }, + ); + let out = run_degnorm(&[a, b], &cfg()).unwrap(); + + assert_eq!( + out.sample_ids, + vec!["ctrl".to_string(), "degraded".to_string()] + ); + assert_eq!(out.gene_ids.len(), 2); + assert!( + out.rho[0][1] > out.rho[0][0] + 0.2, + "DI of the degraded sample ({}) should exceed the control ({})", + out.rho[0][1], + out.rho[0][0] + ); + assert!(out.adjusted_counts[0][1] > out.raw_counts[0][1] * 0.5); + } + + #[test] + fn mismatched_gene_sets_are_rejected() { + let a = synth("a", 2, 400, |_, _| 10); + let mut b = synth("b", 2, 400, |_, _| 10); + b.gene_ids[1] = "OTHER".to_string(); + assert!(run_degnorm(&[a, b], &cfg()).is_err()); + } + + #[test] + fn single_sample_is_rejected() { + let a = synth("a", 2, 400, |_, _| 10); + assert!(run_degnorm(&[a], &cfg()).is_err()); + } + + #[test] + fn outputs_are_written() { + let a = synth("a", 2, 400, |_, _| 20); + let b = synth("b", 2, 400, |_, j| if j >= 200 { 3 } else { 20 }); + let out = run_degnorm(&[a, b], &cfg()).unwrap(); + let dir = tempfile::tempdir().unwrap(); + let target = dir.path().join("DegNorm.out"); + write_outputs(&out, &target, &cfg()).unwrap(); + for name in [ + "DegradationIndex.tab", + "AdjustedCounts.tab", + "RawCounts.tab", + "ScaleFactors.tab", + "Summary.txt", + ] { + assert!(target.join(name).exists(), "{name} missing"); + } + let di = std::fs::read_to_string(target.join("DegradationIndex.tab")).unwrap(); + assert!(di.starts_with("gene\ta\tb\n")); + assert_eq!(di.lines().count(), 3); + } +} diff --git a/src/lib.rs b/src/lib.rs index 5086fcb..adf5d18 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -29,6 +29,7 @@ pub mod bam_dedup; pub mod chimeric; pub mod clip; pub mod cpu; +pub mod degnorm; pub mod genome; pub mod index; pub mod io; @@ -82,6 +83,7 @@ pub fn run(params: &Parameters) -> anyhow::Result<()> { RunMode::InputAlignmentsFromBAM => bam_dedup::run(params), RunMode::LiftOver => liftover::run(params), RunMode::SoloCellFiltering => crate::solo::count::run_cell_filtering(params), + RunMode::DegNorm => degnorm::run_mode(params), } } @@ -300,10 +302,10 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { // Build gene-count context if --quantMode GeneCounts was requested. // GTF requirement is already validated in params.validate(). let quant_ctx: Option> = - if params.quant_gene_counts() { + if params.quant_gene_counts() || params.quant_gene_coverage() { let gtf_path = params.sjdb_gtf_file.as_ref().unwrap(); info!( - "quantMode GeneCounts: building gene annotation from {}", + "quantMode GeneCounts/GeneCoverage: building gene annotation from {}", gtf_path.display() ); let ctx = crate::quant::QuantContext::build( @@ -312,6 +314,7 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { ¶ms.sjdb_gtf_feature_exon, ¶ms.sjdb_gtf_chr_prefix, ¶ms.sjdb_gtf_tag_exon_parent_gene, + params.quant_gene_coverage(), )?; Some(std::sync::Arc::new(ctx)) } else { @@ -425,9 +428,19 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { // Write ReadsPerGene.out.tab if quantMode GeneCounts was requested. if let Some(ref ctx) = quant_ctx { - let quant_path = params.output_path("ReadsPerGene.out.tab"); - ctx.counts.write_output(&quant_path, &ctx.gene_ann)?; - info!("Wrote {}", quant_path.display()); + if params.quant_gene_counts() { + let quant_path = params.output_path("ReadsPerGene.out.tab"); + ctx.counts.write_output(&quant_path, &ctx.gene_ann)?; + info!("Wrote {}", quant_path.display()); + } + // quantMode GeneCoverage (DegNorm phase 1, not a STAR mode). + if let Some(ref cov) = ctx.coverage { + let cov_path = params.output_path("GeneCoverage.out.bin"); + let sample_id = params.deg_norm_sample_id_or_default(); + let is_paired = params.read_files_in.len() == 2 && !params.solo_enabled(); + cov.write_file(&cov_path, &ctx.gene_ann, &sample_id, is_paired)?; + info!("Wrote {}", cov_path.display()); + } } info!("Alignment complete!"); @@ -1881,6 +1894,9 @@ fn align_reads_single_end( if let Some(ref q) = quant { q.counts .count_se_read(&transcripts, n_for_mapq, &q.gene_ann); + if let Some(ref cov) = q.coverage { + cov.count_se_read(&transcripts, &q.gene_ann); + } } // Record junction statistics (per-read dedup, fix A) @@ -3262,6 +3278,9 @@ fn align_reads_paired_end( has_half_mapped, &q.gene_ann, ); + if let Some(ref cov) = q.coverage { + cov.count_pe_read(&bm_deref, &q.gene_ann); + } } // Record junction statistics diff --git a/src/params/mod.rs b/src/params/mod.rs index 3b50731..c27c1c4 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -39,6 +39,8 @@ pub enum RunMode { LiftOver, /// Cell-call an existing raw count matrix, without aligning anything. SoloCellFiltering, + /// rustar-aligner extension: DegNorm degradation normalization merge pass. + DegNorm, } impl std::str::FromStr for RunMode { @@ -50,9 +52,10 @@ impl std::str::FromStr for RunMode { "inputAlignmentsFromBAM" => Ok(Self::InputAlignmentsFromBAM), "liftOver" => Ok(Self::LiftOver), "soloCellFiltering" => Ok(Self::SoloCellFiltering), + "degNorm" => Ok(Self::DegNorm), _ => Err(format!( "unknown runMode '{s}'; expected 'alignReads', 'genomeGenerate', \ - 'inputAlignmentsFromBAM', 'liftOver', or 'soloCellFiltering'" + 'inputAlignmentsFromBAM', 'liftOver', 'soloCellFiltering', or 'degNorm'" )), } } @@ -66,6 +69,7 @@ impl std::fmt::Display for RunMode { Self::InputAlignmentsFromBAM => write!(f, "inputAlignmentsFromBAM"), Self::LiftOver => write!(f, "liftOver"), Self::SoloCellFiltering => write!(f, "soloCellFiltering"), + Self::DegNorm => write!(f, "degNorm"), } } } @@ -1124,6 +1128,50 @@ pub struct Parameters { )] pub quant_transcriptome_sam_output: crate::quant::transcriptome::QuantTranscriptomeSAMoutput, + // ── DegNorm (rustar-aligner extension, not STAR) ──────────────────── + /// Sample name recorded in `GeneCoverage.out.bin` + /// (`--quantMode GeneCoverage`). Defaults to the basename of + /// `--outFileNamePrefix`. + #[arg(long = "degNormSampleId", default_value = "")] + pub deg_norm_sample_id: String, + + /// `GeneCoverage.out.bin` files, one per sample, for `--runMode degNorm`. + /// At least two are required: a degradation index is defined across + /// samples, not within one. + #[arg(long = "degNormCoverageFiles", num_args = 0.., value_name = "FILE")] + pub deg_norm_coverage_files: Vec, + + /// Outer DegNorm iterations (DegNorm's `--iter`). + #[arg(long = "degNormIter", default_value_t = 5)] + pub deg_norm_iter: usize, + + /// NMF-OA iterations per gene per outer iteration (DegNorm's `--nmf-iter`). + #[arg(long = "degNormNmfIter", default_value_t = 100)] + pub deg_norm_nmf_iter: usize, + + /// Systematic "take every" downsample rate over transcript positions + /// (DegNorm's `--downsample-rate`); 1 = no downsampling. + #[arg(long = "degNormDownsampleRate", default_value_t = 1)] + pub deg_norm_downsample_rate: usize, + + /// A gene enters the DegNorm fit only if its maximum coverage reaches this + /// value in every sample (DegNorm's `--minimax-coverage`). + #[arg(long = "degNormMinimaxCoverage", default_value_t = 0)] + pub deg_norm_minimax_coverage: u32, + + /// Skip baseline selection: faster, less accurate DI scores. + #[arg(long = "degNormSkipBaselineSelection", default_value_t = false)] + pub deg_norm_skip_baseline_selection: bool, + + /// Number of bins used by baseline selection. + #[arg(long = "degNormBins", default_value_t = 20)] + pub deg_norm_bins: usize, + + /// Minimum number of high-coverage transcript positions required before a + /// gene is sent through baseline selection. + #[arg(long = "degNormMinHighCoverage", default_value_t = 50)] + pub deg_norm_min_high_coverage: usize, + // ── Two-pass ──────────────────────────────────────────────────────── /// Two-pass mode: None or Basic #[arg(long = "twopassMode", default_value = "None")] @@ -1634,6 +1682,16 @@ impl Parameters { } } + // degNorm merges per-sample GeneCoverage files; a degradation index is + // defined across samples, so one file is never enough. + if params.run_mode() == RunMode::DegNorm && params.deg_norm_coverage_files.len() < 2 { + return Err(command.error( + ErrorKind::MissingRequiredArgument, + "--runMode degNorm requires at least 2 --degNormCoverageFiles (one per sample); \ + a degradation index is defined across samples, not within one", + )); + } + // --outWigType at alignReads: only `bedGraph` (stranded, full-length) is // implemented. `wiggle`, `--outWigStrand Unstranded`, and the 2nd word // (`read1_5p`/`read2`) are STAR features of `--runMode @@ -1677,6 +1735,15 @@ impl Parameters { )); } + // quantMode GeneCoverage (rustar-aligner extension, DegNorm phase 1) + // needs the same gene model. + if params.quant_gene_coverage() && params.sjdb_gtf_file.is_none() { + return Err(command.error( + ErrorKind::MissingRequiredArgument, + "--quantMode GeneCoverage requires --sjdbGTFfile", + )); + } + // Read group: `RG` in outSAMattributes without an RG line is a fatal // error (STAR: Parameters_samAttributes.cpp:206). STAR's "All" preset // does NOT include RG, so only match a literal user-supplied RG flag. @@ -2074,6 +2141,31 @@ impl Parameters { self.quant_mode.iter().any(|m| m == "GeneCounts") } + /// Returns true if `--quantMode GeneCoverage` was requested. + /// + /// Not a STAR mode: this is the rustar-aligner extension that captures + /// per-gene exonic coverage for the DegNorm pipeline (`--runMode degNorm`). + pub fn quant_gene_coverage(&self) -> bool { + self.quant_mode.iter().any(|m| m == "GeneCoverage") + } + + /// Sample id recorded in `GeneCoverage.out.bin`: `--degNormSampleId` when + /// set, otherwise the basename of `--outFileNamePrefix`, otherwise + /// `sample`. + pub fn deg_norm_sample_id_or_default(&self) -> String { + if !self.deg_norm_sample_id.is_empty() { + return self.deg_norm_sample_id.clone(); + } + let trimmed = self.out_file_name_prefix.trim_end_matches('/'); + let base = trimmed.rsplit('/').next().unwrap_or(""); + let base = base.trim_end_matches(['.', '_']); + if base.is_empty() { + "sample".to_string() + } else { + base.to_string() + } + } + /// Returns true if `--quantMode TranscriptomeSAM` was requested. pub fn quant_transcriptome_sam(&self) -> bool { self.quant_mode.iter().any(|m| m == "TranscriptomeSAM") @@ -2185,6 +2277,73 @@ mod tests { Parameters::try_parse_from(&full) } + #[test] + fn quant_mode_gene_coverage_is_recognised() { + let p = try_parse(&[ + "--readFilesIn", + "r.fq", + "--quantMode", + "GeneCoverage", + "--sjdbGTFfile", + "a.gtf", + ]) + .unwrap(); + assert!(p.quant_gene_coverage()); + assert!(!p.quant_gene_counts()); + } + + #[test] + fn quant_mode_gene_coverage_requires_gtf() { + assert!(try_parse(&["--readFilesIn", "r.fq", "--quantMode", "GeneCoverage"]).is_err()); + } + + #[test] + fn deg_norm_sample_id_falls_back_to_output_prefix() { + let p = try_parse(&[ + "--readFilesIn", + "r.fq", + "--outFileNamePrefix", + "out/sampleA_", + ]) + .unwrap(); + assert_eq!(p.deg_norm_sample_id_or_default(), "sampleA"); + let p = try_parse(&[ + "--readFilesIn", + "r.fq", + "--outFileNamePrefix", + "out/sampleA_", + "--degNormSampleId", + "explicit", + ]) + .unwrap(); + assert_eq!(p.deg_norm_sample_id_or_default(), "explicit"); + } + + #[test] + fn run_mode_degnorm_parses_and_validates() { + let p = try_parse(&[ + "--runMode", + "degNorm", + "--degNormCoverageFiles", + "a.bin", + "b.bin", + ]) + .unwrap(); + assert_eq!(p.run_mode(), RunMode::DegNorm); + assert_eq!(p.deg_norm_coverage_files.len(), 2); + assert_eq!(p.deg_norm_iter, 5); + assert_eq!(p.deg_norm_nmf_iter, 100); + assert_eq!(p.deg_norm_bins, 20); + assert_eq!(p.deg_norm_min_high_coverage, 50); + assert_eq!(p.deg_norm_downsample_rate, 1); + } + + #[test] + fn run_mode_degnorm_requires_two_coverage_files() { + assert!(try_parse(&["--runMode", "degNorm", "--degNormCoverageFiles", "a.bin"]).is_err()); + assert!(try_parse(&["--runMode", "degNorm"]).is_err()); + } + #[test] fn defaults() { let p = try_parse(&["--readFilesIn", "reads.fq"]).unwrap(); diff --git a/src/quant/coverage.rs b/src/quant/coverage.rs new file mode 100644 index 0000000..9cf3c42 --- /dev/null +++ b/src/quant/coverage.rs @@ -0,0 +1,488 @@ +//! Per-gene, per-exonic-base coverage accumulation (`--quantMode GeneCoverage`). +//! +//! **Not a STAR feature.** This is the input half of the DegNorm degradation +//! normalization pipeline (see [`crate::degnorm`]), captured during alignment so +//! that no sorted BAM has to be written, indexed, and re-read to obtain gene +//! coverage curves. +//! +//! Coverage lives in *transcript space*: a gene's merged, sorted exons are +//! concatenated, and position `j` of gene `g` is the `j`-th exonic base of that +//! gene. The DegNorm model fits a rank-one envelope over exactly these +//! coordinates. + +use std::io::{BufReader, BufWriter, Read, Write}; +use std::path::Path; +use std::sync::atomic::{AtomicU32, AtomicU64, Ordering}; + +use flate2::Compression; +use flate2::read::GzDecoder; +use flate2::write::GzEncoder; + +use crate::align::read_align::PairedAlignment; +use crate::align::transcript::Transcript; +use crate::error::Error; +use crate::quant::GeneAnnotation; + +const COV_MAGIC: &[u8; 8] = b"RSDGNCOV"; +const COV_VERSION: u32 = 1; + +/// Lock-free per-gene coverage in transcript space. +pub struct GeneCoverage { + /// Prefix sums of per-gene transcript lengths; `len == n_genes + 1`. + offsets: Vec, + /// Flat coverage array, `len == offsets[n_genes]`. + cov: Vec, + /// Per-gene raw unique read/fragment counts (the `GeneCounts` column-1 rule). + counts: Vec, + /// Total reads/fragments counted into any gene (library size). + n_counted: AtomicU64, +} + +impl GeneCoverage { + pub fn new(ann: &GeneAnnotation) -> Self { + let n = ann.n_genes(); + let mut offsets = Vec::with_capacity(n + 1); + let mut acc: u64 = 0; + offsets.push(0); + for exons in ann.gene_exons.iter().take(n) { + acc += exons.iter().map(|&(s, e)| e - s).sum::(); + offsets.push(acc); + } + GeneCoverage { + offsets, + cov: (0..acc).map(|_| AtomicU32::new(0)).collect(), + counts: (0..n).map(|_| AtomicU32::new(0)).collect(), + n_counted: AtomicU64::new(0), + } + } + + pub fn n_genes(&self) -> usize { + self.counts.len() + } + + pub fn gene_len(&self, g: usize) -> u32 { + (self.offsets[g + 1] - self.offsets[g]) as u32 + } + + pub fn total_len(&self) -> u64 { + self.offsets[self.offsets.len() - 1] + } + + pub fn n_counted(&self) -> u64 { + self.n_counted.load(Ordering::Relaxed) + } + + pub fn gene_count(&self, g: usize) -> u32 { + self.counts[g].load(Ordering::Relaxed) + } + + /// Increment gene `g`'s coverage over the genomic interval `[start, end)`, + /// restricted to the gene's merged exons. + pub fn add_block(&self, g: usize, start: u64, end: u64, ann: &GeneAnnotation) { + let mut tx_off: u64 = 0; + for &(es, ee) in &ann.gene_exons[g] { + if end > es && start < ee { + let s = start.max(es); + let e = end.min(ee); + let base = self.offsets[g] + tx_off + (s - es); + for i in 0..(e - s) { + self.cov[(base + i) as usize].fetch_add(1, Ordering::Relaxed); + } + } + tx_off += ee - es; + } + } + + /// Single-end: accumulate coverage for a uniquely mapped read assigned to + /// exactly one gene. + /// + /// Mirrors the `GeneCounts` column-1 rule so counts and coverage always + /// agree; multimappers and gene-ambiguous reads are skipped, matching + /// DegNorm's default (unique alignments only). + pub fn count_se_read(&self, transcripts: &[Transcript], ann: &GeneAnnotation) { + if transcripts.len() != 1 { + return; + } + let t = &transcripts[0]; + let mut genes = Vec::new(); + ann.overlapping_genes_into(t, &mut genes); + if genes.len() != 1 { + return; + } + let g = genes[0]; + for ex in &t.exons { + self.add_block(g, ex.genome_start, ex.genome_end, ann); + } + self.counts[g].fetch_add(1, Ordering::Relaxed); + self.n_counted.fetch_add(1, Ordering::Relaxed); + } + + /// Paired-end: the fragment is one observation. Mate blocks are merged + /// before accumulation, so an overlapping pair contributes 1 per base, which + /// is what DegNorm's paired-read coverage does. + pub fn count_pe_read(&self, both_mapped: &[&PairedAlignment], ann: &GeneAnnotation) { + if both_mapped.len() != 1 { + return; + } + let pair = both_mapped[0]; + let mut genes = Vec::new(); + let mut genes2 = Vec::new(); + ann.overlapping_genes_into(&pair.mate1_transcript, &mut genes); + ann.overlapping_genes_into(&pair.mate2_transcript, &mut genes2); + genes.extend_from_slice(&genes2); + genes.sort_unstable(); + genes.dedup(); + if genes.len() != 1 { + return; + } + let g = genes[0]; + + let mut blocks: Vec<(u64, u64)> = pair + .mate1_transcript + .exons + .iter() + .chain(pair.mate2_transcript.exons.iter()) + .map(|e| (e.genome_start, e.genome_end)) + .collect(); + blocks.sort_unstable(); + let mut merged: Vec<(u64, u64)> = Vec::with_capacity(blocks.len()); + for (s, e) in blocks { + match merged.last_mut() { + Some(last) if s <= last.1 => last.1 = last.1.max(e), + _ => merged.push((s, e)), + } + } + for (s, e) in merged { + self.add_block(g, s, e, ann); + } + self.counts[g].fetch_add(1, Ordering::Relaxed); + self.n_counted.fetch_add(1, Ordering::Relaxed); + } + + /// Read back one gene's coverage vector. + pub fn gene_slice(&self, g: usize) -> Vec { + (self.offsets[g]..self.offsets[g + 1]) + .map(|i| self.cov[i as usize].load(Ordering::Relaxed)) + .collect() + } + + /// Write `GeneCoverage.out.bin`: a gzip stream holding a fixed header, a + /// per-gene table, the gene id block, and the flat coverage array. + pub fn write_file( + &self, + path: &Path, + ann: &GeneAnnotation, + sample_id: &str, + paired: bool, + ) -> Result<(), Error> { + let file = std::fs::File::create(path).map_err(|e| Error::io(e, path))?; + let mut w = GzEncoder::new(BufWriter::new(file), Compression::new(6)); + let io = |e: std::io::Error| Error::io(e, path); + + w.write_all(COV_MAGIC).map_err(io)?; + w.write_all(&COV_VERSION.to_le_bytes()).map_err(io)?; + w.write_all(&u32::from(paired).to_le_bytes()).map_err(io)?; + w.write_all(&(self.n_genes() as u64).to_le_bytes()) + .map_err(io)?; + w.write_all(&self.total_len().to_le_bytes()).map_err(io)?; + w.write_all(&self.n_counted().to_le_bytes()).map_err(io)?; + + let sid = sample_id.as_bytes(); + w.write_all(&(sid.len() as u16).to_le_bytes()).map_err(io)?; + w.write_all(sid).map_err(io)?; + + for g in 0..self.n_genes() { + w.write_all(&self.gene_len(g).to_le_bytes()).map_err(io)?; + w.write_all(&self.gene_count(g).to_le_bytes()).map_err(io)?; + } + for id in ann.gene_ids.iter().take(self.n_genes()) { + let b = id.as_bytes(); + w.write_all(&(b.len() as u16).to_le_bytes()).map_err(io)?; + w.write_all(b).map_err(io)?; + } + // Buffer the coverage block in chunks to avoid a syscall per base. + let mut buf: Vec = Vec::with_capacity(1 << 16); + for c in &self.cov { + buf.extend_from_slice(&c.load(Ordering::Relaxed).to_le_bytes()); + if buf.len() >= (1 << 16) { + w.write_all(&buf).map_err(io)?; + buf.clear(); + } + } + if !buf.is_empty() { + w.write_all(&buf).map_err(io)?; + } + w.finish().map_err(io)?; + Ok(()) + } +} + +/// A loaded `GeneCoverage.out.bin`. +pub struct CoverageFile { + pub sample_id: String, + pub paired: bool, + pub n_counted: u64, + pub gene_ids: Vec, + pub gene_lens: Vec, + pub counts: Vec, + /// Prefix sums of `gene_lens`; `len == n_genes + 1`. + pub offsets: Vec, + pub cov: Vec, +} + +impl CoverageFile { + pub fn read(path: &Path) -> Result { + let file = std::fs::File::open(path).map_err(|e| Error::io(e, path))?; + let mut r = GzDecoder::new(BufReader::new(file)); + let io = |e: std::io::Error| Error::io(e, path); + + let mut magic = [0u8; 8]; + r.read_exact(&mut magic).map_err(io)?; + if &magic != COV_MAGIC { + return Err(Error::Parameter(format!( + "{} is not a rustar-aligner GeneCoverage file", + path.display() + ))); + } + let mut b4 = [0u8; 4]; + let mut b8 = [0u8; 8]; + let mut b2 = [0u8; 2]; + + r.read_exact(&mut b4).map_err(io)?; + let version = u32::from_le_bytes(b4); + if version != COV_VERSION { + return Err(Error::Parameter(format!( + "{} has GeneCoverage version {version}, expected {COV_VERSION}", + path.display() + ))); + } + r.read_exact(&mut b4).map_err(io)?; + let paired = u32::from_le_bytes(b4) != 0; + r.read_exact(&mut b8).map_err(io)?; + let n_genes = u64::from_le_bytes(b8) as usize; + r.read_exact(&mut b8).map_err(io)?; + let total_len = u64::from_le_bytes(b8) as usize; + r.read_exact(&mut b8).map_err(io)?; + let n_counted = u64::from_le_bytes(b8); + + r.read_exact(&mut b2).map_err(io)?; + let mut sid = vec![0u8; u16::from_le_bytes(b2) as usize]; + r.read_exact(&mut sid).map_err(io)?; + let sample_id = String::from_utf8_lossy(&sid).into_owned(); + + let mut gene_lens = Vec::with_capacity(n_genes); + let mut counts = Vec::with_capacity(n_genes); + for _ in 0..n_genes { + r.read_exact(&mut b4).map_err(io)?; + gene_lens.push(u32::from_le_bytes(b4)); + r.read_exact(&mut b4).map_err(io)?; + counts.push(u32::from_le_bytes(b4)); + } + let mut gene_ids = Vec::with_capacity(n_genes); + for _ in 0..n_genes { + r.read_exact(&mut b2).map_err(io)?; + let mut b = vec![0u8; u16::from_le_bytes(b2) as usize]; + r.read_exact(&mut b).map_err(io)?; + gene_ids.push(String::from_utf8_lossy(&b).into_owned()); + } + let mut raw = vec![0u8; total_len * 4]; + r.read_exact(&mut raw).map_err(io)?; + let cov: Vec = raw + .chunks_exact(4) + .map(|c| u32::from_le_bytes([c[0], c[1], c[2], c[3]])) + .collect(); + + let mut offsets = Vec::with_capacity(n_genes + 1); + let mut acc = 0u64; + offsets.push(0); + for &l in &gene_lens { + acc += u64::from(l); + offsets.push(acc); + } + + Ok(CoverageFile { + sample_id, + paired, + n_counted, + gene_ids, + gene_lens, + counts, + offsets, + cov, + }) + } + + pub fn gene(&self, g: usize) -> &[u32] { + &self.cov[self.offsets[g] as usize..self.offsets[g + 1] as usize] + } + + pub fn n_genes(&self) -> usize { + self.gene_ids.len() + } +} + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +#[cfg(test)] +mod tests { + use super::*; + use crate::align::transcript::{Exon, Transcript}; + + /// One two-exon gene on chr1: `[100,110)` and `[200,205)`, transcript length 15. + /// + /// Built through the real GTF path so the overlap segment trees are present. + fn two_exon_ann() -> GeneAnnotation { + let genome = crate::genome::Genome { + transform_blocks: None, + sequence: vec![0u8; 2000].into(), + n_genome: 2000, + n_genome_real: 2000, + n_chr_real: 2, + chr_start: vec![0, 1000, 2000], + chr_length: vec![1000, 1000], + chr_name: vec!["chr1".to_string(), "chr2".to_string()], + }; + let exon = |start: u64, end: u64| { + let mut attrs = std::collections::HashMap::new(); + attrs.insert("gene_id".to_string(), "G1".to_string()); + attrs.insert("transcript_id".to_string(), "T1".to_string()); + crate::junction::gtf::GtfRecord { + seqname: "chr1".to_string(), + feature: "exon".to_string(), + start, + end, + strand: '+', + attributes: attrs, + } + }; + // GTF is 1-based inclusive: 101..110 -> [100, 110), 201..205 -> [200, 205). + GeneAnnotation::from_gtf_exons(&[exon(101, 110), exon(201, 205)], &genome) + } + + fn tr(blocks: &[(u64, u64)]) -> Transcript { + Transcript { + chr_idx: 0, + genome_start: blocks[0].0, + genome_end: blocks[blocks.len() - 1].1, + is_reverse: false, + exons: blocks + .iter() + .map(|&(s, e)| Exon { + genome_start: s, + genome_end: e, + read_start: 0, + read_end: (e - s) as usize, + i_frag: 0, + }) + .collect(), + cigar: Vec::new(), + score: 0, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: Vec::new(), + junction_annotated: Vec::new(), + } + } + + #[test] + fn block_inside_first_exon_maps_to_transcript_prefix() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + assert_eq!(cov.gene_len(0), 15); + cov.add_block(0, 102, 105, &ann); + assert_eq!( + cov.gene_slice(0), + vec![0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + ); + } + + #[test] + fn block_spanning_intron_covers_both_exons_only() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.add_block(0, 105, 202, &ann); + assert_eq!( + cov.gene_slice(0), + vec![0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0] + ); + } + + #[test] + fn block_outside_all_exons_is_ignored() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.add_block(0, 120, 150, &ann); + assert_eq!(cov.gene_slice(0).iter().sum::(), 0); + } + + #[test] + fn multimapper_is_not_counted() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.count_se_read(&[tr(&[(100, 105)]), tr(&[(100, 105)])], &ann); + assert_eq!(cov.gene_slice(0).iter().sum::(), 0); + assert_eq!(cov.gene_count(0), 0); + assert_eq!(cov.n_counted(), 0); + } + + #[test] + fn unique_read_increments_coverage_and_count() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.count_se_read(&[tr(&[(100, 105)])], &ann); + assert_eq!(cov.gene_slice(0)[..5], [1, 1, 1, 1, 1]); + assert_eq!(cov.gene_count(0), 1); + assert_eq!(cov.n_counted(), 1); + } + + #[test] + fn overlapping_mates_are_counted_once_per_base() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + let pair = PairedAlignment { + mate1_transcript: tr(&[(100, 106)]), + mate2_transcript: tr(&[(104, 110)]), + mate1_region: (0, 6), + mate2_region: (0, 6), + is_proper_pair: true, + insert_size: 10, + combined_wt_score: 0, + combined_n_match: 12, + }; + cov.count_pe_read(&[&pair], &ann); + assert_eq!(cov.gene_slice(0)[..10], [1u32; 10]); + assert_eq!(cov.gene_count(0), 1); + assert_eq!(cov.n_counted(), 1); + } + + #[test] + fn coverage_file_round_trips() { + let ann = two_exon_ann(); + let cov = GeneCoverage::new(&ann); + cov.count_se_read(&[tr(&[(100, 105)])], &ann); + + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("GeneCoverage.out.bin"); + cov.write_file(&path, &ann, "sampleA", false).unwrap(); + + let f = CoverageFile::read(&path).unwrap(); + assert_eq!(f.sample_id, "sampleA"); + assert!(!f.paired); + assert_eq!(f.n_counted, 1); + assert_eq!(f.gene_ids, vec!["G1".to_string()]); + assert_eq!(f.gene_lens, vec![15]); + assert_eq!(f.counts, vec![1]); + assert_eq!(f.gene(0), &cov.gene_slice(0)[..]); + } + + #[test] + fn bad_magic_is_an_error_not_a_panic() { + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("junk.bin"); + std::fs::write(&path, b"not a coverage file").unwrap(); + assert!(CoverageFile::read(&path).is_err()); + } +} diff --git a/src/quant/mod.rs b/src/quant/mod.rs index 7df5ddc..5d87678 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -8,6 +8,7 @@ /// Submodules: /// - `transcriptome` — transcript-level alignment projection for /// `--quantMode TranscriptomeSAM` (Salmon / RSEM input). +pub mod coverage; pub mod transcriptome; use std::io::Write as _; @@ -545,6 +546,9 @@ impl GeneCounts { pub struct QuantContext { pub gene_ann: GeneAnnotation, pub counts: GeneCounts, + /// Per-gene exonic coverage, present only when `--quantMode GeneCoverage` + /// was requested (DegNorm phase 1; not a STAR feature). + pub coverage: Option, } impl QuantContext { @@ -555,13 +559,23 @@ impl QuantContext { feature_exon: &str, chr_prefix: &str, gene_tag: &str, + with_coverage: bool, ) -> Result { let exons = crate::junction::gtf::parse_gtf_configured(gtf_path, feature_exon, chr_prefix)?; let gene_ann = GeneAnnotation::from_gtf_exons_configured(&exons, genome, gene_tag); let n = gene_ann.n_genes(); log::info!("quantMode GeneCounts: {n} genes loaded from GTF"); let counts = GeneCounts::new(n); - Ok(QuantContext { gene_ann, counts }) + let coverage = if with_coverage { + Some(coverage::GeneCoverage::new(&gene_ann)) + } else { + None + }; + Ok(QuantContext { + gene_ann, + counts, + coverage, + }) } } diff --git a/tests/degnorm.rs b/tests/degnorm.rs new file mode 100644 index 0000000..0118e4d --- /dev/null +++ b/tests/degnorm.rs @@ -0,0 +1,284 @@ +//! End-to-end DegNorm: align two synthetic samples with +//! `--quantMode GeneCoverage`, then merge them with `--runMode degNorm`. +//! +//! Sample `A` covers both genes evenly. Sample `B` is identical except that +//! gene `G1` is only sequenced over its 5' 40%, which is what 3' transcript +//! degradation looks like in coverage space. The DI score for `G1` in `B` must +//! therefore be clearly higher than in `A`, while `G2` stays undegraded in both. + +use assert_cmd::cargo::cargo_bin_cmd; +use std::collections::HashMap; +use std::fs; +use std::io::Write; +use std::path::{Path, PathBuf}; +use tempfile::TempDir; + +// --------------------------------------------------------------------------- +// Fixtures +// --------------------------------------------------------------------------- + +/// LCG pseudo-random sequence generator (same LCG as `tests/alignment_features.rs`). +fn lcg_seq(seed: u32, length: usize) -> Vec { + let bases: [u8; 4] = *b"ACGT"; + let mut state = seed; + let mut seq = Vec::with_capacity(length); + for _ in 0..length { + state = state.wrapping_mul(1_103_515_245).wrapping_add(12345); + seq.push(bases[((state >> 16) & 3) as usize]); + } + seq +} + +/// Single-exon genes, 0-based half-open: G1 = [2000, 4000), G2 = [6000, 8000). +const G1: (usize, usize) = (2000, 4000); +const G2: (usize, usize) = (6000, 8000); + +fn write_fasta(dir: &Path, genome: &[u8]) -> PathBuf { + let path = dir.join("genome.fa"); + let mut f = fs::File::create(&path).unwrap(); + writeln!(f, ">chr1").unwrap(); + f.write_all(genome).unwrap(); + writeln!(f).unwrap(); + path +} + +fn write_gtf(dir: &Path) -> PathBuf { + let path = dir.join("genes.gtf"); + let mut f = fs::File::create(&path).unwrap(); + for (gene, (start, end)) in [("G1", G1), ("G2", G2)] { + writeln!( + f, + "chr1\ttest\texon\t{}\t{}\t.\t+\t.\tgene_id \"{gene}\"; transcript_id \"{gene}_T1\";", + start + 1, + end + ) + .unwrap(); + } + path +} + +fn build_index(fasta: &Path, genome_dir: &Path) { + fs::create_dir_all(genome_dir).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "genomeGenerate", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--genomeFastaFiles", + fasta.to_str().unwrap(), + "--genomeSAindexNbases", + "7", + ]) + .assert() + .success(); +} + +/// Tile `[from, to)` with 100 bp reads every `stride` bases. +fn tile_reads(f: &mut fs::File, genome: &[u8], label: &str, from: usize, to: usize, stride: usize) { + let read_len = 100; + let mut pos = from; + let mut i = 0; + while pos + read_len <= to { + writeln!(f, "@{label}_{i}").unwrap(); + f.write_all(&genome[pos..pos + read_len]).unwrap(); + writeln!(f).unwrap(); + writeln!(f, "+").unwrap(); + writeln!(f, "{}", "I".repeat(read_len)).unwrap(); + pos += stride; + i += 1; + } +} + +/// FASTQ for one sample. When `degrade_g1` is set, G1's reads only come from +/// its 5' 40%. +fn write_fastq(dir: &Path, genome: &[u8], name: &str, degrade_g1: bool) -> PathBuf { + let path = dir.join(format!("{name}.fq")); + let mut f = fs::File::create(&path).unwrap(); + let g1_end = if degrade_g1 { + G1.0 + (G1.1 - G1.0) * 2 / 5 + } else { + G1.1 + }; + tile_reads(&mut f, genome, &format!("{name}g1"), G1.0, g1_end, 25); + tile_reads(&mut f, genome, &format!("{name}g2"), G2.0, G2.1, 25); + path +} + +fn align(genome_dir: &Path, fastq: &Path, prefix: &str, quant_mode: &[&str]) { + let mut cmd = cargo_bin_cmd!("rustar-aligner"); + cmd.args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + fastq.to_str().unwrap(), + "--outFileNamePrefix", + prefix, + "--quantMode", + ]); + cmd.args(quant_mode); + cmd.arg("--sjdbGTFfile"); + cmd.arg(genome_dir.parent().unwrap().join("genes.gtf")); + cmd.assert().success(); +} + +/// Parse a `genesample...` matrix into `gene -> (sample -> value)`. +fn read_matrix(path: &Path) -> (Vec, HashMap>) { + let text = fs::read_to_string(path).unwrap(); + let mut lines = text.lines(); + let header: Vec = lines + .next() + .unwrap() + .split('\t') + .skip(1) + .map(String::from) + .collect(); + let mut rows = HashMap::new(); + for line in lines { + let mut it = line.split('\t'); + let gene = it.next().unwrap().to_string(); + rows.insert(gene, it.map(|v| v.parse::().unwrap()).collect()); + } + (header, rows) +} + +// --------------------------------------------------------------------------- +// Tests +// --------------------------------------------------------------------------- + +#[test] +fn degnorm_end_to_end_flags_the_degraded_sample() { + let tmpdir = TempDir::new().unwrap(); + let root = tmpdir.path(); + let genome = lcg_seq(88888, 20000); + let fasta = write_fasta(root, &genome); + write_gtf(root); + let genome_dir = root.join("genome"); + build_index(&fasta, &genome_dir); + + let fq_a = write_fastq(root, &genome, "A", false); + let fq_b = write_fastq(root, &genome, "B", true); + + let prefix_a = format!("{}/A_", root.display()); + let prefix_b = format!("{}/B_", root.display()); + align( + &genome_dir, + &fq_a, + &prefix_a, + &["GeneCounts", "GeneCoverage"], + ); + align( + &genome_dir, + &fq_b, + &prefix_b, + &["GeneCounts", "GeneCoverage"], + ); + + let cov_a = root.join("A_GeneCoverage.out.bin"); + let cov_b = root.join("B_GeneCoverage.out.bin"); + assert!(cov_a.exists(), "A_GeneCoverage.out.bin not written"); + assert!(cov_b.exists(), "B_GeneCoverage.out.bin not written"); + + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "degNorm", + "--degNormCoverageFiles", + cov_a.to_str().unwrap(), + cov_b.to_str().unwrap(), + "--degNormIter", + "2", + "--degNormNmfIter", + "50", + "--degNormMinHighCoverage", + "20", + "--outFileNamePrefix", + &format!("{}/", root.display()), + ]) + .assert() + .success(); + + let out_dir = root.join("DegNorm.out"); + for name in [ + "DegradationIndex.tab", + "AdjustedCounts.tab", + "RawCounts.tab", + "ScaleFactors.tab", + "Summary.txt", + ] { + assert!(out_dir.join(name).exists(), "{name} missing"); + } + + let (samples, di) = read_matrix(&out_dir.join("DegradationIndex.tab")); + assert_eq!(samples, vec!["A".to_string(), "B".to_string()]); + + let g1 = di.get("G1").expect("G1 missing from DI table"); + let g2 = di.get("G2").expect("G2 missing from DI table"); + assert!( + g1[1] > g1[0] + 0.2, + "G1 DI should be much higher in the degraded sample: A={}, B={}", + g1[0], + g1[1] + ); + assert!( + g2[0] < 0.1 && g2[1] < 0.1, + "G2 is undegraded in both samples but got DI A={}, B={}", + g2[0], + g2[1] + ); + + // Degradation-adjusted counts scale the degraded sample's G1 count up. + let (_, adjusted) = read_matrix(&out_dir.join("AdjustedCounts.tab")); + let (_, raw) = read_matrix(&out_dir.join("RawCounts.tab")); + assert!(raw["G1"][1] > 0.0, "G1 has no reads in sample B"); + assert!( + adjusted["G1"][1] > raw["G1"][1], + "adjusted count {} should exceed raw count {}", + adjusted["G1"][1], + raw["G1"][1] + ); +} + +#[test] +fn gene_coverage_does_not_change_alignment_output() { + let tmpdir = TempDir::new().unwrap(); + let root = tmpdir.path(); + let genome = lcg_seq(88888, 20000); + let fasta = write_fasta(root, &genome); + write_gtf(root); + let genome_dir = root.join("genome"); + build_index(&fasta, &genome_dir); + let fq = write_fastq(root, &genome, "A", false); + + let with_prefix = format!("{}/with_", root.display()); + let without_prefix = format!("{}/without_", root.display()); + align(&genome_dir, &fq, &without_prefix, &["GeneCounts"]); + align( + &genome_dir, + &fq, + &with_prefix, + &["GeneCounts", "GeneCoverage"], + ); + + // Compare alignment records only: the @PG header records the command line, + // which necessarily differs by the extra --quantMode word. + let records = |p: PathBuf| -> Vec { + fs::read_to_string(p) + .unwrap() + .lines() + .filter(|l| !l.starts_with('@')) + .map(String::from) + .collect() + }; + let a = records(root.join("without_Aligned.out.sam")); + let b = records(root.join("with_Aligned.out.sam")); + assert!(!a.is_empty(), "no alignments produced"); + assert_eq!(a, b, "GeneCoverage changed the SAM output"); + + // ReadsPerGene.out.tab must be identical too. + let counts_a = fs::read(root.join("without_ReadsPerGene.out.tab")).unwrap(); + let counts_b = fs::read(root.join("with_ReadsPerGene.out.tab")).unwrap(); + assert_eq!(counts_a, counts_b, "GeneCoverage changed ReadsPerGene"); +}