limma/

diffsplice.rs

//! Differential exon usage. Port of limma's `diffSplice` and `topSplice`
//! (`diffSplice.R`, `topSplice.R`).
//!
//! `diff_splice` takes an exon-level [`MArrayLM`] fit (one row per exon) plus a
//! gene-id label per exon, and tests whether each exon's log-fold-change departs
//! from its gene's average — i.e. whether exon usage changes between conditions.
//! It returns exon-level moderated t-statistics and gene-level F / Simes /
//! Bonferroni summaries. `top_splice` collates a single coefficient's results
//! into a ranked table, mirroring `topSplice`.
//!
//! The annotation-frame plumbing of R's `diffSplice` (carrying through arbitrary
//! `genes` columns and detecting which are gene-level) is out of scope for this
//! numeric port: gene and exon identity are carried as plain string labels.

use anyhow::{bail, Result};
use ndarray::{Array1, Array2};

use crate::ebayes::squeeze_var_legacy;
use crate::fit::MArrayLM;
use crate::special::{f_sf, t_sf};
use crate::toptable::p_adjust_bh;

/// Result of [`diff_splice`]. Exon-level fields have one row per *kept* exon
/// (exons belonging to genes with at least two exons), in gene-sorted order;
/// gene-level fields have one row per kept gene.
#[derive(Clone, Debug)]
pub struct DiffSplice {
    pub coef_names: Vec<String>,

    // --- Exon level (kept exons, gene-sorted) ---
    /// Exon identifier per kept exon (defaults to the 1-based original row index).
    pub exon_id: Vec<String>,
    /// Gene identifier per kept exon.
    pub exon_geneid: Vec<String>,
    /// 0-based index into the kept-gene arrays for each kept exon.
    pub exon_gene_index: Vec<usize>,
    /// Leverage-adjusted exon log-fold-changes (relative to the gene average).
    pub coefficients: Array2<f64>,
    /// Exon moderated t-statistics.
    pub t: Array2<f64>,
    /// Two-sided exon p-values.
    pub p_value: Array2<f64>,

    // --- Gene level (kept genes) ---
    pub gene_id: Vec<String>,
    pub gene_nexons: Vec<usize>,
    pub gene_df_prior: Array1<f64>,
    pub gene_df_residual: Array1<f64>,
    pub gene_df_total: Array1<f64>,
    pub gene_s2: Array1<f64>,
    pub gene_s2_post: Array1<f64>,
    /// Gene-level F-statistics (one column per coefficient).
    pub gene_f: Array2<f64>,
    pub gene_f_p_value: Array2<f64>,
    pub gene_simes_p_value: Array2<f64>,
    pub gene_bonferroni_p_value: Array2<f64>,
    /// 0-based first/last (inclusive) kept-exon row index for each kept gene.
    pub gene_firstexon: Vec<usize>,
    pub gene_lastexon: Vec<usize>,
}

