fixed missing tool implementations

Author: pbenner

Cargo.toml

@@ -28,6 +28,7 @@ list_comprehension_macro = "0.1"
 mysql = "25.0.1"
 num = "0.4"
 plotters = "0.3.7"
+quick-xml = { version = "0.38", features = ["serialize"] }
 rand = "0.8.5"
 reqwest = { version = "0.12.8", features = ["blocking"] }
 serde = { version = "1.0", features = ["derive"] }

README.md

@@ -25,6 +25,8 @@ The package contains the following command line tools:
 | bam-to-fastq             | reconstruct FASTQ records from a BAM file                                |
 | bam-to-bigwig            | convert bam to bigWig (estimate fragment length if required)             |
 | bam-view                 | print contents of a bam file                                             |
+| bed-remove-overlaps      | remove BED or table rows that overlap inadmissible regions               |
+| bigwig-counts-to-quantiles | convert bigWig counts to empirical quantiles                           |
 | bigwig-edit-chrom-names  | rewrite a bigWig with chromosome names transformed by a regex            |
 | bigwig-extract           | extract bigWig data for BED regions as a table or bigWig                 |
 | bigwig-extract-chroms    | write a new bigWig containing only selected chromosomes                  |
@@ -49,6 +51,7 @@ The package contains the following command line tools:
 | pwm-scan-regions         | score genomic regions with one or more PWMs                              |
 | pwm-scan-sequences       | scan FASTA sequences with a PWM and export a bigWig track                |
 | segmentation-differential| merge and score differential chromatin states across segmentations       |
+| sequence-similarity      | compute sliding-window k-mer similarity to a reference sequence          |
 
 
 ## Examples
@@ -240,6 +243,9 @@ The result is:
 # Apply a custom shared-library function across tracks
 bigwig-map mapper.so result.bw track1.bw track2.bw
 
+# For Rust plugins, export `rustynetics_bigwig_map` and select the Rust ABI
+bigwig-map --abi rust mapper.so result.bw track1.bw track2.bw
+
 # Extract selected regions from a bigWig file as a table
 bigwig-extract signal.bw regions.bed signal.table
 
@@ -260,6 +266,62 @@ bigwig-positive result.table track1.bw:2.0 track2.bw:2.0
 bigwig-edit-chrom-names --dry-run signal.bw '^chr' ''
 ```
 
+`bigwig-map` loads a shared library and calls a mapping function for each bin.
+The mapper receives the sequence name, the genomic position, and the values from
+all input tracks for the current bin, and must return the output value for that bin.
+
+For Rust, the recommended interface is the Rust ABI plugin mode:
+
+- build the mapper as a `cdylib`
+- export the symbol `rustynetics_bigwig_map`
+- accept a pointer to `rustynetics::bigwig_map_plugin::BigWigMapInput`
+
+`bigwig-map` defaults to `--abi auto`, which first looks for the Rust symbol
+`rustynetics_bigwig_map` and then falls back to the legacy shared-library ABI
+(`F` or `bigwig_map`). You can force one mode explicitly with `--abi rust` or
+`--abi legacy`. `--symbol` overrides the default symbol name.
+
+A minimal Rust mapper looks like this:
+
+```rust
+use rustynetics::bigwig_map_plugin::BigWigMapInput;
+
+#[no_mangle]
+pub unsafe extern "C" fn rustynetics_bigwig_map(input: *const BigWigMapInput) -> f64 {
+    let input = &*input;
+    let values = input.values();
+
+    if values.is_empty() {
+        return f64::NAN;
+    }
+
+    values.iter().copied().sum::<f64>() / values.len() as f64
+}
+```
+
+The `BigWigMapInput` helper provides:
+
+- `input.seqname()` for the chromosome name
+- `input.position` for the current genomic position
+- `input.values()` for the per-track input values at that position
+
+The plugin crate needs to be compiled as a shared library:
+
+```toml
+[lib]
+crate-type = ["cdylib"]
+```
+
+Then build and run it like this:
+
+```bash
+cargo build --release
+bigwig-map --abi rust target/release/libmy_mapper.so result.bw track1.bw track2.bw
+```
+
+On macOS the shared library is usually `libmy_mapper.dylib`, and on Windows
+`my_mapper.dll`.
+
 ### Motif XML extraction
 
 ```bash

src/bbi.rs

