rsomics_bed_reldist/

lib.rs

//! Relative distance statistics between two BED interval sets.
//!
//! ## Algorithm
//!
//! For each query interval in file A, finds the nearest upstream and downstream
//! intervals in file B. The relative distance is:
//!
//!   `relDist = min(distLeft, distRight) / (distLeft + distRight)`
//!
//! where `distLeft` is the distance to the B interval whose end is closest
//! upstream, and `distRight` is the distance to the B interval whose start is
//! closest downstream. Only intervals on the same chromosome are considered.
//!
//! Both files must be sorted by chromosome then start position (matching
//! `bedtools sort` / `sort -k1,1 -k2,2n` output).
//!
//! The relative distance is binned into 100 equal bins over [0.0, 0.5]. The
//! output reports, for each bin, the observed fraction and count.
//!
//! ## Reference
//!
//! Favorov et al. (2012). Exploring Massive, Genome Scale Datasets with the
//! `GenometriCorr` Package. PLOS Computational Biology 8(5): e1002529.
//! DOI: 10.1371/journal.pcbi.1002529

use std::io::{BufRead, BufReader, Read, Write};

use rsomics_common::{Result, RsomicsError};

const NUM_BINS: usize = 100;

/// One row of the reldist output table.
#[derive(Debug, Clone)]
pub struct RelDistRow {
    pub rel_dist: f64,
    pub count: u64,
    pub total: u64,
    pub fraction: f64,
}

/// Parsed BED3 record (chrom, start, end).
#[derive(Debug, Clone)]
struct Bed3 {
    chrom: String,
    start: i64,
    end: i64,
}

fn parse_bed3(line: &str) -> Option<Bed3> {
    let t = line.trim_end_matches(['\n', '\r']);
    if t.is_empty() || t.starts_with('#') || t.starts_with("track") || t.starts_with("browser") {
        return None;
    }
    let mut cols = t.splitn(4, '\t');
    let chrom = cols.next()?.to_owned();
    let start: i64 = cols.next()?.parse().ok()?;
    let end: i64 = cols.next()?.parse().ok()?;
    Some(Bed3 { chrom, start, end })
}

fn read_bed<R: Read>(r: R) -> Result<Vec<Bed3>> {
    let mut records = Vec::new();
    for line in BufReader::new(r).lines() {
        let line = line.map_err(RsomicsError::Io)?;
        if let Some(rec) = parse_bed3(&line) {
            records.push(rec);
        }
    }
    Ok(records)
}

/// Compute relative distances between query intervals in `a` and reference
/// intervals in `b`, returning a 100-bin histogram table.
///
/// Both `a` and `b` must be sorted by (chrom, start).
pub fn reldist<RA: Read, RB: Read>(a: RA, b: RB) -> Result<Vec<RelDistRow>> {
    let a_recs = read_bed(a)?;
    let b_recs = read_bed(b)?;

    let mut bins = [0u64; NUM_BINS];
    let mut total: u64 = 0;

    // Group b by chromosome for O(log n) per-record lookup.
    let mut chrom_ranges: Vec<(&str, usize, usize)> = Vec::new();
    {
        let mut i = 0;
        while i < b_recs.len() {
            let chrom = b_recs[i].chrom.as_str();
            let start = i;
            while i < b_recs.len() && b_recs[i].chrom == chrom {
                i += 1;
            }
            chrom_ranges.push((chrom, start, i));
        }
    }

    for a_rec in &a_recs {
        let Some((_ch, b_start, b_end)) = chrom_ranges
            .iter()
            .find(|(c, _, _)| *c == a_rec.chrom.as_str())
        else {
            continue;
        };
        let b_slice = &b_recs[*b_start..*b_end];
        if b_slice.is_empty() {
            continue;
        }

        let a_mid = i64::midpoint(a_rec.start, a_rec.end);

        // Find leftmost b whose start > a_mid → right neighbor.
        let right_idx = b_slice.partition_point(|b| b.start <= a_mid);

        let dist_left = if right_idx > 0 {
            let b_left = &b_slice[right_idx - 1];
            let b_mid = i64::midpoint(b_left.start, b_left.end);
            (a_mid - b_mid).max(0)
        } else {
            i64::MAX
        };

        let dist_right = if right_idx < b_slice.len() {
            let b_right = &b_slice[right_idx];
            let b_mid = i64::midpoint(b_right.start, b_right.end);
            (b_mid - a_mid).max(0)
        } else {
            i64::MAX
        };

        if dist_left == i64::MAX && dist_right == i64::MAX {
            continue;
        }

        let d_total = match (dist_left, dist_right) {
            (i64::MAX, _) => dist_right * 2,
            (_, i64::MAX) => dist_left * 2,
            (l, r) => l + r,
        };
        if d_total == 0 {
            // Overlapping — relative distance is 0 → bin 0.
            bins[0] += 1;
            total += 1;
            continue;
        }

        let min_dist = dist_left.min(dist_right);
        // Both min_dist and d_total are non-negative; rel ∈ [0.0, 0.5].
        #[allow(clippy::cast_precision_loss)]
        let rel = min_dist as f64 / d_total as f64;
        // Multiply by 2 to map [0.0, 0.5] → [0.0, 1.0], then scale to [0, NUM_BINS).
        // rel is always non-negative, so truncation is the only concern (no sign loss).
        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
        let bin = (rel * 2.0 * NUM_BINS as f64) as usize;
        let bin = bin.min(NUM_BINS - 1);
        bins[bin] += 1;
        total += 1;
    }

    let rows: Vec<RelDistRow> = (0..NUM_BINS)
        .map(|i| {
            // Bin centres: domain [0, 0.5], each bin covers 0.5/NUM_BINS.
            #[allow(clippy::cast_precision_loss)]
            let rel_dist = (i as f64 + 0.5) * 0.5 / NUM_BINS as f64;
            let count = bins[i];
            let fraction = if total > 0 {
                count as f64 / total as f64
            } else {
                0.0
            };
            RelDistRow {
                rel_dist,
                count,
                total,
                fraction,
            }
        })
        .collect();

    Ok(rows)
}

