fix(rseqc): honor explicit infer-experiment sampling limits

docs/src/content/docs/rna/rseqc.mdx

@@ -147,15 +147,35 @@ at the end of the run suggesting the correct value:
 
 </details>
 
-The difference is caused by a read-sampling mismatch: the upstream
-`infer_experiment.py` defaults to sampling only 200,000 reads. On a 186M-read BAM,
+The difference in this benchmark is caused by a read-sampling mismatch: the upstream
+`infer_experiment.py` defaults to sampling only 200,000 reads. On this 186M-read BAM,
 this early-file subsample draws predominantly from the first few chromosomes, which
 underrepresents loci where transcripts on both strands overlap, causing the
 "failed-to-determine" fraction to appear lower than it truly is. RustQC processes
-all reads, giving a more representative estimate.
-
-Both tools correctly identify the library as strongly reverse-stranded
-(antisense ~88-92%), so the practical interpretation is identical.
+all reads when the sample-size setting is unset, giving a more representative
+estimate for this dataset.
+
+An explicit
+`--infer-experiment-sample-size N` or `rna.infer_experiment.sample_size: N` selects
+the first N eligible transcript-overlapping records in input-file order,
+independently for each input and across thread counts. Paired-end records each
+consume a slot; an ambiguous transcript overlap also consumes a slot. Existing
+quality and alignment filters still apply.
+
+A finite cap uses a streaming prefix pass for infer-experiment while the other
+analyses retain their normal parallel pass. If too few eligible reads are
+present, this can add one full input scan. Zero samples no reads; an explicit
+`18446744073709551615` (the maximum unsigned 64-bit value) is unlimited.
+Every Unknown report, including no-data and mixed
+library reports, is exactly `Unknown Data type` followed by a newline, without
+leading blank lines. Frame selection and input success/error policy are unchanged.
+
+On this benchmark, both tools correctly identify the library as strongly
+reverse-stranded (antisense ~88-92%), so the practical interpretation is identical
+for this dataset. Different sample populations can change a downstream strand/QC
+decision; this benchmark does not guarantee the same interpretation elsewhere.
+Pipelines need a fixed RustQC artifact and an explicit finite cap to request that
+sampling policy.
 
 ## read_duplication
 

src/cli.rs

