diff --git a/.github/workflows/test_rust.yml b/.github/workflows/test_rust.yml index 15ed31703b..2829c2ffb8 100644 --- a/.github/workflows/test_rust.yml +++ b/.github/workflows/test_rust.yml @@ -18,7 +18,15 @@ jobs: with: persist-credentials: false - run: rustup update stable && rustup default stable - - name: cargo build + - name: cargo build (core only) run: cargo build --verbose - - name: cargo test + - name: cargo test (core only) run: cargo test --verbose + - name: cargo build (python extension) + run: cargo build --verbose --features python + - name: cargo test (python extension) + run: cargo test --verbose --features python + - name: core crate has no pyo3 dependency + run: | + cargo tree -e normal --prefix none > tree.txt + ! grep -q '^pyo3' tree.txt diff --git a/Cargo.lock b/Cargo.lock index 4feec2c49a..d5296e3327 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -249,6 +249,24 @@ version = "0.1.7" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ef8ae57c4978a2acd8b869ce6b9ca1dfe817bff704c220209fdef2c0b75a01b9" +[[package]] +name = "deeptools-core" +version = "4.0.0" +dependencies = [ + "bigtools", + "flate2", + "itertools 0.15.0", + "ndarray", + "ndarray-npy", + "openssl-sys", + "pyo3", + "rayon", + "rust-htslib", + "tempfile", + "thiserror 2.0.20", + "tokio", +] + [[package]] name = "derive-new" version = "0.6.0" @@ -473,23 +491,6 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" -[[package]] -name = "hp" -version = "0.1.0" -dependencies = [ - "bigtools", - "flate2", - "itertools 0.15.0", - "ndarray", - "ndarray-npy", - "openssl-sys", - "pyo3", - "rayon", - "rust-htslib", - "tempfile", - "tokio", -] - [[package]] name = "hts-sys" version = "2.2.1" diff --git a/Cargo.toml b/Cargo.toml index 9f8c063766..2c713e6cb9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,23 +1,36 @@ [package] -name = "hp" -version = "0.1.0" +name = "deeptools-core" +version = "4.0.0" edition = "2024" rust-version = "1.87" +description = "Pure-Rust core of deepTools: coverage normalization scale factors" +license-file = "LICENSE.txt" +repository = "https://github.com/deeptools/deepTools" +publish = false # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [lib] -name = "hp" -crate-type = ["cdylib"] +name = "deeptools_core" +crate-type = ["rlib", "cdylib"] + +[features] +python = [ + "dep:pyo3", "pyo3/extension-module", + "dep:openssl-sys", "dep:bigtools", "dep:rust-htslib", "dep:rayon", "dep:tokio", + "dep:itertools", "dep:flate2", "dep:tempfile", "dep:ndarray-npy", +] +default = [] [dependencies] -openssl-sys = { version = "0.9.117", features = ["vendored"] } -pyo3 = { version = ">= 0.28.2", features = ["extension-module"] } -bigtools = { version = "0.5.8", default-features = false, features = ["read", "write"] } -rust-htslib = "1.0.1" -rayon = "1.12.0" -tokio = "1.53.0" -itertools = "0.15.0" -flate2 = "1.1.9" -tempfile = "3.27.0" +openssl-sys = { version = "0.9.117", features = ["vendored"], optional = true } +pyo3 = { version = ">= 0.28.2", optional = true } +bigtools = { version = "0.5.8", default-features = false, features = ["read", "write"], optional = true } +rust-htslib = { version = "1.0.1", optional = true } +rayon = { version = "1.12.0", optional = true } +tokio = { version = "1.53.0", optional = true } +itertools = { version = "0.15.0", optional = true } +flate2 = { version = "1.1.9", optional = true } +tempfile = { version = "3.27.0", optional = true } ndarray = "0.17.2" -ndarray-npy = "0.10.0" +ndarray-npy = { version = "0.10.0", optional = true } +thiserror = "2.0" diff --git a/pyproject.toml b/pyproject.toml index 56d35d64d7..3cc5c29a81 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,6 +9,7 @@ python-source = "pydeeptools" exclude = ["pydeeptools/deeptools/test/**/*"] module-name = "deeptools.hp" bindings = "pyo3" +features = ["python"] [project] name = "deeptools" diff --git a/release-please-config.json b/release-please-config.json index 138caaa93a..56135f6334 100644 --- a/release-please-config.json +++ b/release-please-config.json @@ -11,6 +11,16 @@ { "type": "generic", "path": "galaxy/wrapper/deepTools_macros.xml" + }, + { + "type": "toml", + "path": "Cargo.toml", + "jsonpath": "$.package.version" + }, + { + "type": "toml", + "path": "Cargo.lock", + "jsonpath": "$.package[?(@.name.value=='deeptools-core')].version" } ] } diff --git a/scripts/check_version_sync.py b/scripts/check_version_sync.py index 5ff64a70d2..5818fc6b06 100644 --- a/scripts/check_version_sync.py +++ b/scripts/check_version_sync.py @@ -99,6 +99,20 @@ def galaxy_tool_version() -> str: return match.group(1) +def cargo_crate_version() -> str: + data = tomllib.loads((REPO_ROOT / "Cargo.toml").read_text()) + return data["package"]["version"] + + +def cargo_lock_crate_version() -> str: + name = tomllib.loads((REPO_ROOT / "Cargo.toml").read_text())["package"]["name"] + lock = tomllib.loads((REPO_ROOT / "Cargo.lock").read_text()) + for package in lock["package"]: + if package["name"] == name: + return package["version"] + sys.exit(f"could not find crate {name!r} in Cargo.lock") + + def main() -> None: cargo_version = cargo_rust_version() mismatches: dict[str, tuple[str, str]] = {} @@ -125,6 +139,13 @@ def main() -> None: galaxy_version, ) + for label, version in ( + ("deeptools: Cargo.toml [package] version", cargo_crate_version()), + ("deeptools: Cargo.lock crate version", cargo_lock_crate_version()), + ): + if version != pkg_version: + mismatches[label] = (pkg_version, version) + if mismatches: print("Version mismatch detected:", file=sys.stderr) for name, (expected, actual) in mismatches.items(): diff --git a/src/bamcompare.rs b/src/bamcompare.rs index 85dd8b8fa8..8293341a63 100644 --- a/src/bamcompare.rs +++ b/src/bamcompare.rs @@ -206,18 +206,24 @@ pub fn r_bamcompare( } (given_sf1, given_sf2) } else { + let to_py = |e: crate::normalization::NormalizationError| { + pyo3::exceptions::PyValueError::new_err(e.to_string()) + }; + let method: crate::normalization::ScaleFactorsMethod = + scalefactorsmethod.parse().map_err(to_py)?; + let norm_method: crate::normalization::Normalization = norm.parse().map_err(to_py)?; + let stats = |i: usize| crate::normalization::SampleStats { + mapped: covcalcs[i].mapped, + readlen: covcalcs[i].readlen, + }; let (sf1, sf2) = scale_factor_bamcompare( - scalefactorsmethod, - covcalcs[0].mapped, - covcalcs[1].mapped, + method, + [stats(0), stats(1)], binsize, effective_genome_size, - norm, - covcalcs[0].readlen, - covcalcs[1].readlen, - covcalcs[0].fraglen, - covcalcs[1].fraglen, - ); + norm_method, + ) + .map_err(to_py)?; if verbose { println!("scale factor1 = {}, scale factor2 = {}", sf1, sf2); } diff --git a/src/bamcoverage.rs b/src/bamcoverage.rs index 02654de209..9f186bf25b 100644 --- a/src/bamcoverage.rs +++ b/src/bamcoverage.rs @@ -66,10 +66,12 @@ pub fn r_bamcoverage( if mnase && !ispe { panic!("Error: MNase-seq requires paired-end data."); } - if norm == "RPGC" && effectivegenomesize == 0 { - panic!("Error: Effective genome size is required for RPGC normalization."); - } - if norm != "None" && scalefactor != 1.0 { + let norm_method: crate::normalization::Normalization = norm + .parse() + .map_err(|e: crate::normalization::NormalizationError| { + pyo3::exceptions::PyValueError::new_err(e.to_string()) + })?; + if norm_method != crate::normalization::Normalization::None && scalefactor != 1.0 { println!( "Warning: You have set a normalization option ({}), but also a scale factor. Only the scale factor will be used", norm @@ -230,16 +232,19 @@ pub fn r_bamcoverage( println!("Readlen: {}, Fraglen: {}", readlen, fraglen); } + let preset = if scalefactor != 1.0 { Some(scalefactor) } else { None }; let sf = scale_factor( - norm, + norm_method, mapped, binsize, effectivegenomesize, readlen, - fraglen, - scalefactor, - &verbose, - ); + preset, + ) + .map_err(|e| pyo3::exceptions::PyValueError::new_err(e.to_string()))?; + if verbose { + println!("Scale factor: {}", sf); + } // Create output stream let lines = bg.into_iter().flat_map(|bg| { diff --git a/src/calc.rs b/src/calc.rs index 0af033bd88..6a66bfbd1f 100644 --- a/src/calc.rs +++ b/src/calc.rs @@ -161,6 +161,15 @@ pub fn calc_ratio( } } +/// DESeq-style median-of-ratios scale factors. +/// +/// `array2` has one row per bin and one column per sample (raw counts). +/// Returns one scale factor per sample column: multiply a sample's counts by +/// its factor to make samples comparable. +/// +/// Entries `<= 0` are ignored, and so are bins where any sample is zero +/// (their log geometric mean is infinite). A sample with no usable bins gets +/// a factor of `1.0`. pub fn deseq_scalefactors(array2: &Array2) -> Array1 { let loggeomeans = array2.mapv(|v| v.ln()).sum_axis(Axis(1)) / array2.shape()[1] as f32; let masked_array = array2.mapv(|x| if x <= 0.0 { f32::NAN } else { x }); diff --git a/src/lib.rs b/src/lib.rs index 9a37b0398e..591e2ca4fe 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,17 +1,47 @@ +//! Core of deepTools: pure-Rust building blocks for coverage normalization. +//! +//! The crate has no dependency on `pyo3` unless the `python` feature is +//! enabled; that feature only exists to build the Python extension. +//! +//! # Public API +//! +//! * [`normalization`]: scale factors (`RPKM`, `CPM`, `BPM`, `RPGC`, read +//! count scaling between two samples). +//! * [`deseq_scalefactors`]: DESeq-style median-of-ratios scale factors. +//! +//! Everything else is an implementation detail of the Python extension and +//! may change without notice. +#![warn(missing_docs)] + +#[cfg(feature = "python")] use pyo3::prelude::*; -mod alignmentsieve; -mod bamcompare; -mod bamcoverage; +#[cfg_attr(not(feature = "python"), allow(dead_code))] mod calc; -mod computematrix; +#[cfg(feature = "python")] mod covcalc; +#[cfg(feature = "python")] mod filehandler; +#[cfg(feature = "python")] mod filtering; +pub mod normalization; + +pub use calc::deseq_scalefactors; + +#[cfg(feature = "python")] +mod alignmentsieve; +#[cfg(feature = "python")] +mod bamcompare; +#[cfg(feature = "python")] +mod bamcoverage; +#[cfg(feature = "python")] +mod computematrix; +#[cfg(feature = "python")] mod multibamsummary; -mod normalization; + #[cfg(test)] mod tests; +#[cfg(feature = "python")] #[pymodule] fn hp(_py: Python<'_>, m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(bamcoverage::r_bamcoverage, m)?)?; diff --git a/src/normalization.rs b/src/normalization.rs index 1fa0907389..d1eb636863 100644 --- a/src/normalization.rs +++ b/src/normalization.rs @@ -1,97 +1,175 @@ +//! Scale factor calculation for coverage normalization. +//! +//! All functions are pure: no I/O, no printing, no panics on bad input. Invalid +//! input is reported through [`NormalizationError`]. + use std::cmp; +use std::str::FromStr; + +/// Errors returned by the scale factor functions. +#[derive(Debug, Clone, PartialEq, thiserror::Error)] +pub enum NormalizationError { + /// A method name that is not one of the accepted values. + #[error("unknown {kind} '{value}', expected one of: {expected}")] + UnknownMethod { + /// What was being parsed (e.g. `normalization method`). + kind: &'static str, + /// The rejected input. + value: String, + /// Comma separated list of accepted values. + expected: &'static str, + }, + /// RPGC divides by the effective genome size, which must be non-zero. + #[error("RPGC normalization requires a non-zero effective genome size")] + MissingEffectiveGenomeSize, + /// A sample has no mapped reads, so its scale factor is undefined. + #[error("cannot compute a scale factor with zero mapped reads")] + ZeroMappedReads, +} + +/// Coverage normalization method (`--normalizeUsing`). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Normalization { + /// Reads per kilobase per million mapped reads. + Rpkm, + /// Counts per million mapped reads. + Cpm, + /// Bins per million mapped reads. + Bpm, + /// Reads per genomic content (1x coverage). + Rpgc, + /// No normalization; scale factor is 1. + None, +} + +impl FromStr for Normalization { + type Err = NormalizationError; + + /// Accepts `RPKM`, `CPM`, `BPM`, `RPGC` and `None` (case-sensitive). + fn from_str(s: &str) -> Result { + match s { + "RPKM" => Ok(Self::Rpkm), + "CPM" => Ok(Self::Cpm), + "BPM" => Ok(Self::Bpm), + "RPGC" => Ok(Self::Rpgc), + "None" => Ok(Self::None), + _ => Err(NormalizationError::UnknownMethod { + kind: "normalization method", + value: s.to_string(), + expected: "RPKM, CPM, BPM, RPGC, None", + }), + } + } +} + +/// How two samples are scaled against each other (`--scaleFactorsMethod`). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum ScaleFactorsMethod { + /// Scale the larger sample down to the read count of the smaller one. + ReadCount, + /// Scale each sample independently with a [`Normalization`] method. + None, +} + +impl FromStr for ScaleFactorsMethod { + type Err = NormalizationError; + /// Accepts `readCount` and `None` (case-sensitive). + fn from_str(s: &str) -> Result { + match s { + "readCount" => Ok(Self::ReadCount), + "None" => Ok(Self::None), + _ => Err(NormalizationError::UnknownMethod { + kind: "scale factors method", + value: s.to_string(), + expected: "readCount, None", + }), + } + } +} + +/// Scale factor for one sample. +/// +/// * `mapped`: number of mapped reads. +/// * `binsize`: bin size in bp. +/// * `effective_genome_size`: only used (and required non-zero) for RPGC. +/// * `readlen`: median read length in bp, only used for RPGC. +/// * `preset`: user supplied scale factor. When `Some`, it is returned as is +/// and `method` is ignored. pub fn scale_factor( - norm_method: &str, + method: Normalization, mapped: u32, binsize: u32, effective_genome_size: u64, readlen: f32, - _fraglen: f32, - scalefactor: f32, - verbose: &bool -) -> f32 { - if scalefactor != 1.0 { - return scalefactor; + preset: Option, +) -> Result { + if let Some(sf) = preset { + return Ok(sf); } - let mut scale_factor = 1.0; - - match norm_method { - "RPKM" => { - // RPKM = # reads per tile / total reads (millions) * tile length (kb) - let mmr = mapped as f32 / 1e6; - let bs_kb = binsize as f32 / 1000.0; - scale_factor *= 1.0 / (mmr * bs_kb); - } - "CPM" => { - // CPM = # reads per tile / total reads (millions) - let mmr = mapped as f32 / 1e6; - scale_factor *= 1.0 / mmr; - } - "BPM" => { - // BPM = bins per million mapped reads - let bs_kb: f32 = binsize as f32 / 1000.0; - let tmp_scalefactor = (mapped as f32 / bs_kb) / 1e6; - scale_factor *= 1.0 / (tmp_scalefactor * bs_kb); + if mapped == 0 && method != Normalization::None { + return Err(NormalizationError::ZeroMappedReads); + } + let mmr = mapped as f32 / 1e6; + let bs_kb = binsize as f32 / 1000.0; + Ok(match method { + // reads per tile / total reads (millions) / tile length (kb) + Normalization::Rpkm => 1.0 / (mmr * bs_kb), + // reads per tile / total reads (millions) + Normalization::Cpm => 1.0 / mmr, + // bins per million mapped reads + Normalization::Bpm => { + let tmp = (mapped as f32 / bs_kb) / 1e6; + 1.0 / (tmp * bs_kb) } - "RPGC" => { - // RPGC = mapped reads * fragment length / effective genome size - let tmp_scalefactor = (mapped as f32 * readlen as f32) / effective_genome_size as f32; - println!("Tmp scale factor: {}", tmp_scalefactor); - scale_factor *= 1.0 / tmp_scalefactor; + // 1x coverage: mapped reads * read length / effective genome size + Normalization::Rpgc => { + if effective_genome_size == 0 { + return Err(NormalizationError::MissingEffectiveGenomeSize); + } + 1.0 / ((mapped as f32 * readlen) / effective_genome_size as f32) } - _ => {} - }; - if *verbose { - println!("Scale factor: {}", scale_factor); - } - return scale_factor; + Normalization::None => 1.0, + }) } +/// Per-sample read statistics used by [`scale_factor_bamcompare`]. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct SampleStats { + /// Number of mapped reads. + pub mapped: u32, + /// Median read length in bp. + pub readlen: f32, +} + +/// Scale factors `(first, second)` for a pair of samples. +/// +/// With [`ScaleFactorsMethod::ReadCount`] the sample with more reads is scaled +/// down to the read count of the other. With [`ScaleFactorsMethod::None`] each +/// sample is scaled independently using `norm`. pub fn scale_factor_bamcompare( - norm_method: &str, - mapped_bam1: u32, - mapped_bam2: u32, + method: ScaleFactorsMethod, + samples: [SampleStats; 2], binsize: u32, effective_genome_size: u64, - norm: &str, - readlen_bam1: f32, - readlen_bam2: f32, - fraglen_bam1: f32, - fraglen_bam2: f32, -) -> (f32, f32) { - return match norm_method { - "readCount" => { - let min = cmp::min(mapped_bam1, mapped_bam2); - let scale_factor1 = min as f32 / mapped_bam1 as f32; - let scale_factor2 = min as f32 / mapped_bam2 as f32; - (scale_factor1, scale_factor2) - } - "None" => { - // Default to scale factor calculation with 'norm'. - let scale_factor1 = scale_factor( - norm, - mapped_bam1, - binsize, - effective_genome_size, - readlen_bam1, - fraglen_bam1, - 1.0, - &false - ); - let scale_factor2 = scale_factor( - norm, - mapped_bam2, - binsize, - effective_genome_size, - readlen_bam2, - fraglen_bam2, - 1.0, - &false - ); - (scale_factor1, scale_factor2) + norm: Normalization, +) -> Result<(f32, f32), NormalizationError> { + match method { + ScaleFactorsMethod::ReadCount => { + if samples[0].mapped == 0 || samples[1].mapped == 0 { + return Err(NormalizationError::ZeroMappedReads); + } + let min = cmp::min(samples[0].mapped, samples[1].mapped); + Ok(( + min as f32 / samples[0].mapped as f32, + min as f32 / samples[1].mapped as f32, + )) } - _ => { - panic!("ScaleFactorsMethod should either be 'readCount' or 'None'."); + ScaleFactorsMethod::None => { + let sf = |s: SampleStats| { + scale_factor(norm, s.mapped, binsize, effective_genome_size, s.readlen, None) + }; + Ok((sf(samples[0])?, sf(samples[1])?)) } } } diff --git a/src/tests/mod.rs b/src/tests/mod.rs index b27ecd53dd..a9d574a15b 100644 --- a/src/tests/mod.rs +++ b/src/tests/mod.rs @@ -1,5 +1,8 @@ pub mod test_calc; +#[cfg(feature = "python")] pub mod test_covcalc; +#[cfg(feature = "python")] pub mod test_filehandler; +#[cfg(feature = "python")] pub mod test_filtering; pub mod test_normalization; diff --git a/src/tests/test_calc.rs b/src/tests/test_calc.rs index 06fce4bb8a..8d47569686 100644 --- a/src/tests/test_calc.rs +++ b/src/tests/test_calc.rs @@ -89,4 +89,13 @@ mod scalefactor_calculations_tests { let expected_sf = Array1::from_vec(vec![1.233106, 0.9864848, 0.82207066]); assert_eq!(sf, expected_sf); } + + #[test] + fn test_deseq_scalefactors_zero_in_every_bin() { + // Every bin has a zero in column 0, so no bin has a usable geometric + // mean: nothing contributes and all factors fall back to 1. + let counts = Array2::from_shape_vec((2, 2), vec![0.0, 5.0, 0.0, 7.0]).unwrap(); + let sf = deseq_scalefactors(&counts); + assert_eq!(sf, Array1::from_vec(vec![1.0, 1.0])); + } } diff --git a/src/tests/test_normalization.rs b/src/tests/test_normalization.rs index d16a3f3c92..5b66d9f0a6 100644 --- a/src/tests/test_normalization.rs +++ b/src/tests/test_normalization.rs @@ -3,82 +3,161 @@ use crate::normalization::*; mod scale_factor_tests { use super::*; + fn sf(method: Normalization, preset: Option) -> Result { + scale_factor(method, 1000, 10, 10000, 20.0, preset) + } + #[test] fn test_scale_factor_rpkm() { - let nm = "RPKM"; - let sf = scale_factor( - nm, - 1000, 10, 10000, - 20.0, 20.0, - 1.0, &false - ); - assert_eq!(sf, 99999.99); + assert_eq!(sf(Normalization::Rpkm, None).unwrap(), 99999.99); } #[test] fn test_scale_factor_cpm() { - let nm = "CPM"; - let sf = scale_factor( - nm, - 1000, 10, 10000, - 20.0, 20.0, - 1.0, &false - ); - assert_eq!(sf, 999.99994); + assert_eq!(sf(Normalization::Cpm, None).unwrap(), 999.99994); } #[test] fn test_scale_factor_bpm() { - let nm = "BPM"; - let sf = scale_factor( - nm, - 1000, 10, 10000, - 20.0, 20.0, - 1.0, &false - ); - assert_eq!(sf, 999.99994); + assert_eq!(sf(Normalization::Bpm, None).unwrap(), 999.99994); } #[test] fn test_scale_factor_rpgc() { - let nm = "RPGC"; - let sf = scale_factor( - nm, - 1000, 10, 10000, - 20.0, 20.0, - 1.0, &false - ); - assert_eq!(sf, 0.5); + assert_eq!(sf(Normalization::Rpgc, None).unwrap(), 0.5); + } + + #[test] + fn test_scale_factor_none() { + assert_eq!(sf(Normalization::None, None).unwrap(), 1.0); } #[test] fn test_preset_scalefactor() { - let nm = "RPKM"; - let sf = scale_factor( - nm, - 1000, 10, 10000, - 20.0, 20.0, - 5.5, &false + assert_eq!(sf(Normalization::Rpkm, Some(5.5)).unwrap(), 5.5); + } + + #[test] + fn test_rpgc_needs_genome_size() { + assert_eq!( + scale_factor(Normalization::Rpgc, 1000, 10, 0, 20.0, None), + Err(NormalizationError::MissingEffectiveGenomeSize) + ); + } + + #[test] + fn test_preset_skips_rpgc_genome_size_check() { + assert_eq!( + scale_factor(Normalization::Rpgc, 1000, 10, 0, 20.0, Some(2.0)), + Ok(2.0) ); - assert_eq!(sf, 5.5); + } + + #[test] + fn test_zero_mapped_reads_is_an_error() { + for m in [Normalization::Rpkm, Normalization::Cpm, Normalization::Bpm, Normalization::Rpgc] { + assert_eq!( + scale_factor(m, 0, 10, 10000, 20.0, None), + Err(NormalizationError::ZeroMappedReads) + ); + } + } + + #[test] + fn test_zero_mapped_reads_ok_for_none_and_preset() { + assert_eq!(scale_factor(Normalization::None, 0, 10, 10000, 20.0, None), Ok(1.0)); + assert_eq!(scale_factor(Normalization::Cpm, 0, 10, 10000, 20.0, Some(3.0)), Ok(3.0)); } } -mod bamcompare_scale_factor_tests{ +mod parse_tests { use super::*; - + + #[test] + fn test_parse_normalization() { + assert_eq!("RPGC".parse::(), Ok(Normalization::Rpgc)); + assert_eq!("None".parse::(), Ok(Normalization::None)); + assert!(matches!( + "rpkm".parse::(), + Err(NormalizationError::UnknownMethod { .. }) + )); + } + + #[test] + fn test_parse_scale_factors_method() { + assert_eq!("readCount".parse::(), Ok(ScaleFactorsMethod::ReadCount)); + assert!(matches!( + "bogus".parse::(), + Err(NormalizationError::UnknownMethod { .. }) + )); + } +} + +mod bamcompare_scale_factor_tests { + use super::*; + + fn samples() -> [SampleStats; 2] { + [ + SampleStats { mapped: 1000, readlen: 20.0 }, + SampleStats { mapped: 2000, readlen: 20.0 }, + ] + } + #[test] fn test_bamcompare_readcount() { - let nm = "readCount"; let (sf1, sf2) = scale_factor_bamcompare( - nm, - 1000, 2000, - 10, 10000, - "RPKM", - 20.0, 20.0, - 1.0, 1.0 - ); + ScaleFactorsMethod::ReadCount, + samples(), + 10, + 10000, + Normalization::Rpkm, + ) + .unwrap(); assert_eq!(sf1, 1.0); assert_eq!(sf2, 0.5); } -} \ No newline at end of file + + #[test] + fn test_bamcompare_independent_norm() { + let (sf1, sf2) = scale_factor_bamcompare( + ScaleFactorsMethod::None, + samples(), + 10, + 10000, + Normalization::Cpm, + ) + .unwrap(); + assert_eq!(sf1, 999.99994); + assert!((sf2 - 500.0).abs() < 1e-3, "sf2 = {sf2}"); + } + + #[test] + fn test_bamcompare_rpgc_needs_genome_size() { + assert_eq!( + scale_factor_bamcompare( + ScaleFactorsMethod::None, + samples(), + 10, + 0, + Normalization::Rpgc, + ), + Err(NormalizationError::MissingEffectiveGenomeSize) + ); + } +} + +mod bamcompare_zero_reads_tests { + use super::*; + + #[test] + fn test_readcount_with_an_empty_sample_is_an_error() { + let samples = [ + SampleStats { mapped: 0, readlen: 20.0 }, + SampleStats { mapped: 2000, readlen: 20.0 }, + ]; + assert_eq!( + scale_factor_bamcompare(ScaleFactorsMethod::ReadCount, samples, 10, 10000, Normalization::None), + Err(NormalizationError::ZeroMappedReads) + ); + } +}