diff --git a/Cargo.lock b/Cargo.lock index da20215e..f84de282 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -385,9 +385,9 @@ dependencies = [ [[package]] name = "crossbeam-epoch" -version = "0.9.18" +version = "0.9.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +checksum = "2d6914041f254d6e9176c01941b21115dcfb7089e55135a35411081bd106ef3f" dependencies = [ "crossbeam-utils", ] @@ -1037,9 +1037,9 @@ dependencies = [ [[package]] name = "libz-sys" -version = "1.1.25" +version = "1.1.31" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d52f4c29e2a68ac30c9087e1b772dc9f44a2b66ed44edf2266cf2be9b03dafc1" +checksum = "ff20308ecabd635a2dc1dd5834a86c462cca20efba1d5d3ddab1f8aea3455cb2" dependencies = [ "cc", "cmake", diff --git a/docs/src/content/docs/rna/rseqc.mdx b/docs/src/content/docs/rna/rseqc.mdx index 0cf64346..81fa8bc4 100644 --- a/docs/src/content/docs/rna/rseqc.mdx +++ b/docs/src/content/docs/rna/rseqc.mdx @@ -147,15 +147,35 @@ at the end of the run suggesting the correct value: -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 diff --git a/src/cli.rs b/src/cli.rs index 6e6459e8..a09f73cd 100644 --- a/src/cli.rs +++ b/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", diff --git a/src/config.rs b/src/config.rs index 4952a1fd..ca33893a 100644 --- a/src/config.rs +++ b/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, } diff --git a/src/main.rs b/src/main.rs index 4c66c173..3bcbd86a 100644 --- a/src/main.rs +++ b/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 diff --git a/src/rna/dupradar/counting.rs b/src/rna/dupradar/counting.rs index 785115c5..45b97fd5 100644 --- a/src/rna/dupradar/counting.rs +++ b/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> { 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, + chrom_prefix: Option<&str>, + model: &GeneModel, + config: &RseqcConfig, +) -> Result { + let chromosomes: Vec = 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, diff --git a/src/rna/rseqc/accumulators.rs b/src/rna/rseqc/accumulators.rs index b91a409e..f11b220f 100644 --- a/src/rna/rseqc/accumulators.rs +++ b/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. diff --git a/src/rna/rseqc/infer_experiment.rs b/src/rna/rseqc/infer_experiment.rs index bb257881..84e985de 100644 --- a/src/rna/rseqc/infer_experiment.rs +++ b/src/rna/rseqc/infer_experiment.rs @@ -201,8 +201,6 @@ pub fn write_infer_experiment>( )?; } _ => { - 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 { diff --git a/tests/README.infer-sampling.md b/tests/README.infer-sampling.md new file mode 100644 index 00000000..8bae0805 --- /dev/null +++ b/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. diff --git a/tests/data/infer_exp_sampling.yaml b/tests/data/infer_exp_sampling.yaml new file mode 100644 index 00000000..37a5a6c8 --- /dev/null +++ b/tests/data/infer_exp_sampling.yaml @@ -0,0 +1,42 @@ +rna: + dupradar: + dup_matrix: false + intercept_slope: false + density_scatter_plot: false + boxplot: false + expression_histogram: false + multiqc_intercept: false + multiqc_curve: false + featurecounts: + counts_file: false + summary_file: false + biotype_summary_file: false + biotype_counts: false + biotype_counts_mqc: false + biotype_rrna_mqc: false + bam_stat: + enabled: false + infer_experiment: + enabled: true + read_duplication: + enabled: false + read_distribution: + enabled: false + junction_annotation: + enabled: false + junction_saturation: + enabled: false + inner_distance: + enabled: false + flagstat: + enabled: true + idxstats: + enabled: false + tin: + enabled: false + samtools_stats: + enabled: true + preseq: + enabled: false + qualimap: + enabled: false diff --git a/tests/infer_sampling.rs b/tests/infer_sampling.rs new file mode 100644 index 00000000..ced3f70a --- /dev/null +++ b/tests/infer_sampling.rs @@ -0,0 +1,543 @@ +//! Native regression controls for explicit, file-order strandedness sampling. +use std::fs; +use std::path::{Path, PathBuf}; +use std::process::Command; +use std::time::{SystemTime, UNIX_EPOCH}; + +use rust_htslib::bam; +use rust_htslib::bam::header::HeaderRecord; +use rust_htslib::bam::HeaderView; + +fn root(label: &str) -> PathBuf { + let nonce = SystemTime::now() + .duration_since(UNIX_EPOCH) + .unwrap() + .as_nanos(); + let path = std::env::temp_dir().join(format!( + "rustqc-sampling-{label}-{}-{nonce}", + std::process::id() + )); + fs::create_dir(&path).unwrap(); + path +} + +fn write_alignment( + path: &Path, + chroms: &[(&str, u64)], + rows: &[(String, u64)], + format: bam::Format, + reference: Option<&Path>, +) { + let mut header = bam::Header::new(); + header.push_record( + HeaderRecord::new(b"HD") + .push_tag(b"VN", &"1.6") + .push_tag(b"SO", &"coordinate"), + ); + for (chrom, length) in chroms { + header.push_record( + HeaderRecord::new(b"SQ") + .push_tag(b"SN", chrom) + .push_tag(b"LN", length), + ); + } + let view = HeaderView::from_header(&header); + let mut writer = bam::Writer::from_path(path, &header, format).unwrap(); + if let Some(reference) = reference { + writer.set_reference(reference).unwrap(); + } + for (fields, count) in rows { + let line = format!("{fields}\t{}\t{}", "A".repeat(50), "I".repeat(50)); + let record = bam::Record::from_sam(&view, line.as_bytes()).unwrap(); + for _ in 0..*count { + writer.write(&record).unwrap(); + } + } + drop(writer); + bam::index::build(path, None, bam::index::Type::Bai, 1).unwrap(); +} + +fn write_bam(path: &Path, chroms: &[(&str, u64)], rows: &[(String, u64)]) { + write_alignment(path, chroms, rows, bam::Format::Bam, None); +} + +fn run( + root: &Path, + threads: usize, + cli_cap: Option, + yaml_cap: Option, + additions: &str, + expected: [f64; 3], +) -> String { + let folder = root.join(format!("t{threads}-cli{cli_cap:?}-yaml{yaml_cap:?}")); + fs::create_dir(&folder).unwrap(); + let mut config = include_str!("data/infer_exp_sampling.yaml").to_string(); + if let Some(cap) = yaml_cap { + config = config.replace( + " infer_experiment:\n enabled: true\n", + &format!(" infer_experiment:\n enabled: true\n sample_size: {cap}\n"), + ); + } + config.push_str(additions); + let config_path = folder.join("config.yaml"); + fs::write(&config_path, config).unwrap(); + let mut command = Command::new(env!("CARGO_BIN_EXE_rustqc")); + command.env_remove("RUSTQC_INFER_EXPERIMENT_SAMPLE_SIZE"); + command + .args(["rna", "-q", "--skip-dup-check", "--config"]) + .arg(&config_path) + .arg("--gtf") + .arg(root.join("a.gtf")) + .arg("--outdir") + .arg(&folder) + .args(["--threads", &threads.to_string()]) + .arg(root.join("s.bam")); + if let Some(cap) = cli_cap { + command.args(["--infer-experiment-sample-size", &cap.to_string()]); + } + let output = command.output().unwrap(); + assert!( + output.status.success(), + "sampling run failed: {}", + String::from_utf8_lossy(&output.stderr) + ); + let text = + fs::read_to_string(folder.join("rseqc/infer_experiment/s.infer_experiment.txt")).unwrap(); + let observed: Vec = text + .lines() + .filter(|line| line.starts_with("Fraction of reads")) + .map(|line| line.rsplit_once(": ").unwrap().1.parse().unwrap()) + .collect(); + assert_eq!( + observed, expected, + "sampling t{threads} CLI={cli_cap:?} YAML={yaml_cap:?}:\n{text}" + ); + fs::read_to_string(folder.join("samtools/s.flagstat")).unwrap() +} + +fn write_gtf(path: &Path, rows: &[(&str, u64, u64, &str, &str)]) { + let mut annotation = String::new(); + for (chrom, start, end, strand, id) in rows { + annotation.push_str(&format!("{chrom}\tt\texon\t{start}\t{end}\t.\t{strand}\t.\tgene_id \"{id}\"; transcript_id \"{id}\";\n")); + } + fs::write(path, annotation).unwrap(); +} + +#[test] +fn test_explicit_caps_file_order_default_and_yaml() { + let root = root("large"); + write_gtf( + &root.join("a.gtf"), + &[("chr1", 1, 100, "+", "one"), ("chr2", 1, 100, "+", "two")], + ); + let rows: Vec<_> = [ + ("sense1", "chr1", 0, 156_000), + ("anti1", "chr1", 16, 44_000), + ("sense2", "chr2", 0, 100_000), + ] + .into_iter() + .map(|(id, chrom, flags, count)| { + ( + format!("{id}\t{flags}\t{chrom}\t11\t60\t50M\t*\t0\t0"), + count, + ) + }) + .collect(); + write_bam( + &root.join("s.bam"), + &[("chr1", 100_000), ("chr2", 200_000)], + &rows, + ); + let mut flagstat = None; + for threads in [1, 2, 4, 8] { + for (cli, yaml, expected) in [ + (None, None, [0.0, 0.8533, 0.1467]), + (Some(u64::MAX), None, [0.0, 0.8533, 0.1467]), + (Some(1), None, [0.0, 1.0, 0.0]), + (Some(200_000), None, [0.0, 0.78, 0.22]), + (Some(400_000), None, [0.0, 0.8533, 0.1467]), + (None, Some(200_000), [0.0, 0.78, 0.22]), + (Some(1), Some(200_000), [0.0, 1.0, 0.0]), + ] { + let observed = run(&root, threads, cli, yaml, "", expected); + assert!( + observed.starts_with("300000 + 0 in total"), + "sampling must not truncate other analyses" + ); + if let Some(ref baseline) = flagstat { + assert_eq!(&observed, baseline); + } else { + flagstat = Some(observed); + } + } + } + let _ = fs::remove_dir_all(root); +} + +#[test] +fn test_explicit_caps_filtering_and_ambiguous_overlap() { + let root = root("eligible"); + write_gtf( + &root.join("a.gtf"), + &[ + ("chr1", 1001, 1400, "+", "plus"), + ("chr1", 1101, 1200, "-", "minus"), + ("chr2", 1001, 2000, "-", "other"), + ], + ); + let rows: Vec<_> = [ + ("no_overlap", 0, "chr1", 11, 60), + ("qcfail", 512, "chr1", 1001, 60), + ("duplicate", 1024, "chr1", 1002, 60), + ("secondary", 256, "chr1", 1003, 60), + ("low_mapq", 0, "chr1", 1004, 29), + ("unmapped", 4, "chr1", 1005, 60), + ("ambiguous", 0, "chr1", 1110, 60), + ("sense", 0, "chr1", 1250, 60), + ("supplementary", 2048, "chr1", 1251, 60), + ("antisense", 0, "chr2", 1001, 60), + ] + .into_iter() + .map(|(id, flags, chrom, pos, mapq)| { + ( + format!("{id}\t{flags}\t{chrom}\t{pos}\t{mapq}\t50M\t*\t0\t0"), + 1, + ) + }) + .collect(); + write_bam( + &root.join("s.bam"), + &[("chr1", 100_000), ("chr2", 200_000)], + &rows, + ); + let baseline = run(&root, 1, None, None, "", [0.25, 0.5, 0.25]); + for threads in [1, 2, 4, 8] { + for (cap, expected) in [ + (1, [1.0, 0.0, 0.0]), + (2, [0.5, 0.5, 0.0]), + (3, [0.3333, 0.6667, 0.0]), + (10, [0.25, 0.5, 0.25]), + ] { + let flagstat = run(&root, threads, Some(cap), None, "", expected); + assert_eq!( + flagstat, baseline, + "sampling must preserve uncapped flagstat" + ); + } + } + let _ = fs::remove_dir_all(root); +} + +#[test] +fn test_explicit_caps_resolve_chromosome_names() { + for additions in [ + " chromosome_prefix: chr\n", + " chromosome_mapping:\n chr1: \"1\"\n chromosome_prefix: ignored\n", + ] { + let root = root("chromosome-names"); + write_gtf(&root.join("a.gtf"), &[("chr1", 1, 100, "+", "one")]); + write_bam( + &root.join("s.bam"), + &[("1", 100_000)], + &[ + ("sense\t0\t1\t11\t60\t50M\t*\t0\t0".into(), 1), + ("anti\t16\t1\t12\t60\t50M\t*\t0\t0".into(), 1), + ], + ); + for threads in [1, 2, 4, 8] { + assert!( + run(&root, threads, Some(1), None, additions, [0.0, 1.0, 0.0]) + .starts_with("2 + 0 in total") + ); + } + let _ = fs::remove_dir_all(root); + } +} + +#[test] +fn test_explicit_caps_keep_cram_reference_configuration() { + let root = root("cram-reference"); + write_gtf( + &root.join("a.gtf"), + &[("chr1", 1, 1000, "+", "one"), ("chr2", 1, 1000, "+", "two")], + ); + let reference = root.join("reference.fa"); + fs::write( + &reference, + format!(">chr1\n{}\n>chr2\n{}\n", "A".repeat(1000), "A".repeat(1000)), + ) + .unwrap(); + rust_htslib::faidx::build(&reference).unwrap(); + let rows: Vec<_> = [ + "first\t0\tchr1\t101\t60\t50M\t*\t0\t0", + "anti1\t16\tchr1\t201\t60\t50M\t*\t0\t0", + "anti2\t16\tchr2\t101\t60\t50M\t*\t0\t0", + "anti3\t16\tchr2\t151\t60\t50M\t*\t0\t0", + "anti4\t16\tchr2\t201\t60\t50M\t*\t0\t0", + "anti5\t16\tchr2\t251\t60\t50M\t*\t0\t0", + ] + .into_iter() + .map(|line| (line.to_owned(), 1)) + .collect(); + let cram = root.join("s.cram"); + write_alignment( + &cram, + &[("chr1", 1000), ("chr2", 1000)], + &rows, + bam::Format::Cram, + Some(&reference), + ); + let config = root.join("config.yaml"); + fs::write(&config, include_str!("data/infer_exp_sampling.yaml")).unwrap(); + for threads in [1, 2, 4, 8] { + let outdir = root.join(format!("t{threads}")); + let output = Command::new(env!("CARGO_BIN_EXE_rustqc")) + .env_remove("RUSTQC_INFER_EXPERIMENT_SAMPLE_SIZE") + .args(["rna", "-q", "--skip-dup-check", "--config"]) + .arg(&config) + .arg("--gtf") + .arg(root.join("a.gtf")) + .arg("--outdir") + .arg(&outdir) + .args([ + "--threads", + &threads.to_string(), + "--infer-experiment-sample-size", + "1", + "--reference", + ]) + .arg(&reference) + .arg(&cram) + .output() + .unwrap(); + assert!( + output.status.success(), + "CRAM sampling failed: {}", + String::from_utf8_lossy(&output.stderr) + ); + let text = fs::read_to_string(outdir.join("rseqc/infer_experiment/s.infer_experiment.txt")) + .unwrap(); + assert!(text.contains("This is SingleEnd Data")); + assert!( + text.contains("Fraction of reads explained by \"++,--\": 1.0000"), + "{text}" + ); + assert!( + text.contains("Fraction of reads explained by \"+-,-+\": 0.0000"), + "{text}" + ); + assert!(fs::read_to_string(outdir.join("samtools/s.flagstat")) + .unwrap() + .starts_with("6 + 0 in total")); + } + let _ = fs::remove_dir_all(root); +} + +#[test] +fn test_explicit_caps_sample_paired_records_in_file_order() { + let root = root("paired-prefix"); + write_gtf( + &root.join("a.gtf"), + &[("chr1", 1, 1000, "+", "one"), ("chr2", 1, 1000, "+", "two")], + ); + let rows: Vec<_> = [ + "a\t99\tchr1\t101\t60\t50M\t=\t151\t100", + "a\t147\tchr1\t151\t60\t50M\t=\t101\t-100", + "b\t163\tchr1\t201\t60\t50M\t=\t251\t100", + "b\t83\tchr1\t251\t60\t50M\t=\t201\t-100", + "c\t163\tchr2\t101\t60\t50M\t=\t151\t100", + "c\t83\tchr2\t151\t60\t50M\t=\t101\t-100", + ] + .into_iter() + .map(|line| (line.to_string(), 1)) + .collect(); + write_bam( + &root.join("s.bam"), + &[("chr1", 100_000), ("chr2", 200_000)], + &rows, + ); + let config = root.join("config.yaml"); + fs::write(&config, include_str!("data/infer_exp_sampling.yaml")).unwrap(); + for threads in [1, 4] { + for (cap, expected) in [ + (1, [0.0, 1.0, 0.0]), + (3, [0.0, 0.6667, 0.3333]), + (10, [0.0, 0.3333, 0.6667]), + ] { + let outdir = root.join(format!("t{threads}-cap{cap}")); + let output = Command::new(env!("CARGO_BIN_EXE_rustqc")) + .env_remove("RUSTQC_INFER_EXPERIMENT_SAMPLE_SIZE") + .args(["rna", "-q", "--skip-dup-check", "--paired", "--config"]) + .arg(&config) + .arg("--gtf") + .arg(root.join("a.gtf")) + .arg("--outdir") + .arg(&outdir) + .args([ + "--threads", + &threads.to_string(), + "--infer-experiment-sample-size", + &cap.to_string(), + ]) + .arg(root.join("s.bam")) + .output() + .unwrap(); + assert!( + output.status.success(), + "paired sampling failed: {}", + String::from_utf8_lossy(&output.stderr) + ); + let text = + fs::read_to_string(outdir.join("rseqc/infer_experiment/s.infer_experiment.txt")) + .unwrap(); + assert!(text.contains("This is PairEnd Data")); + let actual: Vec = text + .lines() + .filter(|line| line.starts_with("Fraction of reads")) + .map(|line| line.rsplit_once(": ").unwrap().1.parse().unwrap()) + .collect(); + assert_eq!(actual, expected, "t{threads} cap{cap}:\n{text}"); + assert!(fs::read_to_string(outdir.join("samtools/s.flagstat")) + .unwrap() + .starts_with("6 + 0 in total")); + } + } + let _ = fs::remove_dir_all(root); +} + +#[test] +fn test_explicit_zero_cap_retains_unknown_frame() { + let root = root("zero-frame"); + write_gtf(&root.join("a.gtf"), &[("chr1", 1, 1000, "+", "one")]); + write_bam( + &root.join("s.bam"), + &[("chr1", 100_000)], + &[("sense\t0\tchr1\t11\t60\t50M\t*\t0\t0".into(), 1)], + ); + for threads in [1, 4] { + for cli in [false, true] { + let folder = root.join(format!("t{threads}-cli{cli}")); + fs::create_dir(&folder).unwrap(); + let mut config = include_str!("data/infer_exp_sampling.yaml").to_string(); + if !cli { + config = config.replace( + " infer_experiment:\n enabled: true\n", + " infer_experiment:\n enabled: true\n sample_size: 0\n", + ); + } + let config_path = folder.join("config.yaml"); + fs::write(&config_path, config).unwrap(); + let mut command = Command::new(env!("CARGO_BIN_EXE_rustqc")); + command.env_remove("RUSTQC_INFER_EXPERIMENT_SAMPLE_SIZE"); + command + .args(["rna", "-q", "--skip-dup-check", "--config"]) + .arg(&config_path) + .arg("--gtf") + .arg(root.join("a.gtf")) + .arg("--outdir") + .arg(&folder) + .args(["--threads", &threads.to_string()]) + .arg(root.join("s.bam")); + if cli { + command.args(["--infer-experiment-sample-size", "0"]); + } + let output = command.output().unwrap(); + assert!( + output.status.success(), + "zero-cap run failed: {}", + String::from_utf8_lossy(&output.stderr) + ); + assert_eq!( + fs::read(folder.join("rseqc/infer_experiment/s.infer_experiment.txt")).unwrap(), + b"Unknown Data type\n".to_vec() + ); + assert!(fs::read_to_string(folder.join("samtools/s.flagstat")) + .unwrap() + .starts_with("1 + 0 in total")); + } + } + let _ = fs::remove_dir_all(root); +} + +#[test] +fn test_explicit_caps_reset_for_each_input_in_either_order() { + let root = root("multiple-inputs"); + write_gtf(&root.join("a.gtf"), &[("chr1", 1, 1000, "+", "one")]); + let sense = root.join("sense_first.bam"); + let anti = root.join("anti_first.bam"); + for (path, first, rest) in [(&sense, 0, 16), (&anti, 16, 0)] { + let rows: Vec<_> = [first, rest, rest, rest] + .into_iter() + .enumerate() + .map(|(index, flags)| { + ( + format!( + "r{index}\t{flags}\tchr1\t{}\t60\t50M\t*\t0\t0", + 101 + index * 100 + ), + 1, + ) + }) + .collect(); + write_bam(path, &[("chr1", 100_000), ("chr2", 200_000)], &rows); + } + let config = root.join("config.yaml"); + fs::write(&config, include_str!("data/infer_exp_sampling.yaml")).unwrap(); + for threads in [1, 4] { + for reverse in [false, true] { + let outdir = root.join(format!("t{threads}-reversed{reverse}")); + let (first, second) = if reverse { + (&anti, &sense) + } else { + (&sense, &anti) + }; + let output = Command::new(env!("CARGO_BIN_EXE_rustqc")) + .env_remove("RUSTQC_INFER_EXPERIMENT_SAMPLE_SIZE") + .args(["rna", "-q", "--skip-dup-check", "--config"]) + .arg(&config) + .arg("--gtf") + .arg(root.join("a.gtf")) + .arg("--outdir") + .arg(&outdir) + .args([ + "--threads", + &threads.to_string(), + "--infer-experiment-sample-size", + "1", + ]) + .arg(first) + .arg(second) + .output() + .unwrap(); + assert!( + output.status.success(), + "multiple-input run failed: {}", + String::from_utf8_lossy(&output.stderr) + ); + for (name, expected) in [ + ("sense_first", [0.0, 1.0, 0.0]), + ("anti_first", [0.0, 0.0, 1.0]), + ] { + let text = fs::read_to_string(outdir.join(format!( + "rseqc/infer_experiment/{name}.infer_experiment.txt" + ))) + .unwrap(); + let actual: Vec = text + .lines() + .filter(|line| line.starts_with("Fraction of reads")) + .map(|line| line.rsplit_once(": ").unwrap().1.parse().unwrap()) + .collect(); + assert_eq!( + actual, expected, + "{name}, t{threads}, reversed{reverse}:\n{text}" + ); + assert!( + fs::read_to_string(outdir.join(format!("samtools/{name}.flagstat"))) + .unwrap() + .starts_with("4 + 0 in total") + ); + } + } + } + let _ = fs::remove_dir_all(root); +}