diff --git a/Cargo.lock b/Cargo.lock index c77b341b..ff396a4e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -375,9 +375,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", ] diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index 6ebc23a3..af86c8e2 100644 --- a/src/rna/bam/io.rs +++ b/src/rna/bam/io.rs @@ -282,7 +282,7 @@ enum BamCursor { chunk_end: Option, reference_sequence_id: usize, }, - /// Streaming unmapped records from a seeked position to end of file. + /// Streaming records with no reference from a seeked position to end of file. Unmapped, } @@ -436,8 +436,10 @@ impl IndexedReader { let source_path = self.source_path.clone(); match &mut self.backend { IndexedBackend::Bam { reader, cursor, .. } => { - // Seek to the start of the unmapped read region and stream from - // there, filtering for unmapped records as we go. + // Seek to the last linear-index window of the last reference and + // stream from there, keeping only records with no reference + // (htslib's HTS_IDX_NOCOOR). Placed unmapped mates in that window + // belong to their reference's fetch. match reader.index().last_first_record_start_position() { Some(pos) => { reader @@ -448,7 +450,7 @@ impl IndexedReader { None => { // The index carries no metadata pseudo-bin, so there is // no recorded unmapped offset. Reopen and stream from the - // first record, relying on the per-record unmapped filter. + // first record, relying on the per-record no-reference filter. let mut fresh = noodles_bam::io::indexed_reader::Builder::default() .build_from_path(&source_path) .with_context(|| { @@ -470,9 +472,25 @@ impl IndexedReader { reader: heads.reader, header: heads.header, iter_builder: |reader, header| -> Result { + // With no unplaced CRAI entry, noodles 0.94 seeks EOF + // and reports a truncated container instead of an empty query. + if reader + .index() + .iter() + .all(|entry| entry.reference_sequence_id().is_some()) + { + return Ok(Box::new(std::iter::empty())); + } + // The unplaced reads' slice can also hold placed records + // (multi-reference slice); keep only those with no reference. let query = reader .query_unmapped(header) - .context("failed to query unmapped CRAM reads")?; + .context("failed to query unmapped CRAM reads")? + .filter(|result| { + result.as_ref().map_or(true, |record| { + record.reference_sequence_id().is_none() + }) + }); Ok(Box::new(query)) }, } @@ -503,13 +521,13 @@ impl Read for IndexedReader { Ok(0) => return None, Ok(_) => { *cursor = BamCursor::Unmapped; - if scratch.flags().is_unmapped() { + if scratch.reference_sequence_id().is_none() { return Some( Record::from_bam(header.noodles_header(), scratch) .and_then(|loaded| store_record(record, loaded)), ); } - // Mapped record in the unmapped tail: skip and continue. + // Placed record before the unplaced tail: skip and continue. } Err(e) => return Some(Err(e.into())), }, @@ -518,22 +536,19 @@ impl Read for IndexedReader { chunk_end, reference_sequence_id, } => match chunk_end { - None => match chunks.next() { - Some(chunk) => { - if let Err(e) = - reader.get_mut().seek_to_virtual_position(chunk.start()) - { - return Some(Err(anyhow::Error::new(e) - .context("failed to seek to BAM index chunk"))); - } - *cursor = BamCursor::Region { - chunks, - chunk_end: Some(chunk.end()), - reference_sequence_id, - }; + None => { + let chunk = chunks.next()?; + if let Err(e) = reader.get_mut().seek_to_virtual_position(chunk.start()) + { + return Some(Err(anyhow::Error::new(e) + .context("failed to seek to BAM index chunk"))); } - None => return None, - }, + *cursor = BamCursor::Region { + chunks, + chunk_end: Some(chunk.end()), + reference_sequence_id, + }; + } Some(end) => { if reader.get_mut().virtual_position() >= end { // Reached the end of this chunk; advance to the next. @@ -613,7 +628,7 @@ mod tests { use noodles::sam::header::record::value::Map; use noodles::{cram, sam}; - use super::{Read, Reader}; + use super::{FetchTarget, IndexedReader, Read, Reader}; use crate::rna::bam::Record; static COUNTER: AtomicU64 = AtomicU64::new(0); @@ -693,4 +708,168 @@ mod tests { assert_eq!(got, expected.to_vec()); } + + // ------------------------------------------------------------------------ + // Unmapped sweep: htslib's HTS_IDX_NOCOOR returns exactly the records with + // no reference (tid -1). A placed unmapped mate (0x4 with a tid) belongs + // to its reference's fetch, so the sweep must not return it again. + // ------------------------------------------------------------------------ + + const CONTIG_LEN: usize = 100_000; + + /// Two mapped reads on chr1; on chr2, a mapped read whose unmapped mate is + /// placed at its position in chr2's last 16 kb linear-index window; two + /// unplaced reads (no reference) at the end. + fn sweep_world() -> (sam::Header, Vec) { + let header: sam::Header = format!( + "@HD\tVN:1.6\tSO:coordinate\n@SQ\tSN:chr1\tLN:{CONTIG_LEN}\n@SQ\tSN:chr2\tLN:{CONTIG_LEN}\n" + ) + .parse() + .unwrap(); + let rec = |name: &str, flags: u16, tid: Option, start: Option| { + let mut builder = RecordBuf::builder() + .set_name(name) + .set_flags(Flags::from(flags)) + .set_sequence(b"AAAA".to_vec().into()) + .set_quality_scores(vec![30u8; 4].into()); + if let (Some(tid), Some(start)) = (tid, start) { + let start = Position::new(start).unwrap(); + builder = builder + .set_reference_sequence_id(tid) + .set_alignment_start(start); + if flags & 0x1 != 0 { + builder = builder + .set_mate_reference_sequence_id(tid) + .set_mate_alignment_start(start); + } + if flags & 0x4 == 0 { + builder = builder + .set_mapping_quality( + sam::alignment::record::MappingQuality::new(60).unwrap(), + ) + .set_cigar([Op::new(Kind::Match, 4)].into_iter().collect()); + } + } + builder.build() + }; + let records = vec![ + rec("r0", 0x0, Some(0), Some(1_001)), + rec("r1", 0x0, Some(0), Some(1_201)), + rec("r2", 0x1 | 0x8 | 0x40, Some(1), Some(99_001)), + rec("r3", 0x1 | 0x4 | 0x80, Some(1), Some(99_001)), + rec("u1", 0x4, None, None), + rec("u2", 0x4, None, None), + ]; + (header, records) + } + + /// Fetch each reference, then the unmapped bucket; return (fetch, read name) + /// for every record returned. + fn fetch_everything( + path: &std::path::Path, + fasta: Option<&std::path::Path>, + ) -> Vec<(String, String)> { + let mut reader = IndexedReader::from_path(path).unwrap(); + if let Some(fasta) = fasta { + reader.set_reference(fasta).unwrap(); + } + let mut got = Vec::new(); + let mut record = Record::new(); + let targets = [ + ("tid 0", FetchTarget::Tid(0)), + ("tid 1", FetchTarget::Tid(1)), + ("unmapped", FetchTarget::Unmapped), + ]; + for (label, target) in targets { + reader.fetch(target).unwrap(); + while let Some(result) = reader.read(&mut record) { + result.unwrap(); + got.push(( + label.to_string(), + String::from_utf8(record.qname().to_vec()).unwrap(), + )); + } + } + got + } + + fn assert_sweep(got: &[(String, String)]) { + let swept: Vec<&str> = got + .iter() + .filter(|(label, _)| label == "unmapped") + .map(|(_, name)| name.as_str()) + .collect(); + assert_eq!( + swept, + ["u1", "u2"], + "unmapped sweep must return only the unplaced records; all fetches: {got:?}" + ); + } + + fn assert_each_record_once(got: &[(String, String)]) { + let mut names: Vec<&str> = got.iter().map(|(_, name)| name.as_str()).collect(); + names.sort_unstable(); + assert_eq!( + names, + ["r0", "r1", "r2", "r3", "u1", "u2"], + "each record must be fetched exactly once; all fetches: {got:?}" + ); + } + + #[test] + fn bam_unmapped_sweep_returns_only_unplaced_records() { + let (header, records) = sweep_world(); + let bam_path = temp_path("sweep.bam"); + { + let mut writer = noodles::bam::io::writer::Builder + .build_from_path(&bam_path) + .unwrap(); + writer.write_header(&header).unwrap(); + for record in &records { + writer.write_alignment_record(&header, record).unwrap(); + } + writer.try_finish().unwrap(); + } + let bai_path = std::path::PathBuf::from(format!("{}.bai", bam_path.display())); + let index = noodles::bam::fs::index(&bam_path).unwrap(); + noodles::bam::bai::fs::write(&bai_path, &index).unwrap(); + + let got = fetch_everything(&bam_path, None); + let _ = std::fs::remove_file(&bam_path); + let _ = std::fs::remove_file(&bai_path); + assert_sweep(&got); + assert_each_record_once(&got); + } + + /// samtools-written CRAMs of the same records (see the fixture README). + #[test] + fn cram_unmapped_sweep_returns_only_unplaced_records() { + let dir = std::path::Path::new("tests/data/unmapped_sweep"); + let fasta = dir.join("sweep.fa"); + + let got = fetch_everything(&dir.join("single_ref_slices.cram"), Some(&fasta)); + assert_sweep(&got); + assert_each_record_once(&got); + + // Only the sweep is checked here: a per-reference fetch on a + // multi-reference slice also returns the other references' records. + let got = fetch_everything(&dir.join("multi_ref_slice.cram"), Some(&fasta)); + assert_sweep(&got); + } + #[test] + fn cram_unmapped_sweep_with_no_unplaced_reads_is_empty() { + let dir = std::path::Path::new("tests/data/unmapped_sweep"); + let got = fetch_everything(&dir.join("placed_only.cram"), Some(&dir.join("sweep.fa"))); + let swept: Vec<_> = got + .iter() + .filter(|(target, _)| target == "unmapped") + .collect(); + assert!( + swept.is_empty(), + "no unplaced records should be returned: {got:?}" + ); + let mut names: Vec<_> = got.iter().map(|(_, name)| name.as_str()).collect(); + names.sort_unstable(); + assert_eq!(names, ["r0", "r1", "r2", "r3"]); + } } diff --git a/src/rna/rseqc/accumulators.rs b/src/rna/rseqc/accumulators.rs index a2ad5a13..6d0c4ad8 100644 --- a/src/rna/rseqc/accumulators.rs +++ b/src/rna/rseqc/accumulators.rs @@ -2326,9 +2326,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 a583e006..c8af662c 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/data/unmapped_sweep/README.md b/tests/data/unmapped_sweep/README.md new file mode 100644 index 00000000..f71d46e0 --- /dev/null +++ b/tests/data/unmapped_sweep/README.md @@ -0,0 +1,19 @@ +# Unmapped-sweep CRAM fixtures + +`sweep.sam` has two mapped reads on chr1, a mapped read on chr2 with its unmapped mate placed at the same position (flag 0x85, tid 1), and two unplaced reads (no reference). The CRAMs are that file written by samtools 1.22.1 against `sweep.fa`: + +```bash +samtools view -C -T sweep.fa -o single_ref_slices.cram sweep.sam +samtools view -C -T sweep.fa --output-fmt-option multi_seq_per_slice=1 -o multi_ref_slice.cram sweep.sam +samtools index single_ref_slices.cram +samtools index multi_ref_slice.cram +``` + +`single_ref_slices.cram` has one slice per reference plus one for the unplaced reads; `multi_ref_slice.cram` holds all six records in one multi-reference slice, so its index points the unplaced-read entry at a slice that also holds placed records. + +`placed_only.sam` drops the two unplaced reads from `sweep.sam`. It pins the empty-unplaced-sweep case; regenerate with samtools 1.22.1: + +```bash +samtools view --no-PG -C -T sweep.fa -o placed_only.cram placed_only.sam +samtools index placed_only.cram +``` diff --git a/tests/data/unmapped_sweep/multi_ref_slice.cram b/tests/data/unmapped_sweep/multi_ref_slice.cram new file mode 100644 index 00000000..9613ea84 Binary files /dev/null and b/tests/data/unmapped_sweep/multi_ref_slice.cram differ diff --git a/tests/data/unmapped_sweep/multi_ref_slice.cram.crai b/tests/data/unmapped_sweep/multi_ref_slice.cram.crai new file mode 100644 index 00000000..4eb5796b Binary files /dev/null and b/tests/data/unmapped_sweep/multi_ref_slice.cram.crai differ diff --git a/tests/data/unmapped_sweep/placed_only.cram b/tests/data/unmapped_sweep/placed_only.cram new file mode 100644 index 00000000..723129eb Binary files /dev/null and b/tests/data/unmapped_sweep/placed_only.cram differ diff --git a/tests/data/unmapped_sweep/placed_only.cram.crai b/tests/data/unmapped_sweep/placed_only.cram.crai new file mode 100644 index 00000000..6b2243b8 Binary files /dev/null and b/tests/data/unmapped_sweep/placed_only.cram.crai differ diff --git a/tests/data/unmapped_sweep/placed_only.sam b/tests/data/unmapped_sweep/placed_only.sam new file mode 100644 index 00000000..2a401fca --- /dev/null +++ b/tests/data/unmapped_sweep/placed_only.sam @@ -0,0 +1,7 @@ +@HD VN:1.6 SO:coordinate +@SQ SN:chr1 LN:1000 +@SQ SN:chr2 LN:1000 +r0 0 chr1 101 60 4M * 0 0 AAAA IIII +r1 0 chr1 201 60 4M * 0 0 AAAA IIII +r2 73 chr2 901 60 4M = 901 0 AAAA IIII +r3 133 chr2 901 0 * = 901 0 AAAA IIII diff --git a/tests/data/unmapped_sweep/single_ref_slices.cram b/tests/data/unmapped_sweep/single_ref_slices.cram new file mode 100644 index 00000000..d9341523 Binary files /dev/null and b/tests/data/unmapped_sweep/single_ref_slices.cram differ diff --git a/tests/data/unmapped_sweep/single_ref_slices.cram.crai b/tests/data/unmapped_sweep/single_ref_slices.cram.crai new file mode 100644 index 00000000..77f378f8 Binary files /dev/null and b/tests/data/unmapped_sweep/single_ref_slices.cram.crai differ diff --git a/tests/data/unmapped_sweep/sweep.fa b/tests/data/unmapped_sweep/sweep.fa new file mode 100644 index 00000000..a108fc96 --- /dev/null +++ b/tests/data/unmapped_sweep/sweep.fa @@ -0,0 +1,4 @@ +>chr1 +AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA +>chr2 +AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA diff --git a/tests/data/unmapped_sweep/sweep.sam b/tests/data/unmapped_sweep/sweep.sam new file mode 100644 index 00000000..2c86967e --- /dev/null +++ b/tests/data/unmapped_sweep/sweep.sam @@ -0,0 +1,9 @@ +@HD VN:1.6 SO:coordinate +@SQ SN:chr1 LN:1000 +@SQ SN:chr2 LN:1000 +r0 0 chr1 101 60 4M * 0 0 AAAA IIII +r1 0 chr1 201 60 4M * 0 0 AAAA IIII +r2 73 chr2 901 60 4M = 901 0 AAAA IIII +r3 133 chr2 901 0 * = 901 0 AAAA IIII +u1 4 * 0 0 * * 0 0 AAAA IIII +u2 4 * 0 0 * * 0 0 AAAA IIII diff --git a/tests/integration_test.rs b/tests/integration_test.rs index 45fc4509..14d7ef40 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -1026,3 +1026,124 @@ fn test_dup_check_parallel_uses_global_duplicate_state() { let _ = fs::remove_dir_all(root); } + +// =================================================================== +// Unmapped sweep on the parallel path +// =================================================================== + +/// A placed unmapped mate (0x4 with a reference) in the last reference's final +/// 16 kb, plus two unplaced reads. The parallel path's unmapped sweep must not +/// count the placed mate a second time. Expected text is `samtools flagstat` +/// 1.22.1 output for this input. +#[test] +fn test_parallel_unmapped_sweep_counts_placed_unmapped_mate_once() { + let expected = "\ +6 + 0 in total (QC-passed reads + QC-failed reads) +6 + 0 primary +0 + 0 secondary +0 + 0 supplementary +0 + 0 duplicates +0 + 0 primary duplicates +3 + 0 mapped (50.00% : N/A) +3 + 0 primary mapped (50.00% : N/A) +2 + 0 paired in sequencing +1 + 0 read1 +1 + 0 read2 +0 + 0 properly paired (0.00% : N/A) +0 + 0 with itself and mate mapped +1 + 0 singletons (50.00% : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + let seq = "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTAC"; + let qual = "??????????????????????????????????????????????????"; + let lines: Vec = [ + "r0\t0\tchr1\t1001\t60\t50M\t*\t0\t0", + "r1\t0\tchr1\t1201\t60\t50M\t*\t0\t0", + "r2\t73\tchr2\t99001\t60\t50M\t=\t99001\t0", + "r3\t133\tchr2\t99001\t0\t*\t=\t99001\t0", + "u1\t4\t*\t0\t0\t*\t*\t0\t0", + "u2\t4\t*\t0\t0\t*\t*\t0\t0", + ] + .iter() + .map(|fields| format!("{fields}\t{seq}\t{qual}")) + .collect(); + let lines: Vec<&str> = lines.iter().map(String::as_str).collect(); + + for threads in [1, 4] { + let root = unique_test_dir(&format!("unmapped-sweep-t{threads}")); + 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, + "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(); + write_bam_fixture(&bam_path, &[("chr1", 100_000), ("chr2", 100_000)], &lines); + + let output = run_rustqc_custom(&outdir, &bam_path, >f_path, threads, true, &[]); + assert!( + output.status.success(), + "rustqc failed with --threads {threads}:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + let flagstat = fs::read_to_string(outdir.join("samtools/s.flagstat")).unwrap(); + assert_eq!(flagstat, expected, "flagstat with --threads {threads}"); + let _ = fs::remove_dir_all(root); + } +} + +/// A CRAM with no unplaced records must succeed on the parallel sweep path too. +#[test] +fn test_parallel_cram_with_no_unplaced_reads() { + let expected = "\ +4 + 0 in total (QC-passed reads + QC-failed reads) +4 + 0 primary +0 + 0 secondary +0 + 0 supplementary +0 + 0 duplicates +0 + 0 primary duplicates +3 + 0 mapped (75.00% : N/A) +3 + 0 primary mapped (75.00% : N/A) +2 + 0 paired in sequencing +1 + 0 read1 +1 + 0 read2 +0 + 0 properly paired (0.00% : N/A) +0 + 0 with itself and mate mapped +1 + 0 singletons (50.00% : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + let data = Path::new("tests/data/unmapped_sweep"); + for threads in [1, 4] { + let root = unique_test_dir(&format!("cram-no-unplaced-t{threads}")); + let gtf_path = root.join("a.gtf"); + let outdir = root.join("out"); + fs::write(>f_path, + "chr1\tt\tgene\t1\t1000\t.\t+\t.\tgene_id \"g1\"; gene_name \"g1\";\n\ + chr1\tt\texon\t1\t1000\t.\t+\t.\tgene_id \"g1\"; transcript_id \"t1\"; exon_number \"1\";\n" + ).unwrap(); + let output = Command::new(env!("CARGO_BIN_EXE_rustqc")) + .args(["rna", "--skip-dup-check", "--threads", &threads.to_string()]) + .arg("--reference") + .arg(data.join("sweep.fa")) + .arg("--gtf") + .arg(>f_path) + .arg("--outdir") + .arg(&outdir) + .arg(data.join("placed_only.cram")) + .output() + .expect("Failed to execute rustqc"); + assert!( + output.status.success(), + "rustqc failed with --threads {threads}:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + let flagstat = fs::read_to_string(outdir.join("samtools/placed_only.flagstat")).unwrap(); + assert_eq!(flagstat, expected, "flagstat with --threads {threads}"); + fs::remove_dir_all(root).unwrap(); + } +}