@@ -185,7 +185,7 @@ pub struct RnaArgs {
     pub verbose: bool,
 
     // ── Tool parameters ─────────────────────────────────────────────────
-    /// infer_experiment: sample size [default: 200000]
+    /// infer_experiment: maximum usable reads [default: all reads]
     #[arg(
         long = "infer-experiment-sample-size",
         value_name = "N",

src/config.rs

@@ -329,7 +329,8 @@ pub struct InferExperimentConfig {
     /// Whether to run infer_experiment analysis. Defaults to true.
     pub enabled: bool,
 
-    /// Maximum number of reads to sample for strandedness inference.
+    /// Maximum number of usable reads to sample for strandedness inference.
+    /// When unset, processes all reads.
     /// Can be overridden by `--infer-experiment-sample-size` CLI flag.
     pub sample_size: Option<u64>,
 }

src/main.rs

@@ -545,7 +545,7 @@ fn run_rna(args: cli::RnaArgs, ui: &Ui) -> Result<()> {
         exon_bitset: exon_bitset.as_ref(),
         transcript_tree: transcript_tree.as_ref(),
         mapq_cut: args.mapq_cut,
-        infer_experiment_sample_size: config.infer_experiment.sample_size.unwrap_or(200_000),
+        infer_experiment_sample_size: config.infer_experiment.sample_size.unwrap_or(u64::MAX),
         min_intron: config.junction_annotation.min_intron.unwrap_or(50),
         junction_saturation_min_coverage: config.junction_saturation.min_coverage.unwrap_or(1),
         junction_saturation_percentile_floor: config

src/rna/dupradar/counting.rs

@@ -12,7 +12,10 @@
 use crate::gtf::Gene;
 use crate::io::format_count;
 use crate::rna::qualimap::QualimapAccum;
-use crate::rna::rseqc::accumulators::{RseqcAccumulators, RseqcAnnotations, RseqcConfig};
+use crate::rna::rseqc::accumulators::{
+    InferExpAccum, RseqcAccumulators, RseqcAnnotations, RseqcConfig,
+};
+use crate::rna::rseqc::infer_experiment::GeneModel;
 use crate::Strandedness;
 use anyhow::{Context, Result};
 use coitrees::{COITree, Interval, IntervalTree};
@@ -1130,6 +1133,60 @@ fn partition_chromosomes(lengths: &[u64], num_workers: usize) -> Vec<Vec<u32>> {
     batches
 }
 
+/// Collect the first usable reads in file order, independent of chromosome batching.
+fn sample_infer_experiment(
+    bam_path: &str,
+    reference: Option<&str>,
+    tid_to_name: &[String],
+    chrom_mapping: &HashMap<String, String>,
+    chrom_prefix: Option<&str>,
+    model: &GeneModel,
+    config: &RseqcConfig,
+) -> Result<InferExpAccum> {
+    let chromosomes: Vec<String> = tid_to_name
+        .iter()
+        .map(|name| {
+            if let Some(mapped) = chrom_mapping.get(name) {
+                mapped.clone()
+            } else if let Some(prefix) = chrom_prefix {
+                format!("{prefix}{name}")
+            } else {
+                name.clone()
+            }
+        })
+        .collect();
+    let mut bam = bam::Reader::from_path(bam_path).with_context(|| {
+        format!("Failed to open alignment file for infer_experiment sampling: {bam_path}")
+    })?;
+    if let Some(path) = reference {
+        bam.set_reference(path)
+            .with_context(|| format!("Failed to set reference FASTA: {path}"))?;
+    }
+    let limit = config.infer_experiment_sample_size;
+    let mut sampled = 0;
+    let mut accum = InferExpAccum::default();
+    let mut record = bam::Record::new();
+    while sampled < limit {
+        let Some(result) = bam.read(&mut record) else {
+            break;
+        };
+        result.context("Error reading alignment record for infer_experiment sampling")?;
+        let chrom = if record.tid() >= 0 {
+            chromosomes
+                .get(record.tid() as usize)
+                .map(String::as_str)
+                .unwrap_or("")
+        } else {
+            ""
+        };
+        if accum.process_sample_read(&record, chrom, model, config.mapq_cut) {
+            sampled += 1;
+        }
+    }
+    debug!("infer_experiment sampled {sampled} usable reads in file order (cap {limit})");
+    Ok(accum)
+}
+
 ///   an index for `threads > 1`; SAM files always use single-threaded mode)
 /// * `genes` - Gene annotation map from GTF parsing
 /// * `stranded` - Library strandedness (unstranded, forward, or reverse)
@@ -1179,6 +1236,34 @@ pub fn count_reads(
         (names, lengths)
     };
 
+    // Aggregate worker maps cannot recover a global file-order prefix. Sample
+    // explicit caps separately; other analyses retain their normal parallel pass.
+    let mut worker_rseqc_config = rseqc_config.cloned();
+    let sampled_infer_exp = match (
+        rseqc_config,
+        rseqc_annotations.and_then(|annotations| annotations.gene_model),
+    ) {
+        (Some(cfg), Some(model))
+            if cfg.infer_experiment_enabled && cfg.infer_experiment_sample_size != u64::MAX =>
+        {
+            let sampled = sample_infer_experiment(
+                bam_path,
+                reference,
+                &tid_to_name,
+                chrom_mapping,
+                chrom_prefix,
+                model,
+                cfg,
+            )?;
+            if let Some(worker_config) = &mut worker_rseqc_config {
+                worker_config.infer_experiment_enabled = false;
+            }
+            Some(sampled)
+        }
+        _ => None,
+    };
+    let rseqc_config = worker_rseqc_config.as_ref();
+
     if skip_dup_check {
         log::info!("Skipping duplicate-marking verification (--skip-dup-check)");
     }
@@ -1685,6 +1770,10 @@ pub fn count_reads(
         })
         .collect();
 
+    if let (Some(accums), Some(infer_exp)) = (&mut merged_rseqc, sampled_infer_exp) {
+        accums.infer_exp = Some(infer_exp);
+    }
+
     Ok(CountResult {
         gene_counts: gene_counts_map,
         stat_total_reads: merged.total_reads,

src/rna/rseqc/accumulators.rs

@@ -61,7 +61,7 @@ pub struct RseqcAnnotations<'a> {
 pub struct RseqcConfig {
     /// MAPQ cutoff for read quality filtering.
     pub mapq_cut: u8,
-    /// Maximum reads to sample for infer_experiment.
+    /// Maximum infer_experiment reads; u64::MAX retains unlimited sampling.
     pub infer_experiment_sample_size: u64,
     /// Minimum intron size for junction filtering.
     pub min_intron: u64,
@@ -1316,10 +1316,8 @@ impl BamStatAccum {
 
 /// infer_experiment accumulator — strand protocol inference.
 ///
-/// Processes ALL reads (no sampling cap). Upstream RSeQC `infer_experiment.py`
-/// samples 200K reads, but since RustQC processes in a single pass, processing
-/// all reads gives equivalent-or-better results on small data and avoids
-/// sampling divergence on large data.
+/// Explicit caps are collected in a file-order prepass before chromosome workers
+/// run; the normal worker accumulators retain the all-read behavior when unset.
 #[derive(Debug, Default)]
 pub struct InferExpAccum {
     /// Paired-end strand class counts (keys: "1++", "1--", "2+-", "2-+", etc.)
@@ -1339,6 +1337,18 @@ impl InferExpAccum {
         model: &GeneModel,
         mapq_cut: u8,
     ) {
+        self.process_sample_read(record, chrom, model, mapq_cut);
+    }
+
+    /// Count one record and report whether it consumed a sample slot.
+    /// Ambiguous transcript overlaps consume a slot, matching native RSeQC.
+    pub(crate) fn process_sample_read(
+        &mut self,
+        record: &bam::Record,
+        chrom: &str,
+        model: &GeneModel,
+        mapq_cut: u8,
+    ) -> bool {
         let flags = record.flags();
 
         // Skip QC-fail, dup, secondary, unmapped.
@@ -1349,11 +1359,11 @@ impl InferExpAccum {
             || flags & BAM_FSECONDARY != 0
             || flags & BAM_FUNMAP != 0
         {
-            return;
+            return false;
         }
 
         if record.mapq() < mapq_cut {
-            return;
+            return false;
         }
 
         let map_strand = if record.is_reverse() { '-' } else { '+' };
@@ -1381,7 +1391,7 @@ impl InferExpAccum {
 
         let strands = model.find_strands(chrom, read_start, read_end);
         if strands.is_empty() {
-            return;
+            return false;
         }
 
         // Build key in reusable buffer to avoid per-read allocation.
@@ -1411,6 +1421,7 @@ impl InferExpAccum {
                 map.insert(self.key_buf.clone(), 1);
             }
         }
+        true
     }
 
     /// Merge another accumulator into this one.

src/rna/rseqc/infer_experiment.rs

@@ -201,8 +201,6 @@ pub fn write_infer_experiment<P: AsRef<Path>>(
             )?;
         }
         _ => {
-            writeln!(writer)?;
-            writeln!(writer)?;
             writeln!(writer, "Unknown Data type")?;
         }
     }
@@ -316,6 +314,32 @@ mod tests {
     use crate::gtf::{Exon, Gene};
     use indexmap::IndexMap;
 
+    #[test]
+    fn test_write_unknown_reports_match_native_bytes() {
+        let nonce = std::time::SystemTime::now()
+            .duration_since(std::time::UNIX_EPOCH)
+            .unwrap()
+            .as_nanos();
+        let root = std::env::temp_dir().join(format!(
+            "rustqc-infer-unknown-{}-{nonce}",
+            std::process::id()
+        ));
+        fs::create_dir(&root).unwrap();
+        for (library_type, total_sampled) in [("Undetermined", 0), ("Mixture", 2)] {
+            let result = InferExperimentResult {
+                total_sampled,
+                library_type: library_type.to_owned(),
+                frac_failed: 0.0,
+                frac_protocol1: 0.5,
+                frac_protocol2: 0.5,
+            };
+            let path = root.join(format!("{library_type}.txt"));
+            write_infer_experiment(&result, &path).unwrap();
+            assert_eq!(fs::read(path).unwrap(), b"Unknown Data type\n".to_vec());
+        }
+        fs::remove_dir_all(root).unwrap();
+    }
+
     /// Helper to create a simple Gene for testing.
     fn make_gene(gene_id: &str, chrom: &str, start: u64, end: u64, strand: char) -> Gene {
         Gene {

tests/README.infer-sampling.md

@@ -0,0 +1,43 @@
+Run the native reproduction with:
+
+```sh
+cargo test --release --test infer_sampling
+```
+
+On unchanged v0.2.1/main, the explicit cap is ignored. The 300,000-record SE
+witness has 156,000 sense and 44,000 antisense records in its first contig,
+followed by 100,000 sense records on a larger second contig. Cap 1 must report
+1.0000/0.0000 and cap 200,000 must report 0.7800/0.2200; the broken program uses
+all records and reports 0.8533/0.1467. That population change can change a
+downstream strandedness/QC decision.
+
+The tests invoke Cargo's compile-time binary path. They generate indexed BAMs
+and a reference-backed indexed CRAM with rust-htslib, including the FASTA index;
+no noodles-only fixture, external executable or test dependency is required.
+Expected positive-cap fractions and the exact zero frame were checked against
+native RSeQC 5.0.4. Paired cap 3 stops midpair after three records. Other controls
+cover filtered and ambiguous reads, supplementary eligibility, CLI/YAML
+precedence, prefix/mapping resolution, zero, maximum cap, current unset/all-read
+behavior and threads 1/2/4/8. A two-input control repeats both argument orders
+to check that each file independently gets its cap. New commands clear ambient
+RUSTQC_INFER_EXPERIMENT_SAMPLE_SIZE. Full flagstat equality protects other analyses.
+
+The default is intentionally unchanged. A production pipeline benefits only
+after using a fixed artifact and requesting the intended finite cap; this PR
+does not change pipeline artifacts, flags or thresholds. Unknown-strand GTF
+eligibility and whole-pipeline equivalence are separate concerns.
+
+The fix collects an explicitly capped prefix in a streaming prepass, then runs
+the normal parallel pipeline with only infer-experiment disabled. Each input
+gets its own counter. The prepass retains one alignment and the small strand
+counters; annotation/reference memory follows the existing backend. If usable
+reads are scarce or the cap exceeds the population, it can add one full scan.
+Unset and maximum caps keep the existing all-read behavior. Explicit zero
+samples no reads. Every Unknown report, including mixed/no-data reports,
+now writes native RSeQC's exact 18-byte line without the previous two leading
+blank lines; a direct writer test covers those cases. Input success/error
+policy and frame selection are unchanged.
+
+Counter-correctness fixes are excluded. In particular, the filtered-read test
+compares complete flagstat bytes to this branch's uncapped baseline rather
+than requiring an unrelated QC-column correction.