RustDESeq2

A pure Rust implementation of the DESeq2 algorithm for differential expression analysis of RNA-seq count data. Produces results identical to R DESeq2 while running 8-12x faster with ~90x less memory.

Features

• Complete DESeq2 pipeline: size factor estimation → dispersion estimation → GLM fitting → statistical testing
• Statistical tests: Wald test, Likelihood Ratio Test (LRT)
• Three dispersion fit types: parametric (Gamma GLM), local (LOESS), mean
• LFC shrinkage: normal, apeglm, ashr methods
• Outlier handling: Cook's distance detection and replacement
• Multi-factor designs: batch correction with categorical and continuous covariates
• t-distribution p-values: more conservative testing for small samples (--use-t)
• Transformations: VST (Variance Stabilizing Transformation), rlog (Regularized Log)
• Normalization: median-of-ratios, positive counts, iterative methods
• Convenience functions: fpkm, fpm, collapseReplicates
• 100% significance agreement with R DESeq2 across all validated datasets
• 8-12x faster, ~90x less memory than R DESeq2

Installation

# From crates.io
cargo install rust_deseq2

# From source
git clone https://github.com/necoli1822/rust_deseq2.git
cd rust_deseq2
cargo build --release
# Binary will be at target/release/rust_deseq2

Quick Start

# Basic differential expression analysis
rust_deseq2 run -c counts.csv -m metadata.csv -d condition \
  --numerator treated --denominator control -o results.csv

# With batch correction
rust_deseq2 run -c counts.csv -m metadata.csv -d treatment \
  --covariate batch --numerator drug --denominator placebo

# Likelihood Ratio Test
rust_deseq2 run -c counts.csv -m metadata.csv -d treatment \
  --test lrt --reduced "~1" -o results.csv

# With local dispersion fit and t-distribution
rust_deseq2 run -c counts.csv -m metadata.csv -d condition \
  --numerator treated --denominator control --fit-type local --use-t

# LFC shrinkage (apeglm)
rust_deseq2 run -c counts.csv -m metadata.csv -d condition \
  --numerator treated --denominator control --shrinkage --shrinkage-method apeglm

# Variance Stabilizing Transformation (for PCA, heatmaps)
rust_deseq2 vst -c counts.csv -m metadata.csv -d condition -o vst.tsv --blind

# Regularized log transformation
rust_deseq2 rlog -c counts.csv -m metadata.csv -d condition -o rlog.tsv

# Normalize counts only
rust_deseq2 normalize -c counts.csv -o normalized.tsv

CLI Reference

run — Full DESeq2 analysis

normalize — Normalize counts

vst — Variance Stabilizing Transformation

rlog — Regularized Log Transformation

Input File Format

Count matrix (CSV/TSV, auto-detected)

gene_id,sample1,sample2,sample3,...
GENE001,100,200,150,...
GENE002,50,75,60,...

Metadata (CSV)

sample,condition
sample1,control
sample2,control
sample3,treated

Benchmark

Environment

All benchmarks are median of 5 runs on Linux-native filesystem. R timings are proc.time() internal (exclude package loading overhead).

Execution Time

Memory Usage

R memory is dominated by the R runtime, S4 object overhead, and garbage collector. Rust allocates only the data structures needed for computation.

Accuracy

Validation Against R DESeq2

All 8 validated datasets achieve 100% significance agreement — every gene is classified identically (significant or not at padj < 0.05) between R and Rust.

| Dataset | Genes | Pearson r (log2FC) | Max |padj diff| | Significance Agreement | |---|---:|---:|---:|---:| | A. baumannii GSE151925 | 3,851 | 1.00000000 | ~1e-7 | 3,851 / 3,851 (100%) | | M. tuberculosis GSE100097 | 5,220 | 1.00000000 | ~1e-7 | 5,220 / 5,220 (100%) | | Salmonella GSE46391 | 4,418 | 1.00000000 | ~1e-7 | 4,418 / 4,418 (100%) | | P. aeruginosa GSE55197 | 6,014 | 1.00000000 | 5.00e-7 | 6,014 / 6,014 (100%) | | E. coli GSE220559 | 4,498 | 1.00000000 | 5.01e-7 | 4,498 / 4,498 (100%) | | S. aureus GSE130777 | 2,648 | 1.00000000 | 3.27e-5 | 2,648 / 2,648 (100%) | | S. pneumoniae GSE137447 | 2,280 | 0.99999200 | 3.80e-2 | 2,280 / 2,280 (100%) | | B. subtilis (simulated) | 4,200 | 1.00000000 | 1.21e-5 | 4,200 / 4,200 (100%) |

Across all 33,129 genes tested, there are zero discordant significance calls.

Why Are There Small Numerical Differences?

The minor differences in log2 fold-change values (typically < 1e-06) are inherent to floating-point arithmetic and do not indicate a bug. They arise because:

1. Different linear algebra backends. R DESeq2 calls LAPACK/BLAS routines written in Fortran. RustDESeq2 uses pure Rust implementations (ndarray). Even with identical algorithms, different compilers emit different instruction orderings and SIMD vectorizations.

2. Floating-point non-associativity. IEEE 754 arithmetic does not obey the associative law: (a + b) + c may differ from a + (b + c) by the least significant bit. Each matrix multiplication, QR decomposition, or dot product accumulates these rounding differences.

3. IRLS convergence paths. DESeq2 fits a negative binomial GLM via Iteratively Reweighted Least Squares. When two implementations reach the convergence threshold at slightly different iteration counts, the final coefficient estimates can diverge at the last few decimal places.

These differences are a fundamental property of numerical computing — even recompiling the same code with a different compiler version or on a different CPU architecture can change the last few bits of a floating-point result.

R Function Mapping

Key equivalences between R DESeq2 and RustDESeq2:

Roadmap: Single-Cell RNA-seq Support

RustDESeq2 currently targets bulk RNA-seq. The following plan covers full single-cell support, implementing all scRNA-seq features from R DESeq2 in pure Rust (no R/glmGamPoi dependency).

Phase 1: Parameter Exposure (Easy)

Already-implemented algorithms that just need CLI/API configurability.

Phase 2: Observation Weights (Medium)

Per-gene-per-sample weight matrix for zero-inflation models (zinbwave integration).

Phase 3: glmGamPoi Equivalent (Hard)

Pure Rust reimplementation of the glmGamPoi algorithms, optimized for datasets with many cells and sparse counts.

Implementation Notes

• All features will be implemented in pure Rust with no external R/Python dependencies
• Phase 1 can be released independently as it only requires parameter plumbing
• Phase 2 enables zinbwave-style workflows where weights are pre-computed externally
• Phase 3 is the most complex but provides the largest performance benefit for scRNA-seq
• glmGamPoi and observation weights are mutually exclusive (same as R)

License

MIT

Acknowledgments

This project is an independent Rust reimplementation of the DESeq2 algorithm originally developed by Michael Love, Wolfgang Huber, and Simon Anders.

Love, M.I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. https://doi.org/10.1186/s13059-014-0550-8

The original R implementation is available at https://bioconductor.org/packages/release/bioc/html/DESeq2.html (LGPL >= 3).

The LFC shrinkage methods reference:

• apeglm: Zhu, A., Ibrahim, J.G., Love, M.I. (2019). Heavy-tailed prior distributions for sequence count data. Bioinformatics, 35:2084-2092.
• ashr: Stephens, M. (2017). False discovery rates: a new deal. Biostatistics, 18:275-294.