rsomics_hommola/

lib.rs

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

use rayon::prelude::*;
use rsomics_common::{Result, RsomicsError};

pub mod dm;
pub mod interaction;
mod rng;

pub use dm::DistanceMatrix;
pub use interaction::Interaction;

pub struct HommolaResult {
    pub corr_coeff: f64,
    pub p_value: f64,
    pub num_hosts: usize,
    pub num_pars: usize,
    pub num_interactions: usize,
    pub num_pairs: usize,
    pub permutations: usize,
    pub perm_stats: Vec<f64>,
}

/// Hommola et al. (2009) host-parasite cospeciation test.
///
/// Builds, for every pair of interaction edges, the host patristic distance and
/// the parasite patristic distance, then takes the Pearson correlation between
/// the two resulting vectors. The coefficient is deterministic and matches
/// scikit-bio's `hommola_cospeciation` to floating-point tolerance; the p-value
/// is a seeded permutation (Monte-Carlo) estimate of the one-sided greater
/// alternative and does not reproduce numpy's Generator stream bit-for-bit.
pub fn hommola(
    host: &DistanceMatrix,
    par: &DistanceMatrix,
    inter: &Interaction,
    permutations: usize,
    seed: u64,
) -> HommolaResult {
    let num_hosts = host.n;
    let num_pars = par.n;
    let edges = inter.edges();
    let pairs = edges * (edges - 1) / 2;

    // (i,j) edge-pair list, i<j, in skbio's row-major upper-triangle order.
    let mut ei = Vec::with_capacity(pairs);
    let mut ej = Vec::with_capacity(pairs);
    for i in 0..edges {
        for j in (i + 1)..edges {
            ei.push(i);
            ej.push(j);
        }
    }

    let identity_h: Vec<usize> = (0..num_hosts).collect();
    let identity_p: Vec<usize> = (0..num_pars).collect();
    let corr_coeff = pearson_over_pairs(host, par, inter, &ei, &ej, &identity_h, &identity_p);

    let (p_value, perm_stats) = if permutations == 0 || corr_coeff.is_nan() {
        (f64::NAN, vec![f64::NAN; permutations])
    } else {
        let stats: Vec<f64> = (0..permutations)
            .into_par_iter()
            .map(|k| {
                let (mh, mp) = rng::host_par_perms(num_hosts, num_pars, seed, k as u64);
                pearson_over_pairs(host, par, inter, &ei, &ej, &mh, &mp)
            })
            .collect();
        let extreme = stats.iter().filter(|&&s| s >= corr_coeff).count();
        let p = (extreme + 1) as f64 / (permutations + 1) as f64;
        (p, stats)
    };

    HommolaResult {
        corr_coeff,
        p_value,
        num_hosts,
        num_pars,
        num_interactions: edges,
        num_pairs: pairs,
        permutations,
        perm_stats,
    }
}

/// Pearson r of (host-dist, par-dist) over every edge pair, with hosts and
/// parasites relabelled through `mh`/`mp`. Streaming sums, no per-call vector.
fn pearson_over_pairs(
    host: &DistanceMatrix,
    par: &DistanceMatrix,
    inter: &Interaction,
    ei: &[usize],
    ej: &[usize],
    mh: &[usize],
    mp: &[usize],
) -> f64 {
    let nh = host.n;
    let np = par.n;
    let hi = &inter.host_idx;
    let pi = &inter.par_idx;

    let mut sx = 0.0;
    let mut sy = 0.0;
    let mut sxx = 0.0;
    let mut syy = 0.0;
    let mut sxy = 0.0;
    let count = ei.len() as f64;

    for (&a, &b) in ei.iter().zip(ej) {
        let hx = host.data[mh[hi[a]] * nh + mh[hi[b]]];
        let py = par.data[mp[pi[a]] * np + mp[pi[b]]];
        sx += hx;
        sy += py;
        sxx += hx * hx;
        syy += py * py;
        sxy += hx * py;
    }

    let cov = sxy - sx * sy / count;
    let vx = sxx - sx * sx / count;
    let vy = syy - sy * sy / count;
    let denom = (vx * vy).sqrt();
    if denom == 0.0 {
        f64::NAN
    } else {
        (cov / denom).clamp(-1.0, 1.0)
    }
}

pub fn read_matrix<R: BufRead>(reader: R, source: &str) -> Result<DistanceMatrix> {
    DistanceMatrix::read(reader, source)
}

pub fn read_interaction<R: BufRead>(reader: R, source: &str) -> Result<Interaction> {
    Interaction::read(reader, source)
}