@@ -1140,6 +1140,10 @@ impl Default for BData {
 /* -------------------------------------------------------------------------- */
 
 impl BData {
+    pub fn key_offsets(&self) -> &[i64] {
+        &self.ptr_keys
+    }
+
     pub fn add(&mut self, key: Vec<u8>, value: Vec<u8>) -> Result<(), String> {
         if key.len() as u32 != self.key_size {
             return Err("BData.Add(): key has invalid length".to_string());

src/bigwig_map_plugin.rs

@@ -0,0 +1,100 @@
+use std::ffi::c_char;
+use std::slice;
+use std::str::{self, Utf8Error};
+
+/// Default symbol name for Rust-oriented `bigwig-map` plugins.
+pub const BIGWIG_MAP_RUST_SYMBOL: &str = "rustynetics_bigwig_map";
+
+/// Compatibility symbol used by older shared-library mappers.
+pub const BIGWIG_MAP_LEGACY_SYMBOL: &str = "bigwig_map";
+
+/// Legacy `bigwig-map` ABI:
+/// `(seqname_ptr, position, values_ptr, values_len) -> f64`.
+pub type BigWigMapLegacyFn = unsafe extern "C" fn(*const c_char, usize, *const f64, usize) -> f64;
+
+/// Rust-oriented `bigwig-map` ABI:
+/// `(&BigWigMapInput) -> f64`.
+pub type BigWigMapRustFn = unsafe extern "C" fn(*const BigWigMapInput) -> f64;
+
+/// FFI payload passed to Rust-oriented `bigwig-map` plugins.
+///
+/// A plugin can depend on `rustynetics` and export:
+///
+/// ```ignore
+/// use rustynetics::bigwig_map_plugin::{BigWigMapInput, BIGWIG_MAP_RUST_SYMBOL};
+///
+/// #[no_mangle]
+/// pub unsafe extern "C" fn rustynetics_bigwig_map(input: *const BigWigMapInput) -> f64 {
+///     let input = &*input;
+///     let values = input.values();
+///     values.iter().copied().sum()
+/// }
+/// ```
+#[repr(C)]
+#[derive(Clone, Copy, Debug)]
+pub struct BigWigMapInput {
+    pub seqname_ptr: *const c_char,
+    pub seqname_len: usize,
+    pub position: usize,
+    pub values_ptr: *const f64,
+    pub values_len: usize,
+}
+
+impl BigWigMapInput {
+    /// Returns the chromosome name as raw bytes.
+    ///
+    /// # Safety
+    /// `seqname_ptr` must point to a valid buffer of length `seqname_len`.
+    pub unsafe fn seqname_bytes(&self) -> &[u8] {
+        if self.seqname_len == 0 {
+            &[]
+        } else {
+            slice::from_raw_parts(self.seqname_ptr as *const u8, self.seqname_len)
+        }
+    }
+
+    /// Returns the chromosome name as UTF-8 text.
+    ///
+    /// # Safety
+    /// `seqname_ptr` must point to a valid UTF-8 buffer of length `seqname_len`.
+    pub unsafe fn seqname(&self) -> Result<&str, Utf8Error> {
+        str::from_utf8(self.seqname_bytes())
+    }
+
+    /// Returns the input values for the current genomic position.
+    ///
+    /// # Safety
+    /// `values_ptr` must point to a valid buffer of length `values_len`.
+    pub unsafe fn values(&self) -> &[f64] {
+        if self.values_len == 0 {
+            &[]
+        } else {
+            slice::from_raw_parts(self.values_ptr, self.values_len)
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use std::ffi::CString;
+
+    #[test]
+    fn input_exposes_seqname_and_values() {
+        let seqname = CString::new("chr1").unwrap();
+        let values = [1.0, 2.5, 3.0];
+        let input = BigWigMapInput {
+            seqname_ptr: seqname.as_ptr(),
+            seqname_len: 4,
+            position: 10,
+            values_ptr: values.as_ptr(),
+            values_len: values.len(),
+        };
+
+        let parsed_seqname = unsafe { input.seqname().unwrap() };
+        let parsed_values = unsafe { input.values() };
+
+        assert_eq!(parsed_seqname, "chr1");
+        assert_eq!(parsed_values, &values);
+    }
+}

src/bin/bam-to-bigwig.rs

@@ -314,6 +314,10 @@ fn main() {
             .long("filter-chromosomes")
             .num_args(1)
             .help("Remove all reads on the given chromosomes [comma separated list]"))
+        .arg(Arg::new("remove-filtered-chromosomes")
+            .long("remove-filtered-chromosomes")
+            .action(ArgAction::SetTrue)
+            .help("Remove all chromosomes that have been filtered out"))
         // track options
         .arg(Arg::new("binning-method")
             .long("binning-method")
@@ -331,7 +335,7 @@ fn main() {
         .arg(Arg::new("pseudocounts")
             .long("pseudocounts")
             .num_args(1)
-            .help("Pseudocounts added to treatment and control signal [default: `1.0,1.0']"))
+            .help("Pseudocounts added to treatment and control signal [default: `0.0,0.0']"))
         .arg(Arg::new("smoothen-control")
             .long("smoothen-control")
             .action(ArgAction::SetTrue)
@@ -591,6 +595,9 @@ fn main() {
             opt_filter_chroms.split(',').map(String::from).collect(),
         ));
     }
+    options_list.push(OptionCoverage::RemoveFilteredChroms(
+        matches.get_flag("remove-filtered-chromosomes"),
+    ));
 
     if let Some(opt_initial_value) = matches.get_one::<String>("initial-value") {
         let value: f64 = opt_initial_value.parse().unwrap_or_else(|_| {

src/bin/bed-remove-overlaps.rs

@@ -0,0 +1,126 @@
+use std::process;
+
+use clap::{Arg, ArgAction, Command};
+
+use rustynetics::granges::GRanges;
+use rustynetics::granges_table::OptionPrintStrand;
+
+mod common;
+
+#[derive(Clone, Copy, Debug, Eq, PartialEq)]
+enum InputFormat {
+    Bed3,
+    Bed6,
+    Bed9,
+    Table,
+}
+
+fn read_table_auto(path: Option<&str>) -> Result<GRanges, Box<dyn std::error::Error>> {
+    let mut reader = common::open_reader(path)?;
+    let mut granges = GRanges::default();
+    granges.read_table_all(&mut reader)?;
+    Ok(granges)
+}
+
+fn import_input(
+    path: Option<&str>,
+    format: InputFormat,
+) -> Result<GRanges, Box<dyn std::error::Error>> {
+    match format {
+        InputFormat::Bed3 => {
+            let mut reader = common::open_reader(path)?;
+            let mut granges = GRanges::default();
+            granges.read_bed3(&mut reader)?;
+            Ok(granges)
+        }
+        InputFormat::Bed6 => {
+            let mut reader = common::open_reader(path)?;
+            let mut granges = GRanges::default();
+            granges.read_bed6(&mut reader)?;
+            Ok(granges)
+        }
+        InputFormat::Bed9 => {
+            let mut reader = common::open_reader(path)?;
+            let mut granges = GRanges::default();
+            granges.read_bed9(&mut reader)?;
+            Ok(granges)
+        }
+        InputFormat::Table => read_table_auto(path),
+    }
+}
+
+fn export_output(
+    granges: &GRanges,
+    path: Option<&str>,
+    format: InputFormat,
+) -> Result<(), Box<dyn std::error::Error>> {
+    let mut writer = common::open_writer(path)?;
+    match format {
+        InputFormat::Bed3 => granges.write_bed3(&mut writer)?,
+        InputFormat::Bed6 => granges.write_bed6(&mut writer)?,
+        InputFormat::Bed9 => granges.write_bed9(&mut writer)?,
+        InputFormat::Table => granges.write_table(&mut writer, &[&OptionPrintStrand(true)])?,
+    }
+    Ok(())
+}
+
+fn main() {
+    let matches = Command::new("bed-remove-overlaps")
+        .about("Remove regions that overlap a set of inadmissible intervals")
+        .arg(
+            Arg::new("input-format")
+                .long("input-format")
+                .default_value("bed3")
+                .value_parser(["bed3", "bed6", "bed9", "table"]),
+        )
+        .arg(
+            Arg::new("verbose")
+                .short('v')
+                .long("verbose")
+                .action(ArgAction::Count),
+        )
+        .arg(Arg::new("inadmissible").required(true).index(1))
+        .arg(Arg::new("input").index(2))
+        .arg(Arg::new("output").index(3))
+        .get_matches();
+
+    let format = match matches.get_one::<String>("input-format").unwrap().as_str() {
+        "bed3" => InputFormat::Bed3,
+        "bed6" => InputFormat::Bed6,
+        "bed9" => InputFormat::Bed9,
+        "table" => InputFormat::Table,
+        _ => unreachable!(),
+    };
+    let inadmissible = matches.get_one::<String>("inadmissible").unwrap();
+    let input = matches.get_one::<String>("input").map(String::as_str);
+    let output = matches.get_one::<String>("output").map(String::as_str);
+    let verbose = matches.get_count("verbose") > 0;
+
+    if verbose {
+        eprintln!("Reading inadmissible regions `{inadmissible}`...");
+    }
+    let remove = import_input(Some(inadmissible), InputFormat::Bed3).unwrap_or_else(|error| {
+        eprintln!("reading inadmissible regions failed: {error}");
+        process::exit(1);
+    });
+    if verbose {
+        if let Some(path) = input {
+            eprintln!("Reading input `{path}`...");
+        } else {
+            eprintln!("Reading input from stdin...");
+        }
+    }
+    let input_ranges = import_input(input, format).unwrap_or_else(|error| {
+        eprintln!("reading input failed: {error}");
+        process::exit(1);
+    });
+    let before = input_ranges.num_rows();
+    let filtered = input_ranges.remove_overlaps_with(&remove);
+    if verbose {
+        eprintln!("Removed {} regions", before - filtered.num_rows());
+    }
+    if let Err(error) = export_output(&filtered, output, format) {
+        eprintln!("writing output failed: {error}");
+        process::exit(1);
+    }
+}