/// `diffSplice` — test for differential exon usage from an exon-level fit.
///
/// * `fit` — an [`MArrayLM`] with one row per exon (from `lmFit` on exon data).
/// * `geneid` — gene label for each exon (length = `fit` rows).
/// * `exonid` — optional exon label per exon; defaults to 1-based row indices.
/// * `robust` / `legacy` — forwarded to `squeezeVar` (R defaults: `FALSE`).
pub fn diff_splice(
    fit: &MArrayLM,
    geneid: &[String],
    exonid: Option<&[String]>,
    robust: bool,
    legacy: bool,
) -> Result<DiffSplice> {
    let n_exons = fit.coefficients.nrows();
    let n_coef = fit.coefficients.ncols();
    if geneid.len() != n_exons {
        bail!(
            "geneid has length {} but fit has {} exons",
            geneid.len(),
            n_exons
        );
    }
    if let Some(eid) = exonid {
        if eid.len() != n_exons {
            bail!(
                "exonid has length {} but fit has {} exons",
                eid.len(),
                n_exons
            );
        }
    }

    // Exon labels default to the 1-based original row index (R: ExonID=1:nrow).
    let exon_label: Vec<String> = match exonid {
        Some(eid) => eid.to_vec(),
        None => (1..=n_exons).map(|i| i.to_string()).collect(),
    };

    // Sort exons by gene (then exon label), stably — matching R's stable radix
    // `order(geneid)` / `order(geneid, exonid)`.
    let mut order: Vec<usize> = (0..n_exons).collect();
    match exonid {
        Some(eid) => {
            order.sort_by(|&a, &b| geneid[a].cmp(&geneid[b]).then_with(|| eid[a].cmp(&eid[b])))
        }
        None => order.sort_by(|&a, &b| geneid[a].cmp(&geneid[b])),
    }

    let gid: Vec<&str> = order.iter().map(|&i| geneid[i].as_str()).collect();
    let elabel: Vec<String> = order.iter().map(|&i| exon_label[i].clone()).collect();

    let mut exon_s2: Vec<f64> = order.iter().map(|&i| fit.sigma[i].powi(2)).collect();
    let exon_dfres: Vec<f64> = order.iter().map(|&i| fit.df_residual[i]).collect();
    // Zero residual variance where residual df vanishes (R: exon.s2[df<1e-6] <- 0).
    if exon_dfres.iter().cloned().fold(f64::INFINITY, f64::min) < 1e-6 {
        for i in 0..n_exons {
            if exon_dfres[i] < 1e-6 {
                exon_s2[i] = 0.0;
            }
        }
    }

    // Contiguous gene groups (sorted order => first-appearance == sorted order,
    // matching rowsum(..., reorder=FALSE)).
    let mut group_start: Vec<usize> = Vec::new();
    let mut group_len: Vec<usize> = Vec::new();
    let mut group_id: Vec<&str> = Vec::new();
    let mut i = 0;
    while i < n_exons {
        let g = gid[i];
        let start = i;
        while i < n_exons && gid[i] == g {
            i += 1;
        }
        group_id.push(g);
        group_start.push(start);
        group_len.push(i - start);
    }
    let n_genes_all = group_id.len();

    // Gene-wise exon counts, residual df and pooled residual variances.
    let gene_dfres_all: Vec<f64> = (0..n_genes_all)
        .map(|gi| {
            let (s, l) = (group_start[gi], group_len[gi]);
            (s..s + l).map(|k| exon_dfres[k]).sum()
        })
        .collect();
    let gene_s2_all: Vec<f64> = (0..n_genes_all)
        .map(|gi| {
            let (s, l) = (group_start[gi], group_len[gi]);
            let num: f64 = (s..s + l).map(|k| exon_dfres[k] * exon_s2[k]).sum();
            let den: f64 = (s..s + l).map(|k| exon_dfres[k]).sum();
            num / den
        })
        .collect();

    // Empirical-Bayes posterior gene variances (legacy=FALSE forces unequal-df).
    let squeeze = squeeze_var_legacy(
        &Array1::from(gene_s2_all.clone()),
        &Array1::from(gene_dfres_all.clone()),
        None,
        robust,
        Some(legacy),
    )?;

    // Keep genes with more than one exon.
    let kept_idx: Vec<usize> = (0..n_genes_all).filter(|&gi| group_len[gi] > 1).collect();
    let ngenes = kept_idx.len();
    if ngenes == 0 {
        bail!("no genes with more than one exon");
    }

    let gene_nexons: Vec<usize> = kept_idx.iter().map(|&gi| group_len[gi]).collect();
    let gene_df_test: Vec<f64> = gene_nexons.iter().map(|&n| (n - 1) as f64).collect();
    let gene_df_residual: Vec<f64> = kept_idx.iter().map(|&gi| gene_dfres_all[gi]).collect();
    let gene_s2_kept: Vec<f64> = kept_idx.iter().map(|&gi| gene_s2_all[gi]).collect();
    let gene_s2_post: Vec<f64> = kept_idx.iter().map(|&gi| squeeze.var_post[gi]).collect();
    let df_prior_kept: Vec<f64> = if squeeze.df_prior.len() > 1 {
        kept_idx.iter().map(|&gi| squeeze.df_prior[gi]).collect()
    } else {
        vec![squeeze.df_prior[0]; ngenes]
    };
    // df.total = df.residual + df.prior, capped at the total kept residual df.
    let sum_dfres_kept: f64 = gene_df_residual.iter().sum();
    let gene_df_total: Vec<f64> = (0..ngenes)
        .map(|gi| (gene_df_residual[gi] + df_prior_kept[gi]).min(sum_dfres_kept))
        .collect();

    // Kept exons + their kept-gene index (contiguous per gene).
    let mut kept_rows: Vec<usize> = Vec::new();
    let mut g: Vec<usize> = Vec::new();
    let mut gene_firstexon: Vec<usize> = Vec::with_capacity(ngenes);
    let mut gene_lastexon: Vec<usize> = Vec::with_capacity(ngenes);
    for (new_gi, &gi) in kept_idx.iter().enumerate() {
        let (s, l) = (group_start[gi], group_len[gi]);
        gene_firstexon.push(kept_rows.len());
        for k in s..s + l {
            kept_rows.push(k);
            g.push(new_gi);
        }
        gene_lastexon.push(kept_rows.len() - 1);
    }
    let n_kept = kept_rows.len();

    let mut ecoef = Array2::<f64>::zeros((n_kept, n_coef));
    let mut esdu = Array2::<f64>::zeros((n_kept, n_coef));
    for (r, &k) in kept_rows.iter().enumerate() {
        // Gather the kept exon's row straight from the source matrices — no
        // per-row temporary Vec.
        let orig = order[k];
        for c in 0..n_coef {
            ecoef[[r, c]] = fit.coefficients[[orig, c]];
            esdu[[r, c]] = fit.stdev_unscaled[[orig, c]];
        }
    }
    let exon_geneid: Vec<String> = kept_rows.iter().map(|&k| gid[k].to_string()).collect();
    let exon_id_kept: Vec<String> = kept_rows.iter().map(|&k| elabel[k].clone()).collect();

    // Inverse-variance weights and per-gene weighted mean coefficient.
    let mut u2 = Array2::<f64>::zeros((n_kept, n_coef));
    let mut u2_rowsum = Array2::<f64>::zeros((ngenes, n_coef));
    let mut cu2_rowsum = Array2::<f64>::zeros((ngenes, n_coef));
    for r in 0..n_kept {
        let gg = g[r];
        for c in 0..n_coef {
            let w = 1.0 / (esdu[[r, c]] * esdu[[r, c]]);
            u2[[r, c]] = w;
            u2_rowsum[[gg, c]] += w;
            cu2_rowsum[[gg, c]] += ecoef[[r, c]] * w;
        }
    }

    // Center exon coefficients on the gene mean; moderated t before leverage.
    // Every cell of `exon_coef` is assigned in the loop below, so start from
    // zeros rather than cloning `ecoef`.
    let mut exon_coef = Array2::<f64>::zeros((n_kept, n_coef));
    let mut exon_t = Array2::<f64>::zeros((n_kept, n_coef));
    for r in 0..n_kept {
        let gg = g[r];
        let sp = gene_s2_post[gg].sqrt();
        for c in 0..n_coef {
            let centered = ecoef[[r, c]] - cu2_rowsum[[gg, c]] / u2_rowsum[[gg, c]];
            exon_coef[[r, c]] = centered;
            exon_t[[r, c]] = centered / esdu[[r, c]] / sp;
        }
    }
    // Gene F = mean of squared exon t over the gene's test df.
    let mut gene_f = Array2::<f64>::zeros((ngenes, n_coef));
    for r in 0..n_kept {
        let gg = g[r];
        for c in 0..n_coef {
            gene_f[[gg, c]] += exon_t[[r, c]] * exon_t[[r, c]];
        }
    }
    for gg in 0..ngenes {
        for c in 0..n_coef {
            gene_f[[gg, c]] /= gene_df_test[gg];
        }
    }
    // Leverage rescaling of exon coefficients / t, then p-values.
    let mut exon_p = Array2::<f64>::zeros((n_kept, n_coef));
    for r in 0..n_kept {
        let gg = g[r];
        let df = gene_df_total[gg];
        for c in 0..n_coef {
            let one_m_lev = 1.0 - u2[[r, c]] / u2_rowsum[[gg, c]];
            exon_coef[[r, c]] /= one_m_lev;
            exon_t[[r, c]] /= one_m_lev.sqrt();
            exon_p[[r, c]] = 2.0 * t_sf(exon_t[[r, c]].abs(), df);
        }
    }
    let mut gene_f_p = Array2::<f64>::zeros((ngenes, n_coef));
    for gg in 0..ngenes {
        for c in 0..n_coef {
            gene_f_p[[gg, c]] = f_sf(gene_f[[gg, c]], gene_df_test[gg], gene_df_total[gg]);
        }
    }

    // Simes and Bonferroni combination of exon p-values per gene per coef.
    let mut gene_simes_p = Array2::<f64>::zeros((ngenes, n_coef));
    let mut gene_bonf_p = Array2::<f64>::zeros((ngenes, n_coef));
    for gg in 0..ngenes {
        let m = gene_nexons[gg] as f64;
        let (lo, hi) = (gene_firstexon[gg], gene_lastexon[gg]);
        for c in 0..n_coef {
            let mut ps: Vec<f64> = (lo..=hi).map(|r| exon_p[[r, c]]).collect();
            ps.sort_by(|a, b| a.partial_cmp(b).unwrap());
            let mut simes = f64::INFINITY;
            for (k, &p) in ps.iter().enumerate() {
                let adj = p * m / ((k + 1) as f64);
                if adj < simes {
                    simes = adj;
                }
            }
            gene_simes_p[[gg, c]] = simes;
            gene_bonf_p[[gg, c]] = (ps[0] * m).min(1.0);
        }
    }

    Ok(DiffSplice {
        coef_names: fit.coef_names.clone(),
        exon_id: exon_id_kept,
        exon_geneid,
        exon_gene_index: g,
        coefficients: exon_coef,
        t: exon_t,
        p_value: exon_p,
        gene_id: kept_idx
            .iter()
            .map(|&gi| group_id[gi].to_string())
            .collect(),
        gene_nexons,
        gene_df_prior: Array1::from(df_prior_kept),
        gene_df_residual: Array1::from(gene_df_residual),
        gene_df_total: Array1::from(gene_df_total),
        gene_s2: Array1::from(gene_s2_kept),
        gene_s2_post: Array1::from(gene_s2_post),
        gene_f,
        gene_f_p_value: gene_f_p,
        gene_simes_p_value: gene_simes_p,
        gene_bonferroni_p_value: gene_bonf_p,
        gene_firstexon,
        gene_lastexon,
    })
}