pub fn write_result<W: Write>(out: &mut W, res: &HommolaResult) -> Result<()> {
    writeln!(
        out,
        "statistic\tp_value\tnum_hosts\tnum_parasites\tnum_interactions\tnum_pairs\tpermutations"
    )
    .map_err(RsomicsError::Io)?;
    writeln!(
        out,
        "{:.12}\t{}\t{}\t{}\t{}\t{}\t{}",
        res.corr_coeff,
        fmt_p(res.p_value),
        res.num_hosts,
        res.num_pars,
        res.num_interactions,
        res.num_pairs,
        res.permutations,
    )
    .map_err(RsomicsError::Io)?;
    Ok(())
}

fn fmt_p(p: f64) -> String {
    if p.is_nan() {
        "nan".to_string()
    } else {
        format!("{p:.12}")
    }
}

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

    fn dm(rows: &[&[f64]]) -> DistanceMatrix {
        let n = rows.len();
        let mut data = vec![0.0; n * n];
        for (i, r) in rows.iter().enumerate() {
            for (j, &v) in r.iter().enumerate() {
                data[i * n + j] = v;
            }
        }
        DistanceMatrix {
            ids: (0..n).map(|i| i.to_string()).collect(),
            data,
            n,
        }
    }

    /// interaction edges as (parasite_row, host_col) pairs.
    fn inter(edges: &[(usize, usize)], num_hosts: usize, num_pars: usize) -> Interaction {
        Interaction {
            host_ids: (0..num_hosts).map(|i| i.to_string()).collect(),
            par_ids: (0..num_pars).map(|i| i.to_string()).collect(),
            par_idx: edges.iter().map(|&(p, _)| p).collect(),
            host_idx: edges.iter().map(|&(_, h)| h).collect(),
        }
    }

    /// skbio docstring example: corr_coeff == 0.832 (3dp).
    #[test]
    fn skbio_doc_example() {
        let hdist = dm(&[
            &[0.0, 3.0, 8.0, 8.0, 9.0],
            &[3.0, 0.0, 7.0, 7.0, 8.0],
            &[8.0, 7.0, 0.0, 6.0, 7.0],
            &[8.0, 7.0, 6.0, 0.0, 3.0],
            &[9.0, 8.0, 7.0, 3.0, 0.0],
        ]);
        let pdist = dm(&[
            &[0.0, 5.0, 8.0, 8.0, 8.0],
            &[5.0, 0.0, 7.0, 7.0, 7.0],
            &[8.0, 7.0, 0.0, 4.0, 4.0],
            &[8.0, 7.0, 4.0, 0.0, 2.0],
            &[8.0, 7.0, 4.0, 2.0, 0.0],
        ]);
        // interaction rows=parasites, cols=hosts, nonzero row-major:
        // [[1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0],[0,0,0,1,0],[0,0,0,1,1]]
        let edges = [(0, 0), (1, 1), (2, 2), (3, 3), (4, 3), (4, 4)];
        let it = inter(&edges, 5, 5);
        let res = hommola(&hdist, &pdist, &it, 0, 42);
        assert!(
            (res.corr_coeff - 0.832).abs() < 5e-4,
            "corr={}",
            res.corr_coeff
        );
        assert_eq!(res.num_interactions, 6);
        assert_eq!(res.num_pairs, 15);
        assert!(res.p_value.is_nan());
    }

    #[test]
    fn strong_signal_is_significant() {
        let hdist = dm(&[
            &[0.0, 3.0, 8.0, 8.0, 9.0],
            &[3.0, 0.0, 7.0, 7.0, 8.0],
            &[8.0, 7.0, 0.0, 6.0, 7.0],
            &[8.0, 7.0, 6.0, 0.0, 3.0],
            &[9.0, 8.0, 7.0, 3.0, 0.0],
        ]);
        let pdist = dm(&[
            &[0.0, 5.0, 8.0, 8.0, 8.0],
            &[5.0, 0.0, 7.0, 7.0, 7.0],
            &[8.0, 7.0, 0.0, 4.0, 4.0],
            &[8.0, 7.0, 4.0, 0.0, 2.0],
            &[8.0, 7.0, 4.0, 2.0, 0.0],
        ]);
        let edges = [(0, 0), (1, 1), (2, 2), (3, 3), (4, 3), (4, 4)];
        let it = inter(&edges, 5, 5);
        let res = hommola(&hdist, &pdist, &it, 999, 42);
        assert!(res.p_value <= 0.05, "p={}", res.p_value);
        assert_eq!(res.perm_stats.len(), 999);
    }
}