/// Write reldist results to `w` in tab-separated format (matching bedtools reldist output).
pub fn write_reldist<W: Write>(rows: &[RelDistRow], w: &mut W) -> Result<()> {
    writeln!(w, "rel_dist\tcount\ttotal\tfraction").map_err(RsomicsError::Io)?;
    for row in rows {
        writeln!(
            w,
            "{:.2}\t{}\t{}\t{:.4}",
            row.rel_dist, row.count, row.total, row.fraction
        )
        .map_err(RsomicsError::Io)?;
    }
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::io::Cursor;

    fn run(a: &str, b: &str) -> Vec<RelDistRow> {
        reldist(Cursor::new(a), Cursor::new(b)).unwrap()
    }

    #[test]
    fn output_has_100_bins() {
        let a = "chr1\t100\t200\n";
        let b = "chr1\t50\t80\nchr1\t250\t300\n";
        let rows = run(a, b);
        assert_eq!(rows.len(), 100, "must have 100 bins");
    }

    #[test]
    fn empty_a_returns_zero_counts() {
        let rows = run("", "chr1\t100\t200\n");
        let total: u64 = rows.iter().map(|r| r.count).sum();
        assert_eq!(total, 0);
    }

    #[test]
    fn empty_b_returns_zero_counts() {
        let rows = run("chr1\t100\t200\n", "");
        let total: u64 = rows.iter().map(|r| r.count).sum();
        assert_eq!(total, 0);
    }

    #[test]
    fn total_equals_num_a_records_with_neighbors() {
        let a = "chr1\t100\t200\nchr1\t400\t500\n";
        let b = "chr1\t50\t80\nchr1\t250\t300\nchr1\t600\t700\n";
        let rows = run(a, b);
        let total = rows.iter().map(|r| r.count).sum::<u64>();
        assert_eq!(total, 2);
    }

    #[test]
    fn fractions_sum_to_one_when_records_exist() {
        let a = "chr1\t100\t200\nchr1\t400\t500\nchr1\t700\t800\n";
        let b = "chr1\t50\t80\nchr1\t300\t350\nchr1\t600\t650\nchr1\t900\t950\n";
        let rows = run(a, b);
        let total_count: u64 = rows.iter().map(|r| r.count).sum();
        if total_count > 0 {
            let frac_sum: f64 = rows.iter().map(|r| r.fraction).sum();
            assert!(
                (frac_sum - 1.0).abs() < 1e-9,
                "fractions must sum to 1.0 (got {frac_sum})"
            );
        }
    }

    #[test]
    fn different_chromosomes_ignored() {
        let a = "chr1\t100\t200\n";
        let b = "chr2\t50\t80\n";
        let rows = run(a, b);
        let total: u64 = rows.iter().map(|r| r.count).sum();
        assert_eq!(total, 0, "different chromosomes must not be matched");
    }
}