/// Which statistic `top_splice` ranks on. Mirrors `topSplice`'s `test` argument.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum SpliceTest {
    /// Gene-level Simes combination of exon p-values.
    Simes,
    /// Gene-level F-test across exons.
    F,
    /// Exon-level moderated t-test.
    T,
    /// Exon-level TREAT test relative to a log-fold-change threshold.
    Treat,
}

/// `top_splice` row ordering. Mirrors `topSplice`'s `sort.by` argument.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum SpliceSort {
    /// Ascending p-value.
    P,
    /// Original (diffSplice) order.
    None,
    /// Descending absolute logFC (exon t-test only).
    LogFC,
    /// Descending exon count, ties broken by ascending p-value (gene tests only).
    NExons,
}

/// One row of a [`top_splice`] table.
#[derive(Clone, Debug)]
pub struct TopSpliceRow {
    /// Exon id (t / treat tests) or gene id (F / simes tests).
    pub id: String,
    /// Gene id (equals `id` for gene-level tests).
    pub geneid: String,
    pub nexons: Option<usize>,
    pub logfc: Option<f64>,
    pub t: Option<f64>,
    pub f: Option<f64>,
    pub p_value: f64,
    pub fdr: f64,
}

/// `topSplice` — collate one coefficient's `diffSplice` results into a ranked
/// table. `coef` is a 0-based column index. `treat_lfc > 0` with
/// [`SpliceTest::T`] switches to the TREAT variant (as R does).
pub fn top_splice(
    ds: &DiffSplice,
    coef: usize,
    test: SpliceTest,
    number: usize,
    fdr: f64,
    sort_by: SpliceSort,
    treat_lfc: f64,
) -> Result<Vec<TopSpliceRow>> {
    let n_coef = ds.gene_f.ncols();
    if coef >= n_coef {
        bail!("coef index {} out of range (have {})", coef, n_coef);
    }
    if matches!(sort_by, SpliceSort::LogFC) && !matches!(test, SpliceTest::T) {
        bail!("sorting by logFC only available with the exon t-test");
    }
    if matches!(sort_by, SpliceSort::NExons) && matches!(test, SpliceTest::T) {
        bail!("sorting by NExons only available with gene-level tests");
    }
    // R promotes test="t" with a positive treat.lfc to the TREAT variant.
    let test = if matches!(test, SpliceTest::T) && treat_lfc > 0.0 {
        SpliceTest::Treat
    } else {
        test
    };

    let mut rows: Vec<TopSpliceRow> = match test {
        SpliceTest::T | SpliceTest::Treat => (0..ds.coefficients.nrows())
            .map(|r| {
                let logfc = ds.coefficients[[r, coef]];
                let t = ds.t[[r, coef]];
                let p = if matches!(test, SpliceTest::Treat) {
                    diff_splice_treat(ds, r, coef, treat_lfc)
                } else {
                    ds.p_value[[r, coef]]
                };
                TopSpliceRow {
                    id: ds.exon_id[r].clone(),
                    geneid: ds.exon_geneid[r].clone(),
                    nexons: None,
                    logfc: Some(logfc),
                    t: Some(t),
                    f: None,
                    p_value: p,
                    fdr: f64::NAN,
                }
            })
            .collect(),
        SpliceTest::F => (0..ds.gene_id.len())
            .map(|gg| TopSpliceRow {
                id: ds.gene_id[gg].clone(),
                geneid: ds.gene_id[gg].clone(),
                nexons: Some(ds.gene_nexons[gg]),
                logfc: None,
                t: None,
                f: Some(ds.gene_f[[gg, coef]]),
                p_value: ds.gene_f_p_value[[gg, coef]],
                fdr: f64::NAN,
            })
            .collect(),
        SpliceTest::Simes => (0..ds.gene_id.len())
            .map(|gg| TopSpliceRow {
                id: ds.gene_id[gg].clone(),
                geneid: ds.gene_id[gg].clone(),
                nexons: Some(ds.gene_nexons[gg]),
                logfc: None,
                t: None,
                f: None,
                p_value: ds.gene_simes_p_value[[gg, coef]],
                fdr: f64::NAN,
            })
            .collect(),
    };

    // BH FDR over the full set, then optional FDR filtering.
    let pv: Vec<f64> = rows.iter().map(|r| r.p_value).collect();
    let adj = p_adjust_bh(&pv);
    for (r, a) in rows.iter_mut().zip(adj) {
        r.fdr = a;
    }
    if fdr < 1.0 {
        rows.retain(|r| r.fdr <= fdr);
    }

    // R: `number <- min(number, nrow)`; `if(number <= 1) return(out)` unsorted.
    let number = number.min(rows.len());
    if number <= 1 {
        return Ok(rows);
    }

    // Stable sort by the requested key (matching R's stable `order`).
    let mut idx: Vec<usize> = (0..rows.len()).collect();
    match sort_by {
        SpliceSort::P => {
            idx.sort_by(|&a, &b| rows[a].p_value.partial_cmp(&rows[b].p_value).unwrap())
        }
        SpliceSort::LogFC => idx.sort_by(|&a, &b| {
            rows[b]
                .logfc
                .unwrap()
                .abs()
                .partial_cmp(&rows[a].logfc.unwrap().abs())
                .unwrap()
        }),
        SpliceSort::NExons => idx.sort_by(|&a, &b| {
            rows[b]
                .nexons
                .unwrap()
                .cmp(&rows[a].nexons.unwrap())
                .then_with(|| rows[a].p_value.partial_cmp(&rows[b].p_value).unwrap())
        }),
        SpliceSort::None => {}
    }
    idx.truncate(number);
    Ok(idx.into_iter().map(|i| rows[i].clone()).collect())
}

