diff --git a/modkit-core/src/bedmethyl_util/subcommands.rs b/modkit-core/src/bedmethyl_util/subcommands.rs index 96cf1b5..5059691 100644 --- a/modkit-core/src/bedmethyl_util/subcommands.rs +++ b/modkit-core/src/bedmethyl_util/subcommands.rs @@ -182,25 +182,62 @@ pub struct EntryMergeBedMethyl { min_sample_coverage: Option, } -type BedMethylChunk = Vec; +#[derive(Debug)] +struct MergeBedMethylLine { + record: BedMethylLine, + name: String, +} + +impl ParseBedLine for MergeBedMethylLine { + fn parse(line: &str) -> Result { + let record = BedMethylLine::parse(line)?; + let name = line + .split_whitespace() + .nth(3) + .ok_or_else(|| { + MkError::InvalidBedMethyl(format!( + "missing name field in bedMethyl record:\n{line}" + )) + })? + .to_string(); + Ok(Self { record, name }) + } + + fn overlaps(&self, strand_rule: StrandRule) -> bool { + self.record.strand.overlaps(&strand_rule) + } + + fn to_line(&self) -> String { + let parsed_line = self.record.to_line(); + let mut fields = parsed_line.splitn(5, '\t'); + let chrom = fields.next().unwrap(); + let start = fields.next().unwrap(); + let stop = fields.next().unwrap(); + let _parsed_name = fields.next().unwrap(); + let remaining_fields = fields.next().unwrap(); + format!("{chrom}\t{start}\t{stop}\t{}\t{remaining_fields}", self.name) + } +} + +type BedMethylChunk = Vec; fn merge_data( - readers: &[HtsTabixHandler], + readers: &[HtsTabixHandler], chrom_coordinates: ChromCoordinates, tid_to_name: &FxHashMap, io_threads: usize, min_samples: usize, min_sample_coverage: u64, ) -> anyhow::Result { - type Key = (u64, ModCodeRepr, StrandRule); + type Key = (u64, ModCodeRepr, StrandRule, String); // this is safe because of how we constructed this let contig = tid_to_name.get(&chrom_coordinates.chrom_tid).unwrap(); let range = (chrom_coordinates.start_pos as u64) ..(chrom_coordinates.end_pos as u64); - // value is the merged record plus a tally of how many inputs contributed to - // it (each input has at most one record per key), used for the - // --min-samples (inner-join) filter below. - let mut merged_data = FxHashMap::::default(); + // value is the merged record plus a tally of how many requested inputs + // contributed to it, used for the --min-samples filter below. + let mut merged_data = + FxHashMap::::default(); // rationale: // iterate over every possible contig @@ -210,30 +247,45 @@ fn merge_data( // lines in a hashmap write the hashmap to a new bedmethyl // recreate hashmap and repeat process for next contig/regions for index in readers.iter() { - let lines = index.read_bedmethyl(&contig, &range, io_threads)?; + let lines = index.fetch_region( + &contig, + &range, + StrandRule::Both, + io_threads, + )?; + let mut contributed_keys = FxHashSet::::default(); for line in lines { - let line = line?; - // an input only contributes to a position when its record has at // least the requested valid coverage - if line.valid_coverage < min_sample_coverage { + if line.record.valid_coverage < min_sample_coverage { continue; } + let key = ( + line.record.start(), + line.record.raw_mod_code, + line.record.strand, + line.name.clone(), + ); + let first_for_input = contributed_keys.insert(key.clone()); merged_data - .entry((line.start(), line.raw_mod_code, line.strand)) + .entry(key) // modify the methyl data if an entry is found .and_modify(|(methyl, n_samples)| { - methyl.count_methylated += line.count_methylated; - methyl.valid_coverage += line.valid_coverage; - methyl.count_canonical += line.count_canonical; - methyl.count_other += line.count_other; - methyl.count_delete += line.count_delete; - methyl.count_fail += line.count_fail; - methyl.count_diff += line.count_diff; - methyl.count_nocall += line.count_nocall; - *n_samples += 1; + methyl.record.count_methylated += + line.record.count_methylated; + methyl.record.valid_coverage += line.record.valid_coverage; + methyl.record.count_canonical += + line.record.count_canonical; + methyl.record.count_other += line.record.count_other; + methyl.record.count_delete += line.record.count_delete; + methyl.record.count_fail += line.record.count_fail; + methyl.record.count_diff += line.record.count_diff; + methyl.record.count_nocall += line.record.count_nocall; + if first_for_input { + *n_samples += 1; + } }) .or_insert((line, 1)); } @@ -246,12 +298,21 @@ fn merge_data( .filter(|(_, n_samples)| *n_samples >= min_samples) .map(|(methyl, _)| methyl) .sorted_by(|a, b| { - debug_assert_eq!(a.chrom, b.chrom); - match a.start().cmp(&b.start()) { - Ordering::Equal => match a.strand.cmp(&b.strand) { - Ordering::Equal => a.raw_mod_code.cmp(&b.raw_mod_code), - o @ _ => o, - }, + debug_assert_eq!(a.record.chrom, b.record.chrom); + match a.record.start().cmp(&b.record.start()) { + Ordering::Equal => { + match a.record.strand.cmp(&b.record.strand) { + Ordering::Equal => match a + .record + .raw_mod_code + .cmp(&b.record.raw_mod_code) + { + Ordering::Equal => a.name.cmp(&b.name), + o @ _ => o, + }, + o @ _ => o, + } + } o @ _ => o, } }) @@ -282,6 +343,30 @@ impl EntryMergeBedMethyl { } let min_sample_coverage: u64 = self.min_sample_coverage.unwrap_or(0); + // Validate every requested input before opening the output. In + // particular, do not silently shrink the reader list: --min-samples + // is resolved against the requested input cardinality above. + let readers = self + .in_bedmethyl + .iter() + .map(|bedmethyl| { + File::open(bedmethyl).with_context(|| { + format!( + "failed to open input bedMethyl file {}", + bedmethyl.display() + ) + })?; + + HtsTabixHandler::::from_path(bedmethyl) + .with_context(|| { + format!( + "failed to read indexed input bedMethyl file {}", + bedmethyl.display() + ) + }) + }) + .collect::>>()?; + let pool = rayon::ThreadPoolBuilder::new() .num_threads(self.threads) .build()?; @@ -309,22 +394,6 @@ impl EntryMergeBedMethyl { writer.write(bedmethyl_header().as_bytes())?; } - let readers = self - .in_bedmethyl - .iter() - .filter_map(|bedmethyl| { - let index: HtsTabixHandler = - match HtsTabixHandler::from_path(&bedmethyl) { - Ok(reader) => reader, - Err(_) => { - return None; - } - }; - - Some(index) - }) - .collect::>>(); - // get set of contigs from all files // done this way in case one file has a set of contigs that the other // bedmethyl files do not have diff --git a/modkit/tests/test_bedmethyl_util.rs b/modkit/tests/test_bedmethyl_util.rs index 63f64a9..653f9e6 100644 --- a/modkit/tests/test_bedmethyl_util.rs +++ b/modkit/tests/test_bedmethyl_util.rs @@ -1,6 +1,9 @@ use std::{ - fs::File, + ffi::CString, + fs::{self, File}, io::{BufRead, BufReader, BufWriter, Read, Write}, + path::Path, + process::{Command, Output}, }; use common::run_modkit; @@ -54,6 +57,26 @@ fn test_bedmethyl_merge() { .lines() .map(|line| BedMethylLine::parse(&line.unwrap()).unwrap()) .collect::>(); + let expected_output = input_records + .iter() + .cloned() + .map(|mut line| { + line.count_methylated *= 2; + line.valid_coverage *= 2; + line.count_canonical *= 2; + line.count_other *= 2; + line.count_delete *= 2; + line.count_fail *= 2; + line.count_diff *= 2; + line.count_nocall *= 2; + line.to_string() + }) + .collect::(); + assert_eq!( + fs::read(&out_bed).unwrap(), + expected_output.as_bytes(), + "valid merge output must remain byte-identical" + ); let merged_records = BufReader::new(File::open(out_bed).unwrap()) .lines() .map(|line| BedMethylLine::parse(&line.unwrap()).unwrap()) @@ -90,10 +113,160 @@ fn min_samples_test_sizes() -> std::path::PathBuf { sizes_fp } -const MIN_SAMPLES_BED_FP: &str = - "../tests/resources/\ +const MIN_SAMPLES_BED_FP: &str = "../tests/resources/\ lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz"; +const PREFLIGHT_OUTPUT_SENTINEL: &[u8] = b"existing output must survive\n"; + +fn copy_bedmethyl_fixture(data_dst: &Path) { + fs::copy(MIN_SAMPLES_BED_FP, data_dst).unwrap(); +} + +fn run_merge_for_preflight( + inputs: &[&Path], + sizes_fp: &Path, + out_fp: &Path, + extra_args: &[&str], +) -> Output { + let mut command = Command::new(env!("CARGO_BIN_EXE_modkit")); + command.args(["bedmethyl", "merge"]); + for input in inputs { + command.arg(input); + } + command + .arg("-g") + .arg(sizes_fp) + .arg("-o") + .arg(out_fp) + .arg("--force") + .args(extra_args) + .output() + .unwrap() +} + +fn assert_preflight_failure(output: Output, bad_input: &Path, out_fp: &Path) { + let stderr = String::from_utf8_lossy(&output.stderr); + assert!( + !output.status.success(), + "merge unexpectedly accepted invalid input {}\nstderr: {stderr}", + bad_input.display() + ); + let bad_name = bad_input.file_name().unwrap().to_string_lossy(); + assert!( + stderr.contains(bad_name.as_ref()), + "error must identify invalid input {}\nstderr: {stderr}", + bad_input.display() + ); + assert_eq!( + fs::read(out_fp).unwrap(), + PREFLIGHT_OUTPUT_SENTINEL, + "input preflight failure must precede output truncation" + ); +} + +fn preflight_case_paths( + tempdir: &tempfile::TempDir, +) -> (std::path::PathBuf, std::path::PathBuf) { + let sizes_fp = tempdir.path().join("sizes.tsv"); + fs::write(&sizes_fp, b"chr20\t64444167\n").unwrap(); + let out_fp = tempdir.path().join("merged.bed"); + fs::write(&out_fp, PREFLIGHT_OUTPUT_SENTINEL).unwrap(); + (sizes_fp, out_fp) +} + +#[test] +fn test_bedmethyl_merge_rejects_missing_input_before_output() { + let tempdir = tempfile::tempdir().unwrap(); + let (sizes_fp, out_fp) = preflight_case_paths(&tempdir); + let missing = tempdir.path().join("missing.bed.gz"); + let output = run_merge_for_preflight( + &[Path::new(MIN_SAMPLES_BED_FP), &missing], + &sizes_fp, + &out_fp, + &[], + ); + assert_preflight_failure(output, &missing, &out_fp); +} + +#[test] +fn test_bedmethyl_merge_rejects_absent_index_before_output() { + let tempdir = tempfile::tempdir().unwrap(); + let (sizes_fp, out_fp) = preflight_case_paths(&tempdir); + let unindexed = tempdir.path().join("unindexed.bed.gz"); + copy_bedmethyl_fixture(&unindexed); + let output = run_merge_for_preflight( + &[Path::new(MIN_SAMPLES_BED_FP), &unindexed], + &sizes_fp, + &out_fp, + &[], + ); + assert_preflight_failure(output, &unindexed, &out_fp); +} + +#[test] +fn test_bedmethyl_merge_rejects_corrupt_index_before_output() { + let tempdir = tempfile::tempdir().unwrap(); + let (sizes_fp, out_fp) = preflight_case_paths(&tempdir); + let corrupt = tempdir.path().join("corrupt-index.bed.gz"); + copy_bedmethyl_fixture(&corrupt); + fs::write(format!("{}.tbi", corrupt.display()), b"not a tabix index") + .unwrap(); + let output = run_merge_for_preflight( + &[Path::new(MIN_SAMPLES_BED_FP), &corrupt], + &sizes_fp, + &out_fp, + &[], + ); + assert_preflight_failure(output, &corrupt, &out_fp); +} + +#[test] +fn test_bedmethyl_merge_rejects_one_bad_input_with_min_samples_all() { + let tempdir = tempfile::tempdir().unwrap(); + let (sizes_fp, out_fp) = preflight_case_paths(&tempdir); + let missing = tempdir.path().join("missing-middle.bed.gz"); + let valid = Path::new(MIN_SAMPLES_BED_FP); + let output = run_merge_for_preflight( + &[valid, &missing, valid], + &sizes_fp, + &out_fp, + &["--min-samples", "all"], + ); + assert_preflight_failure(output, &missing, &out_fp); +} + +#[test] +fn test_bedmethyl_merge_accepts_csi_index() { + let tempdir = tempfile::tempdir().unwrap(); + let (sizes_fp, out_fp) = preflight_case_paths(&tempdir); + let csi_input = tempdir.path().join("csi-indexed.bed.gz"); + copy_bedmethyl_fixture(&csi_input); + + let csi_input_cstr = CString::new(csi_input.to_str().unwrap()).unwrap(); + let result = unsafe { + rust_htslib::htslib::tbx_index_build( + csi_input_cstr.as_ptr(), + 14, + &rust_htslib::htslib::tbx_conf_bed, + ) + }; + assert_eq!(result, 0, "failed to build CSI fixture"); + assert!(Path::new(&format!("{}.csi", csi_input.display())).exists()); + + let output = run_merge_for_preflight( + &[&csi_input, &csi_input], + &sizes_fp, + &out_fp, + &[], + ); + assert!( + output.status.success(), + "merge rejected a valid CSI-backed input: {}", + String::from_utf8_lossy(&output.stderr) + ); + assert!(count_bed_records(&out_fp) > 0); +} + fn count_bed_records(fp: &std::path::Path) -> usize { BufReader::new(File::open(fp).unwrap()) .lines() @@ -196,3 +369,99 @@ fn test_bedmethyl_merge_min_sample_coverage() { .unwrap(); assert_eq!(count_bed_records(&out_bed), 0); } + +fn write_indexed_bedmethyl(path: &Path, rows: &[&str]) { + let mut writer = rust_htslib::bgzf::Writer::from_path(path).unwrap(); + for row in rows { + writeln!(writer, "{row}").unwrap(); + } + drop(writer); + + let c_path = CString::new(path.to_str().unwrap()).unwrap(); + let result = unsafe { + rust_htslib::htslib::tbx_index_build( + c_path.as_ptr(), + 0, + &rust_htslib::htslib::tbx_conf_bed, + ) + }; + assert_eq!(result, 0, "failed to index {}", path.display()); +} + +#[test] +fn test_bedmethyl_merge_preserves_full_name_and_counts_distinct_inputs() { + let temp_dir = tempfile::tempdir().unwrap(); + let root = temp_dir.path(); + let sizes = root.join("sizes.tsv"); + std::fs::write(&sizes, "chr1\t100\n").unwrap(); + + // Repeating a requested path remains a distinct sample, as in the + // existing merge CLI. These fixtures instead duplicate records within + // one input, which must not inflate that input's --min-samples tally. + let input_a_rows = [ + "chr1\t10\t11\tm,CG,0\t4\t+\t10\t11\t255,0,0\t4\t25.00\t1\t2\t1\t1\t1\t1\t1", + "chr1\t10\t11\tm,DRACH,2\t5\t+\t10\t11\t255,0,0\t5\t40.00\t2\t2\t1\t2\t0\t1\t0", + "chr1\t10\t11\tm,CG,0\t6\t+\t10\t11\t255,0,0\t6\t33.33\t2\t3\t1\t2\t2\t2\t2", + "chr1\t10\t11\th,CG,0\t3\t+\t10\t11\t255,0,0\t3\t33.33\t1\t2\t0\t3\t3\t3\t3", + "chr1\t10\t11\th,CG,0\t4\t+\t10\t11\t255,0,0\t4\t50.00\t2\t2\t0\t4\t4\t4\t4", + "chr1\t10\t11\ta\t2\t+\t10\t11\t255,0,0\t2\t50.00\t1\t1\t0\t1\t1\t1\t1", + "chr1\t10\t11\t76792,CG,0\t1\t+\t10\t11\t255,0,0\t1\t100.00\t1\t0\t0\t1\t1\t1\t1", + ]; + let input_b_rows = [ + "chr1\t10\t11\tm,CG,0\t5\t+\t10\t11\t255,0,0\t5\t80.00\t4\t1\t0\t4\t4\t4\t4", + "chr1\t10\t11\tm,DRACH,2\t5\t+\t10\t11\t255,0,0\t5\t60.00\t3\t2\t0\t3\t3\t3\t3", + "chr1\t10\t11\ta\t3\t+\t10\t11\t255,0,0\t3\t66.67\t2\t1\t0\t2\t2\t2\t2", + "chr1\t10\t11\t76792,CG,0\t2\t+\t10\t11\t255,0,0\t2\t50.00\t1\t1\t0\t2\t2\t2\t2", + ]; + + let input_a = root.join("a.bed.gz"); + let input_b = root.join("b.bed.gz"); + let input_a_reversed = root.join("a-reversed.bed.gz"); + let input_b_reversed = root.join("b-reversed.bed.gz"); + write_indexed_bedmethyl(&input_a, &input_a_rows); + write_indexed_bedmethyl(&input_b, &input_b_rows); + write_indexed_bedmethyl( + &input_a_reversed, + &input_a_rows.iter().copied().rev().collect::>(), + ); + write_indexed_bedmethyl( + &input_b_reversed, + &input_b_rows.iter().copied().rev().collect::>(), + ); + + let expected = concat!( + "chr1\t10\t11\ta\t5\t+\t10\t11\t255,0,0\t5\t60.00\t3\t2\t0\t3\t3\t3\t3\n", + "chr1\t10\t11\tm,CG,0\t15\t+\t10\t11\t255,0,0\t15\t46.67\t7\t6\t2\t7\t7\t7\t7\n", + "chr1\t10\t11\tm,DRACH,2\t10\t+\t10\t11\t255,0,0\t10\t50.00\t5\t4\t1\t5\t3\t4\t3\n", + "chr1\t10\t11\t76792,CG,0\t3\t+\t10\t11\t255,0,0\t3\t66.67\t2\t1\t0\t3\t3\t3\t3\n", + ); + + for (run, inputs) in + [[&input_a, &input_b], [&input_b_reversed, &input_a_reversed]] + .into_iter() + .enumerate() + { + let output = root.join(format!("merged-{run}.bed")); + run_modkit(&[ + "bedmethyl", + "merge", + inputs[0].to_str().unwrap(), + inputs[1].to_str().unwrap(), + "--genome-sizes", + sizes.to_str().unwrap(), + "--out-bed", + output.to_str().unwrap(), + "--threads", + "1", + "--io-threads", + "1", + "--min-samples", + "all", + ]) + .unwrap(); + + let observed = std::fs::read_to_string(output).unwrap(); + assert_eq!(observed, expected); + assert!(!observed.contains("h,CG,0")); + } +}