diff --git a/modkit-core/src/entropy/mod.rs b/modkit-core/src/entropy/mod.rs index 5de1ee4b..e88c0038 100644 --- a/modkit-core/src/entropy/mod.rs +++ b/modkit-core/src/entropy/mod.rs @@ -1386,8 +1386,13 @@ impl DescriptiveStats { "measurements and n_reads should be the same length" ); let mean_entropy = Self::mean(measurements); - let median_entropy = - percentile_linear_interp(measurements, 0.5f32)?; + let median_entropy = if measurements.len() == 1 { + measurements[0] + } else { + let mut sorted_measurements = measurements.to_vec(); + sorted_measurements.sort_unstable_by(f32::total_cmp); + percentile_linear_interp(&sorted_measurements, 0.5f32)? + }; // safe because of above check let (min_entropy, max_entropy) = match measurements.iter().minmax() { @@ -1671,7 +1676,118 @@ impl BedRegion { #[cfg(test)] mod entropy_mod_tests { - use crate::entropy::BedRegion; + use crate::entropy::{BedRegion, DescriptiveStats}; + use itertools::Itertools; + + #[test] + fn singleton_entropy_summary_uses_its_value_as_the_median() { + let stats = DescriptiveStats::new(&[0.25], &[7], 2, 0, &(10..11)) + .expect("a singleton is a valid regional entropy summary"); + + assert_eq!(stats.mean_entropy, 0.25); + assert_eq!(stats.median_entropy, 0.25); + assert_eq!(stats.min_entropy, 0.25); + assert_eq!(stats.max_entropy, 0.25); + assert_eq!(stats.mean_num_reads, 7.0); + assert_eq!(stats.min_num_reads, 7); + assert_eq!(stats.max_num_reads, 7); + assert_eq!(stats.successful_count, 1); + assert_eq!(stats.failed_count, 2); + } + + #[test] + fn multi_window_entropy_summary_keeps_existing_exact_statistics() { + let stats = DescriptiveStats::new( + &[0.25, 0.5, 0.75], + &[2, 4, 6], + 1, + 0, + &(10..20), + ) + .unwrap(); + + assert_eq!(stats.mean_entropy, 0.5); + assert_eq!(stats.median_entropy, 0.5); + assert_eq!(stats.min_entropy, 0.25); + assert_eq!(stats.max_entropy, 0.75); + assert_eq!(stats.mean_num_reads, 4.0); + assert_eq!(stats.min_num_reads, 2); + assert_eq!(stats.max_num_reads, 6); + assert_eq!(stats.successful_count, 3); + assert_eq!(stats.failed_count, 1); + } + + #[test] + fn region_median_is_independent_of_window_encounter_order() { + let observed = DescriptiveStats::new( + &[0.9, 0.1, 0.4, 0.2], + &[9, 1, 4, 2], + 0, + 0, + &(10..20), + ) + .unwrap(); + let sorted = DescriptiveStats::new( + &[0.1, 0.2, 0.4, 0.9], + &[1, 2, 4, 9], + 0, + 0, + &(10..20), + ) + .unwrap(); + + assert_eq!(observed.median_entropy, sorted.median_entropy); + assert_eq!(observed.median_entropy, 0.3); + } + + #[test] + fn odd_region_median_is_independent_of_window_encounter_order() { + let observed = DescriptiveStats::new( + &[0.9, 0.1, 0.4], + &[9, 1, 4], + 0, + 0, + &(10..20), + ) + .unwrap(); + let sorted = DescriptiveStats::new( + &[0.1, 0.4, 0.9], + &[1, 4, 9], + 0, + 0, + &(10..20), + ) + .unwrap(); + + assert_eq!(observed.median_entropy, sorted.median_entropy); + assert_eq!(observed.median_entropy, 0.4); + } + + #[test] + fn regional_medians_are_stable_across_all_small_encounter_orders() { + for (values, expected_median) in [ + (vec![0.75, 0.25, 0.5], 0.5), + (vec![0.75, 0.0, 0.25, 0.5], 0.375), + ] { + for measurements in + values.iter().copied().permutations(values.len()) + { + let original_order = measurements.clone(); + let reads = vec![1; measurements.len()]; + let stats = DescriptiveStats::new( + &measurements, + &reads, + 0, + 0, + &(10..20), + ) + .unwrap(); + + assert_eq!(stats.median_entropy, expected_median); + assert_eq!(measurements, original_order); + } + } + } #[test] fn test_bed_region_parsing() { diff --git a/modkit/tests/test_entropy_region_statistics.rs b/modkit/tests/test_entropy_region_statistics.rs new file mode 100644 index 00000000..bfd0f2a2 --- /dev/null +++ b/modkit/tests/test_entropy_region_statistics.rs @@ -0,0 +1,106 @@ +use rust_htslib::bam::{ + self, + header::HeaderRecord, + record::{Aux, Cigar, CigarString}, +}; +use std::fs; +use std::path::{Path, PathBuf}; +use std::process::Command; + +fn write_reference(root: &Path) -> PathBuf { + let reference = root.join("reference.fa"); + fs::write(&reference, ">chr1\nAACGAA\n").unwrap(); + fs::write(root.join("reference.fa.fai"), "chr1\t6\t6\t6\t7\n").unwrap(); + reference +} + +fn write_bam(root: &Path) -> PathBuf { + let bam_path = root.join("reads.bam"); + let mut header = bam::Header::new(); + let mut sq = HeaderRecord::new(b"SQ"); + sq.push_tag(b"SN", "chr1").push_tag(b"LN", 6); + header.push_record(&sq); + + let cigar = CigarString(vec![Cigar::Match(6)]); + let mut record = bam::Record::new(); + record.set(b"read-0", Some(&cigar), b"AACGAA", &[30; 6]); + record.set_tid(0); + record.set_pos(0); + record.set_mapq(60); + record.push_aux(b"MM", Aux::String("C+m?,0;")).unwrap(); + record.push_aux(b"ML", Aux::ArrayU8((&[255][..]).into())).unwrap(); + record.push_aux(b"MN", Aux::U32(6)).unwrap(); + record.push_aux(b"NM", Aux::U32(0)).unwrap(); + + let mut writer = + bam::Writer::from_path(&bam_path, &header, bam::Format::Bam).unwrap(); + writer.write(&record).unwrap(); + drop(writer); + bam::index::build(&bam_path, None, bam::index::Type::Bai, 1).unwrap(); + bam_path +} + +#[test] +fn singleton_region_emits_exact_summary_statistics() { + let temp_dir = tempfile::tempdir().unwrap(); + let reference = write_reference(temp_dir.path()); + let bam = write_bam(temp_dir.path()); + let regions = temp_dir.path().join("regions.bed"); + fs::write(®ions, "chr1\t2\t3\tsingleton\n").unwrap(); + let output_dir = temp_dir.path().join("output"); + + let result = Command::new(env!("CARGO_BIN_EXE_modkit")) + .args([ + "entropy", + "--in-bam", + bam.to_str().unwrap(), + "--out-bed", + output_dir.to_str().unwrap(), + "--ref", + reference.to_str().unwrap(), + "--base", + "C", + "--num-positions", + "1", + "--window-size", + "1", + "--min-coverage", + "1", + "--max-filtered-positions", + "0", + "--filter-threshold", + "0", + "--threads", + "1", + "--io-threads", + "1", + "--regions", + regions.to_str().unwrap(), + "--prefix", + "singleton", + "--suppress-progress", + ]) + .output() + .unwrap(); + + assert!( + result.status.success(), + "{}", + String::from_utf8_lossy(&result.stderr) + ); + let windows = fs::read_to_string( + output_dir.join("singleton_windows.bedgraph"), + ) + .unwrap(); + let window_fields = windows.trim_end().split('\t').collect::>(); + // Motif/window interval normalization belongs to issue #681. This test + // isolates the regional summary and only requires one successful window + // with the expected biological value and coverage. + assert_eq!(window_fields[0], "chr1"); + assert_eq!(window_fields[1], "2"); + assert_eq!(&window_fields[3..], &["0", "+", "1"]); + assert_eq!( + fs::read(output_dir.join("singleton_regions.bed")).unwrap(), + b"chr1\t2\t3\tsingleton\t0\t+\t0\t0\t0\t1\t1\t1\t1\t0\n" + ); +}