/// `.diffSpliceTreat`: exon-level TREAT p-value relative to `lfc`.
fn diff_splice_treat(ds: &DiffSplice, r: usize, coef: usize, lfc: f64) -> f64 {
    let lfc = lfc.abs();
    let acoef = ds.coefficients[[r, coef]].abs();
    let se = acoef / ds.t[[r, coef]].abs();
    let df = ds.gene_df_total[ds.exon_gene_index[r]];
    let p = t_sf((acoef - lfc) / se, df) + t_sf((acoef + lfc) / se, df);
    if p.is_nan() {
        1.0
    } else {
        p
    }
}

#[cfg(test)]
// Reference literals are copied verbatim from R's 17-significant-digit output
// (scratch/diffsplice_ref.R); the redundant final digit is intentional.
#[allow(clippy::excessive_precision, clippy::approx_constant)]
mod tests {
    use super::*;
    use crate::fit::lmfit;

    // Reference values from scratch/diffsplice_ref.R (limma 3.68.3, R 4.6.0).
    fn fixture() -> DiffSplice {
        #[rustfmt::skip]
        let e = Array2::from_shape_vec((12, 6), vec![
            5.1, 4.8, 5.3, 7.2, 7.0, 7.5,
            3.2, 3.5, 3.0, 3.1, 3.4, 3.3,
            6.0, 6.2, 5.8, 4.1, 4.0, 4.3,
            8.0, 8.1, 7.9, 8.2, 8.0, 8.1,
            2.0, 2.2, 1.9, 4.0, 4.1, 3.8,
            4.4, 4.6, 4.2, 5.0, 5.2, 4.9,
            7.1, 7.0, 7.3, 7.0, 6.9, 7.2,
            3.3, 3.1, 3.5, 6.0, 6.2, 5.8,
            5.5, 5.7, 5.3, 5.4, 5.6, 5.2,
            6.6, 6.4, 6.8, 6.5, 6.7, 6.3,
            4.0, 4.2, 3.8, 7.0, 7.1, 6.9,
            9.0, 9.1, 8.9, 6.0, 5.9, 6.1,
        ]).unwrap();
        #[rustfmt::skip]
        let design = Array2::from_shape_vec((6, 2), vec![
            1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
        ]).unwrap();
        let geneid: Vec<String> = [
            "G1", "G1", "G1", "G2", "G2", "G3", "G3", "G3", "G3", "G4", "G5", "G5",
        ]
        .iter()
        .map(|s| s.to_string())
        .collect();
        let names: Vec<String> = (1..=12).map(|i| format!("ex{i}")).collect();
        let coefs = vec!["Intercept".to_string(), "grpB".to_string()];
        let fit = lmfit(&e, &design, names, coefs).unwrap();
        diff_splice(&fit, &geneid, None, false, false).unwrap()
    }

