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/src/rna/dupradar/counting.rs b/src/rna/dupradar/counting.rs index 785115c5..11535877 100644 --- a/src/rna/dupradar/counting.rs +++ b/src/rna/dupradar/counting.rs @@ -1650,6 +1650,33 @@ pub fn count_reads( .filter(|c| c.all_multi > 0) .count(); if merged.total_mapped > 0 && genes_with_reads == 0 { + // Only a naming mismatch if no alignment name, after prefix/mapping, + // is a GTF chromosome. + let shared_chroms: Vec<&str> = tid_to_name + .iter() + .filter(|bam_name| { + let gtf_name = match chrom_mapping.get(bam_name.as_str()) { + Some(mapped) => mapped.clone(), + None => format!("{}{}", chrom_prefix.unwrap_or(""), bam_name), + }; + index.contains_key(>f_name) + }) + .take(5) + .map(|s| s.as_str()) + .collect(); + if !shared_chroms.is_empty() { + anyhow::bail!( + "No reads could be assigned to any gene under the configured gene-assignment rules \ + ({} mapped reads).\n\ + \n\ + Chromosomes in both the alignment and the GTF (first 5): {}\n\ + \n\ + The chromosome names match. Check the annotation, genome assembly, and \ + strandedness setting for this library.", + merged.total_mapped, + shared_chroms.join(", ") + ); + } let bam_chroms: Vec<&str> = tid_to_name.iter().take(5).map(|s| s.as_str()).collect(); let gtf_chroms: Vec<&str> = index.keys().take(5).map(|s| s.as_str()).collect(); anyhow::bail!( diff --git a/src/rna/rseqc/accumulators.rs b/src/rna/rseqc/accumulators.rs index b91a409e..d49fae43 100644 --- a/src/rna/rseqc/accumulators.rs +++ b/src/rna/rseqc/accumulators.rs @@ -2336,9 +2336,9 @@ fn merge_vec_arrays(target: &mut Vec<[u64; N]>, source: Vec<[u64 if source.len() > target.len() { target.resize(source.len(), [0u64; N]); } - for (i, arr) in source.into_iter().enumerate() { - for j in 0..N { - target[i][j] += arr[j]; + for (target_row, source_row) in target.iter_mut().zip(source) { + for (target_count, source_count) in target_row.iter_mut().zip(source_row) { + *target_count += source_count; } } } diff --git a/src/rna/rseqc/stats.rs b/src/rna/rseqc/stats.rs index 20abf8c0..96d6444a 100644 --- a/src/rna/rseqc/stats.rs +++ b/src/rna/rseqc/stats.rs @@ -525,13 +525,13 @@ fn write_gcc( for cycle in 0..max_cycles { let mut combined = [0u64; 6]; if cycle < fbc.len() { - for j in 0..6 { - combined[j] += fbc[cycle][j]; + for (count, &value) in combined.iter_mut().zip(&fbc[cycle]) { + *count += value; } } if cycle < lbc.len() { - for j in 0..6 { - combined[j] += lbc[cycle][j]; + for (count, &value) in combined.iter_mut().zip(&lbc[cycle]) { + *count += value; } } // Denominator is A+C+G+T only (first 4 elements), matching upstream diff --git a/tests/integration_test.rs b/tests/integration_test.rs index c7ea3acd..2177d78e 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -1026,3 +1026,166 @@ fn test_dup_check_parallel_uses_global_duplicate_state() { let _ = fs::remove_dir_all(root); } + +// =================================================================== +// No-genic-reads diagnostic +// =================================================================== + +/// Run `rustqc rna` on one mapped read at `bam_chrom:11`, outside the only +/// gene (`gtf_chrom:900-3000`), with optional YAML config. Returns the output. +fn run_no_genic_reads( + label: &str, + bam_chrom: &str, + gtf_chrom: &str, + config: Option<&str>, +) -> std::process::Output { + let root = unique_test_dir(label); + let bam_path = root.join("s.bam"); + let gtf_path = root.join("a.gtf"); + let outdir = root.join("out"); + fs::create_dir_all(&outdir).unwrap(); + fs::write( + >f_path, + format!( + "{gtf_chrom}\tt\tgene\t900\t3000\t.\t+\t.\tgene_id \"g1\"; gene_name \"g1\";\n\ + {gtf_chrom}\tt\texon\t900\t3000\t.\t+\t.\tgene_id \"g1\"; transcript_id \"t1\"; exon_number \"1\";\n" + ), + ) + .unwrap(); + let read = format!( + "r0\t0\t{bam_chrom}\t11\t60\t50M\t*\t0\t0\tACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTAC\t??????????????????????????????????????????????????" + ); + write_bam_fixture(&bam_path, &[(bam_chrom, 100_000)], &[&read]); + + let mut cmd = Command::new(rustqc_binary()); + cmd.args(["rna", "--skip-dup-check", "--gtf"]) + .arg(>f_path) + .arg("--outdir") + .arg(&outdir) + .arg(&bam_path); + if let Some(yaml) = config { + let config_path = root.join("config.yaml"); + fs::write(&config_path, yaml).unwrap(); + cmd.arg("--config").arg(&config_path); + } + let output = cmd.output().expect("Failed to execute rustqc"); + let _ = fs::remove_dir_all(root); + output +} + +/// When BAM and GTF share chromosome names but no read lands in a gene, +/// the error must not claim a naming mismatch. It still fails. +#[test] +fn test_no_genic_reads_with_shared_names_is_not_a_name_mismatch() { + let output = run_no_genic_reads("no-genic-shared", "chr1", "chr1", None); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!( + !output.status.success(), + "rustqc should still fail:\n{stderr}" + ); + assert!( + !stderr.contains("Chromosome name mismatch"), + "BAM and GTF both name chr1, so this is not a name mismatch:\n{stderr}" + ); + assert!( + stderr.contains("No reads could be assigned to any gene"), + "missing no-assigned-reads diagnostic:\n{stderr}" + ); +} + +/// Names that only match after `chromosome_prefix` are not a mismatch either. +#[test] +fn test_no_genic_reads_with_prefix_mapped_names_is_not_a_name_mismatch() { + let output = run_no_genic_reads( + "no-genic-prefix", + "1", + "chr1", + Some("rna:\n chromosome_prefix: \"chr\"\n"), + ); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!( + !output.status.success(), + "rustqc should still fail:\n{stderr}" + ); + assert!( + !stderr.contains("Chromosome name mismatch"), + "chromosome_prefix maps 1 to chr1, so this is not a name mismatch:\n{stderr}" + ); + assert!( + stderr.contains("No reads could be assigned to any gene"), + "missing no-assigned-reads diagnostic:\n{stderr}" + ); +} + +/// A real naming mismatch (BAM `1`, GTF `chr1`) keeps the mismatch diagnostic. +#[test] +fn test_no_genic_reads_with_disjoint_names_is_a_name_mismatch() { + let output = run_no_genic_reads("no-genic-disjoint", "1", "chr1", None); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!(!output.status.success(), "rustqc should fail:\n{stderr}"); + assert!( + stderr.contains("Chromosome name mismatch"), + "missing chromosome-name diagnostic:\n{stderr}" + ); +} + +/// A physical gene overlap rejected by strandedness is an assignment failure. +#[test] +fn test_no_genic_reads_with_wrong_strand_is_assignment_failure() { + let root = unique_test_dir("no-genic-wrong-strand"); + let bam_path = root.join("s.bam"); + let gtf_path = root.join("a.gtf"); + fs::write( + >f_path, + "chr1\tt\tgene\t900\t3000\t.\t+\t.\tgene_id \"g1\"; gene_name \"g1\";\n\ + chr1\tt\texon\t900\t3000\t.\t+\t.\tgene_id \"g1\"; transcript_id \"t1\"; exon_number \"1\";\n", + ) + .unwrap(); + let read = format!( + "r0\t16\tchr1\t1001\t60\t50M\t*\t0\t0\t{}\t{}", + "A".repeat(50), + "I".repeat(50) + ); + write_bam_fixture(&bam_path, &[("chr1", 100_000)], &[&read]); + let run = |stranded: &str| { + Command::new(env!("CARGO_BIN_EXE_rustqc")) + .args(["rna", "--skip-dup-check", "--skip-preseq", "--gtf"]) + .arg(>f_path) + .args(["--stranded", stranded, "--outdir"]) + .arg(root.join(stranded)) + .arg(&bam_path) + .output() + .expect("Failed to execute rustqc") + }; + let control = run("unstranded"); + assert!( + control.status.success(), + "the overlapping read must be assigned in the unstranded control:\n{}", + String::from_utf8_lossy(&control.stderr) + ); + let output = run("forward"); + let stderr = String::from_utf8_lossy(&output.stderr); + assert!( + !output.status.success(), + "assignment failure must remain fatal:\n{stderr}" + ); + assert!( + !stderr.contains("Chromosome name mismatch"), + "chromosome names match:\n{stderr}" + ); + assert!( + !stderr.contains("overlaps a gene"), + "the read overlaps chr1:1001-1050 within g1 at chr1:900-3000; the unstranded control passed:\n{stderr}" + ); + assert!( + stderr.contains("configured gene-assignment rules"), + "missing assignment-rule diagnostic:\n{stderr}" + ); + assert!( + stderr.contains("annotation") + && stderr.contains("genome assembly") + && stderr.contains("strandedness"), + "missing diagnostic guidance:\n{stderr}" + ); + fs::remove_dir_all(root).unwrap(); +}