diff --git a/.github/workflows/test_rust.yml b/.github/workflows/test_rust.yml index 15ed31703b..325bb13283 100644 --- a/.github/workflows/test_rust.yml +++ b/.github/workflows/test_rust.yml @@ -22,3 +22,7 @@ jobs: run: cargo build --verbose - name: cargo test run: cargo test --verbose + - name: cargo build (core only) + run: cargo build --verbose --no-default-features + - name: cargo test (core only) + run: cargo test --verbose --no-default-features diff --git a/Cargo.lock b/Cargo.lock index 4feec2c49a..e7e8c0f827 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -249,6 +249,23 @@ 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", + "tokio", +] + [[package]] name = "derive-new" version = "0.6.0" @@ -473,23 +490,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..7ade608558 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,17 +1,21 @@ [package] -name = "hp" -version = "0.1.0" +name = "deeptools-core" +version = "4.0.0" edition = "2024" rust-version = "1.87" # 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"] +default = ["python"] [dependencies] openssl-sys = { version = "0.9.117", features = ["vendored"] } -pyo3 = { version = ">= 0.28.2", features = ["extension-module"] } +pyo3 = { version = ">= 0.28.2", optional = true } bigtools = { version = "0.5.8", default-features = false, features = ["read", "write"] } rust-htslib = "1.0.1" rayon = "1.12.0" diff --git a/galaxy/wrapper/test-data/alignmentSieve2.bam b/galaxy/wrapper/test-data/alignmentSieve2.bam index aebd0399f4..a9ebc760c8 100644 Binary files a/galaxy/wrapper/test-data/alignmentSieve2.bam and b/galaxy/wrapper/test-data/alignmentSieve2.bam differ diff --git a/galaxy/wrapper/test-data/alignmentSieve3.bam b/galaxy/wrapper/test-data/alignmentSieve3.bam index 7dd134440a..cb3448c99b 100644 Binary files a/galaxy/wrapper/test-data/alignmentSieve3.bam and b/galaxy/wrapper/test-data/alignmentSieve3.bam differ diff --git a/pydeeptools/deeptools/test/test_data/alignmentSieve4.bam b/pydeeptools/deeptools/test/test_data/alignmentSieve4.bam index 7cb306a71c..11326f9664 100644 Binary files a/pydeeptools/deeptools/test/test_data/alignmentSieve4.bam and b/pydeeptools/deeptools/test/test_data/alignmentSieve4.bam differ diff --git a/pydeeptools/deeptools/test/test_data/alignmentSieve5.bam b/pydeeptools/deeptools/test/test_data/alignmentSieve5.bam index 72c1aff9d5..9bea0a0e76 100644 Binary files a/pydeeptools/deeptools/test/test_data/alignmentSieve5.bam and b/pydeeptools/deeptools/test/test_data/alignmentSieve5.bam differ diff --git a/pydeeptools/deeptools/test/test_data/alignmentSieve6.bam b/pydeeptools/deeptools/test/test_data/alignmentSieve6.bam index c56bbb9ece..28fdb77d74 100644 Binary files a/pydeeptools/deeptools/test/test_data/alignmentSieve6.bam and b/pydeeptools/deeptools/test/test_data/alignmentSieve6.bam differ diff --git a/pydeeptools/deeptools/test/test_data/alignmentSieve7.bam b/pydeeptools/deeptools/test/test_data/alignmentSieve7.bam index ec948f28c9..231389b804 100644 Binary files a/pydeeptools/deeptools/test/test_data/alignmentSieve7.bam and b/pydeeptools/deeptools/test/test_data/alignmentSieve7.bam differ diff --git a/pydeeptools/deeptools/test/test_data/alignmentSieve8.bam b/pydeeptools/deeptools/test/test_data/alignmentSieve8.bam index 66057a1dd1..97c0bc3181 100644 Binary files a/pydeeptools/deeptools/test/test_data/alignmentSieve8.bam and b/pydeeptools/deeptools/test/test_data/alignmentSieve8.bam differ 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/alignmentsieve.rs b/src/alignmentsieve.rs index 70e917077a..8d16f1b75c 100644 --- a/src/alignmentsieve.rs +++ b/src/alignmentsieve.rs @@ -1,6 +1,8 @@ use crate::filehandler::{is_bed_or_gtf, read_bedfile}; use crate::filtering::Alignmentfilters; +#[cfg(feature = "python")] use pyo3::prelude::*; +#[cfg(feature = "python")] use pyo3::types::PyList; use rust_htslib::bam::record::CigarString; use rust_htslib::bam::{self, Header, Read, Reader, Writer}; @@ -9,6 +11,7 @@ use std::fs::File; use std::io::{BufWriter, Write}; use std::path::Path; +#[cfg(feature = "python")] #[pyfunction] pub fn r_alignmentsieve( py: Python, @@ -194,7 +197,7 @@ pub fn r_alignmentsieve( } let tid = record.tid() as usize; - if let Some(shifted) = apply_shift(&record, &readshift, chrom_sizes[tid]) { + if let Ok(shifted) = apply_shift(&record, &readshift, chrom_sizes[tid]) { if let Some(ref mut w) = obam { w.write(&shifted).unwrap_or_else(|e| { panic!( @@ -489,12 +492,90 @@ fn query_alignment_end(record: &bam::Record) -> i64 { read_pos } -fn apply_shift(record: &bam::Record, shift: &[i32], chrom_len: u64) -> Option { +fn cigar_op_len(op: &bam::record::Cigar) -> u32 { + use bam::record::Cigar::*; + match *op { + Match(n) | Ins(n) | Del(n) | RefSkip(n) | SoftClip(n) | HardClip(n) | Pad(n) | Equal(n) | Diff(n) => n, + } +} + +fn cigar_with_len(op: bam::record::Cigar, n: u32) -> bam::record::Cigar { + use bam::record::Cigar::*; + match op { + Match(_) => Match(n), + Ins(_) => Ins(n), + Del(_) => Del(n), + RefSkip(_) => RefSkip(n), + SoftClip(_) => SoftClip(n), + HardClip(_) => HardClip(n), + Pad(_) => Pad(n), + Equal(_) => Equal(n), + Diff(_) => Diff(n), + } +} + +fn cigar_consumes_ref(op: &bam::record::Cigar) -> bool { + use bam::record::Cigar::*; + matches!(op, Match(_) | Equal(_) | Diff(_) | Del(_) | RefSkip(_)) +} + +/// Remove `k` reference bases from the left of a CIGAR. Query bases removed from +/// M/=/X (and any I/S before them) become one leading soft clip, so the query-consuming +/// length, and therefore SEQ/QUAL, is unchanged. +fn trim_cigar_left(ops: &[bam::record::Cigar], mut k: u32) -> Vec { + use bam::record::Cigar::*; + let mut q = 0u32; + let mut rest: Vec = Vec::new(); + for (i, op) in ops.iter().enumerate() { + if k == 0 { + rest.extend_from_slice(&ops[i..]); + break; + } + let n = cigar_op_len(op); + if !cigar_consumes_ref(op) { + if matches!(op, Ins(_) | SoftClip(_)) { + q += n; + } + continue; + } + let t = n.min(k); + if matches!(op, Match(_) | Equal(_) | Diff(_)) { + q += t; + } + k -= t; + if t < n { + rest.push(cigar_with_len(*op, n - t)); + rest.extend_from_slice(&ops[i + 1..]); + break; + } + } + let mut out = Vec::with_capacity(rest.len() + 1); + if q > 0 { + out.push(SoftClip(q)); + } + out.extend(rest); + out +} + +/// Why `apply_shift` refused a record. Lets callers count each reason separately. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum ShiftRefusal { + /// Record is not part of a proper pair. + NotProperPair, + /// Fewer than 4 shift values were given. + ShortShift, + /// Shift leaves less than 1 bp of span after clamping to the chromosome. + ZeroSpan, + /// Shift would move an alignment end past its original span (start left or end right). + LeftExtension, +} + +pub fn apply_shift(record: &bam::Record, shift: &[i32], chrom_len: u64) -> Result { if !record.is_proper_pair() { - return None; + return Err(ShiftRefusal::NotProperPair); } if shift.len() < 4 { - return None; + return Err(ShiftRefusal::ShortShift); } let tlen = record.insert_size(); @@ -537,7 +618,7 @@ fn apply_shift(record: &bam::Record, shift: &[i32], chrom_len: u64) -> Option Option> = newrec - .aux_iter() - .filter_map(|r| r.ok().map(|(t, _)| t.to_vec())) - .collect(); - for t in &tags { - let _ = newrec.remove_aux(t); + // Keep SEQ, QUAL and aux tags. Only the alignment moves: trim reference bases + // from the left (shift start right) and/or the right (shift end left). + // Extending a read would need reference bases we do not have, so it is refused. + let old_start = record.pos(); + let old_end = old_start + query_alignment_end(record); + let d_start = start - old_start; + let d_end = end - old_end; + if d_start < 0 || d_end > 0 { + return Err(ShiftRefusal::LeftExtension); + } + let ops: Vec = record.cigar().iter().cloned().collect(); + let left = trim_cigar_left(&ops, d_start as u32); + let mut rev: Vec = left.into_iter().rev().collect(); + rev = trim_cigar_left(&rev, (-d_end) as u32); + let new_ops: Vec = rev.into_iter().rev().collect(); + if !new_ops.iter().any(|op| matches!(op, bam::record::Cigar::Match(_) | bam::record::Cigar::Equal(_) | bam::record::Cigar::Diff(_))) { + return Err(ShiftRefusal::ZeroSpan); } + let seq: Vec = record.seq().as_bytes(); + let qual: Vec = record.qual().to_vec(); + newrec.set(record.qname(), Some(&CigarString(new_ops)), &seq, &qual); newrec.set_pos(start); newrec.set_insert_size(new_tlen); @@ -578,7 +666,7 @@ fn apply_shift(record: &bam::Record, shift: &[i32], chrom_len: u64) -> Option String { @@ -596,3 +684,91 @@ fn smart_label(label: &str) -> String { without_ext.to_string() } } + +#[cfg(test)] +mod tests { + use super::*; + use rust_htslib::bam::record::{Aux, Cigar}; + + fn proper_pair(cigar: &[Cigar], pos: i64) -> bam::Record { + let qlen: u32 = cigar + .iter() + .map(|c| match c { + Cigar::Match(n) | Cigar::SoftClip(n) | Cigar::Ins(n) | Cigar::Equal(n) | Cigar::Diff(n) => *n, + _ => 0, + }) + .sum(); + let seq = vec![b'A'; qlen as usize]; + let qual = vec![30u8; qlen as usize]; + let mut r = bam::Record::new(); + r.set(b"read1", Some(&CigarString(cigar.to_vec())), &seq, &qual); + r.set_pos(pos); + r.set_tid(0); + r.set_mtid(0); + r.set_mpos(pos + 200); + r.set_insert_size(300); + r.set_flags(0x1 | 0x2 | 0x40); + r.push_aux(b"CB", Aux::String("AAAC-1")).unwrap(); + r + } + + #[test] + fn shift_keeps_cell_barcode_tag() { + let r = proper_pair(&[Cigar::Match(50)], 1000); + let out = apply_shift(&r, &[4, -5, 5, -4], 1_000_000).expect("shifted"); + assert!(out.aux(b"CB").is_ok(), "CB tag dropped by apply_shift"); + } + + #[test] + fn shift_keeps_soft_clip_in_cigar() { + let r = proper_pair(&[Cigar::SoftClip(5), Cigar::Match(50)], 1000); + let out = apply_shift(&r, &[4, -5, 5, -4], 1_000_000).expect("shifted"); + assert_eq!(out.pos(), 1004); + assert_eq!(out.cigar().end_pos() - out.pos(), 46, "ref span must equal end - start"); + let q: u32 = out.cigar().iter().filter(|c| !matches!(c, Cigar::Del(_) | Cigar::RefSkip(_) | Cigar::HardClip(_) | Cigar::Pad(_))).map(|c| match c { Cigar::Match(n) | Cigar::Ins(n) | Cigar::SoftClip(n) | Cigar::Equal(n) | Cigar::Diff(n) => *n, _ => 0 }).sum(); + assert_eq!(q as usize, out.seq_len(), "CIGAR query length must equal SEQ length"); + } +} + +#[cfg(test)] +mod refusal_tests { + use super::*; + use rust_htslib::bam::record::{Aux, Cigar, CigarString}; + + fn fwd_read1(pos: i64) -> bam::Record { + let seq = vec![b'A'; 50]; + let qual = vec![30u8; 50]; + let mut r = bam::Record::new(); + r.set(b"read1", Some(&CigarString(vec![Cigar::Match(50)])), &seq, &qual); + r.set_pos(pos); + r.set_tid(0); + r.set_mtid(0); + r.set_mpos(pos + 200); + r.set_insert_size(300); + r.set_flags(0x1 | 0x2 | 0x40); + r.push_aux(b"CB", Aux::String("AAAC-1")).unwrap(); + r + } + + #[test] + fn forward_read1_negative_start_shift_is_left_extension() { + let r = fwd_read1(1000); + assert_eq!( + apply_shift(&r, &[-3, 0, 0, 0], 1_000_000).unwrap_err(), + ShiftRefusal::LeftExtension + ); + } + + #[test] + fn shift_consuming_every_aligned_base_is_zero_span() { + // 1M5I1M: shift start by 2 ref bases -> only soft clip would remain + let mut r = fwd_read1(1000); + let seq = vec![b'A'; 7]; + let qual = vec![30u8; 7]; + r.set(b"read1", Some(&CigarString(vec![Cigar::Match(1), Cigar::Ins(5), Cigar::Match(1)])), &seq, &qual); + assert_eq!( + apply_shift(&r, &[2, 0, 0, 0], 1_000_000).unwrap_err(), + ShiftRefusal::ZeroSpan + ); + } +} diff --git a/src/lib.rs b/src/lib.rs index 9a37b0398e..31d1651230 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,17 +1,25 @@ +#[cfg(feature = "python")] use pyo3::prelude::*; -mod alignmentsieve; +pub mod calc; +pub mod covcalc; +pub mod filehandler; +pub mod filtering; +pub mod normalization; + +pub mod alignmentsieve; +#[cfg(feature = "python")] mod bamcompare; +#[cfg(feature = "python")] mod bamcoverage; -mod calc; +#[cfg(feature = "python")] mod computematrix; -mod covcalc; -mod filehandler; -mod filtering; +#[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)?)?;