    fn rclose(got: f64, want: f64) -> bool {
        if want.is_nan() {
            return got.is_nan();
        }
        (got - want).abs() <= 1e-6 * want.abs().max(1e-300) + 1e-12
    }

    fn assert_mat(got: &Array2<f64>, want: &[[f64; 2]]) {
        assert_eq!(got.nrows(), want.len());
        for (r, row) in want.iter().enumerate() {
            for c in 0..2 {
                assert!(
                    rclose(got[[r, c]], row[c]),
                    "mismatch at [{r},{c}]: got {} want {}",
                    got[[r, c]],
                    row[c]
                );
            }
        }
    }

    #[test]
    fn diff_splice_gene_level_matches_r() {
        let ds = fixture();
        assert_eq!(ds.gene_id, ["G1", "G2", "G3", "G5"]);
        assert_eq!(ds.gene_nexons, [3, 2, 4, 2]);
        assert_eq!(ds.gene_df_residual.to_vec(), [12.0, 8.0, 16.0, 8.0]);
        assert_eq!(ds.gene_df_total.to_vec(), [44.0, 44.0, 44.0, 44.0]);
        assert!(rclose(ds.gene_df_prior[0], 8134.8447803826621));

        let s2 = [
            0.046111111111111179,
            0.016666666666666594,
            0.033750000000000002,
            0.01749999999999995,
        ];
        let s2p = [
            0.032952104211060054,
            0.032916712208519362,
            0.032934297253120325,
            0.032917530923212625,
        ];
        for i in 0..4 {
            assert!(rclose(ds.gene_s2[i], s2[i]));
            assert!(rclose(ds.gene_s2_post[i], s2p[i]));
        }

        assert_mat(
            &ds.gene_f,
            &[
                [180.36278640252982, 185.3363349152192],
                [1622.3268693537839, 76.582172525751389],
                [243.61635141837323, 79.805356499122354],
                [1139.2105953352645, 820.2316286413901],
            ],
        );
        assert_mat(
            &ds.gene_f_p_value,
            &[
                [6.2867278756513315e-22, 3.6850195711806162e-22],
                [2.2936973401206709e-36, 3.4576638739475986e-11],
                [2.0512282816641372e-27, 7.91897383774075e-18],
                [4.3073010662264811e-33, 4.3519465964385832e-30],
            ],
        );
        assert_mat(
            &ds.gene_simes_p_value,
            &[
                [2.3179315147636661e-21, 1.8221109981110334e-20],
                [2.2936973401206709e-36, 3.4576638739475986e-11],
                [3.08052592995522e-25, 3.0664863233974733e-18],
                [4.3073010662264811e-33, 4.3519465964385832e-30],
            ],
        );
        assert_mat(
            &ds.gene_bonferroni_p_value,
            &[
                [2.3179315147636661e-21, 1.8221109981110334e-20],
                [4.5873946802413418e-36, 6.9153277478951973e-11],
                [3.08052592995522e-25, 3.0664863233974733e-18],
                [8.6146021324529622e-33, 8.7038931928770445e-30],
            ],
        );
    }

    #[test]
    fn diff_splice_exon_level_matches_r() {
        let ds = fixture();
        assert_eq!(
            ds.exon_geneid,
            ["G1", "G1", "G1", "G2", "G2", "G3", "G3", "G3", "G3", "G5", "G5"]
        );
        assert_eq!(
            ds.exon_id,
            ["1", "2", "3", "4", "5", "6", "7", "8", "9", "11", "12"]
        );
        assert_mat(
            &ds.coefficients,
            &[
                [0.44999999999999707, 3.0833333333333366],
                [-2.2999999999999985, -0.1166666666666679],
                [1.8500000000000014, -2.966666666666669],
                [5.9666666666666686, -1.8333333333333348],
                [-5.9666666666666694, 1.8333333333333348],
                [-0.9111111111111091, -0.20000000000000048],
                [2.7333333333333307, -1.1777777777777767],
                [-2.3777777777777787, 2.5555555555555558],
                [0.55555555555555591, -1.1777777777777783],
                [-5.0000000000000036, 6.0000000000000018],
                [5.0000000000000036, -6.0000000000000027],
            ],
        );
        assert_mat(
            &ds.t,
            &[
                [3.5057903027150417, 16.985522142464994],
                [-17.918483769432541, -0.64269543241760041],
                [14.412693466717498, -16.342826710047394],
                [40.278118989766433, -8.7511240721264709],
                [-40.27811898976644, 8.7511240721264709],
                [-7.5307529659810823, -1.1689126440774047],
                [22.592258897943278, -6.883596681789137],
                [-19.653428472194584, 14.93610600765569],
                [4.5919225402323809, -6.8835966817891467],
                [-33.75219393365807, 28.639686252495675],
                [33.75219393365807, -28.639686252495682],
            ],
        );
        assert_mat(
            &ds.p_value,
            &[
                [0.0010604632227702825, 6.0737033270367781e-21],
                [7.7264383825455539e-22, 0.52375685669844319],
                [2.8381673681066256e-18, 2.6414727327492922e-20],
                [2.2936973401206709e-36, 3.4576638739475986e-11],
                [2.2936973401206709e-36, 3.4576638739475986e-11],
                [1.9210853211775688e-09, 0.2487328076651891],
                [7.7013148248880499e-26, 1.6935989640586146e-08],
                [2.0608265277598819e-23, 7.6662158084936832e-19],
                [3.6630572076199251e-05, 1.6935989640585597e-08],
                [4.3073010662264811e-33, 4.3519465964385832e-30],
                [4.3073010662264811e-33, 4.3519465964385222e-30],
            ],
        );
    }

    #[test]
    fn top_splice_f_matches_r() {
        let ds = fixture();
        let top = top_splice(&ds, 1, SpliceTest::F, 100, 1.0, SpliceSort::P, 0.0).unwrap();
        let ids: Vec<&str> = top.iter().map(|r| r.id.as_str()).collect();
        assert_eq!(ids, ["G5", "G1", "G3", "G2"]);
        let f = [
            820.2316286413901,
            185.3363349152192,
            79.805356499122354,
            76.582172525751389,
        ];
        let p = [
            4.3519465964385832e-30,
            3.6850195711806162e-22,
            7.91897383774075e-18,
            3.4576638739475986e-11,
        ];
        let fdr = [
            1.7407786385754333e-29,
            7.3700391423612324e-22,
            1.0558631783654333e-17,
            3.4576638739475986e-11,
        ];
        for (i, row) in top.iter().enumerate() {
            assert!(rclose(row.f.unwrap(), f[i]));
            assert!(rclose(row.p_value, p[i]));
            assert!(rclose(row.fdr, fdr[i]));
        }
    }

    #[test]
    fn top_splice_simes_matches_r() {
        let ds = fixture();
        let top = top_splice(&ds, 1, SpliceTest::Simes, 100, 1.0, SpliceSort::P, 0.0).unwrap();
        let ids: Vec<&str> = top.iter().map(|r| r.id.as_str()).collect();
        assert_eq!(ids, ["G5", "G1", "G3", "G2"]);
        let p = [
            4.3519465964385832e-30,
            1.8221109981110334e-20,
            3.0664863233974733e-18,
            3.4576638739475986e-11,
        ];
        let fdr = [
            1.7407786385754333e-29,
            3.6442219962220669e-20,
            4.0886484311966305e-18,
            3.4576638739475986e-11,
        ];
        for (i, row) in top.iter().enumerate() {
            assert!(rclose(row.p_value, p[i]));
            assert!(rclose(row.fdr, fdr[i]));
        }
    }

    #[test]
    fn top_splice_t_matches_r() {
        let ds = fixture();
        let top = top_splice(&ds, 1, SpliceTest::T, 100, 1.0, SpliceSort::P, 0.0).unwrap();
        let ids: Vec<&str> = top.iter().map(|r| r.id.as_str()).collect();
        // Exons 11 and 12 (gene G5) tie to 15 significant digits (p ~= 4.35e-30);
        // their relative order is decided by sub-ULP rounding in `pt`, so only
        // assert they occupy the top two slots. The remaining order is exact.
        let mut head: Vec<&str> = ids[..2].to_vec();
        head.sort_unstable();
        assert_eq!(head, ["11", "12"]);
        assert_eq!(&ids[2..], &["1", "3", "8", "4", "5", "9", "7", "6", "2"]);
        for row in &top[..2] {
            assert!(rclose(row.p_value, 4.3519465964385832e-30));
            assert!(rclose(row.fdr, 2.3935706280412207e-29));
            assert!(rclose(row.logfc.unwrap().abs(), 6.0));
        }
        assert!(rclose(top[2].p_value, 6.0737033270367781e-21));
        assert!(rclose(top[10].p_value, 0.52375685669844319));
    }
}