Differential analysis of genomics count data with edge*

This paper introduces edgePython, a Python port of the widely used edgeR package that extends its capabilities with a negative binomial-gamma mixed model and empirical Bayes shrinkage to facilitate differential expression analysis of single-cell genomics data within the Python ecosystem.

The Gist

Imagine you are a detective trying to solve a mystery: Which genes are acting differently in different situations? Maybe you want to know which genes turn on when a jellyfish eats, or which ones change when a mouse is pregnant versus nursing.

For a long time, the best detective tool for this job was a software package called edgeR. It was built by experts in the R programming language and was the gold standard for analyzing "count data" (basically, tallying up how many times a gene was seen).

However, the world of biology has changed. Today, most scientists work with Python because it's the language of choice for analyzing single-cell data (looking at genes inside individual cells, not just a soup of cells). But edgeR only spoke R. This created a language barrier: scientists using Python had to translate their data back and forth to R, which was slow, clunky, and prone to errors.

This paper introduces edgePython, a new tool that solves this problem. Here is the breakdown in simple terms:

1. The Great Translation (The Port)

Think of edgeR as a famous, highly detailed cookbook written in French (R). The scientists wanted to translate this exact cookbook into English (Python) so everyone could use it without needing a French dictionary.

• What they did: They didn't just guess; they meticulously translated every recipe (function) from the original French book into English.
• The Result: edgePython. It does exactly what the original edgeR does. If you run the same experiment in both, you get the exact same answer. It's like having a perfect twin of the original tool that speaks the language the modern community uses.

2. The New Superpower (Single-Cell Analysis)

The original edgeR was great at looking at a "smoothie" of cells (bulk data), but it struggled when looking at individual cells (single-cell data).

• The Problem: In single-cell data, you have many cells from the same person (or animal). These cells are related, like siblings. If you treat them all as completely independent strangers, you get false alarms (thinking a gene changed when it didn't).
• The Solution: The authors added a new "smart filter" to edgePython. They used a statistical model (a Negative Binomial–Gamma mixed model) that understands family relationships. It knows that cells from the same mouse are related and adjusts the math accordingly.
• The "Shrinkage" Trick: Imagine you are trying to guess the height of a tree, but you only have a tiny, shaky ruler. Your guess might be wild. But if you look at 1,000 other trees and see a pattern, you can "shrink" your wild guess toward the average to make it more reliable. edgePython does this mathematically to make sure the results are stable, even when you don't have many cells to study.

3. The Speed Boost

The original edgeR had some parts written in a very fast, low-level language called C. The new edgePython is written in Python, which is usually slower. However, the authors used a special "turbocharger" (called Numba) that compiles the Python code to run at the speed of C.

• The Result: For complex single-cell analyses, edgePython is actually faster than the original R version.

4. The AI Assistant

Here is the most futuristic part of the paper. The authors didn't just write the code; they used an AI (Large Language Model) to help translate the complex statistical code from R to Python.

• The Analogy: It's like having a master translator who can read a 20-year-old, handwritten technical manual and instantly rewrite it in a modern language, fixing typos and optimizing the flow.
• The Implication: This suggests that in the future, we might not need to spend years learning to code complex tools. We might just ask an AI to "port" a tool to a new language, making advanced science accessible to anyone who can speak English.

Summary

edgePython is a bridge.

1. It bridges the gap between R (the old guard) and Python (the new guard).
2. It bridges the gap between simple gene analysis and complex single-cell analysis.
3. It bridges the gap between human coding and AI-assisted development.

It allows scientists to use the most powerful statistical tools in the world without getting stuck in a language barrier, making it easier to discover how life works at the cellular level.

Technical Summary

1. Problem Statement

The edgeR Bioconductor package is a gold-standard tool for differential expression analysis of count-based genomics data, utilizing a negative binomial framework and empirical Bayes moderation. However, it faces two critical limitations in the modern genomics landscape:

• Ecosystem Fragmentation: edgeR is implemented exclusively in R. The single-cell genomics field has largely shifted to a Python-centric ecosystem (dominated by Scanpy, AnnData, and scverse). This forces researchers to perform complex data conversions between R and Python or rely on fragile inter-language bridges, creating a significant barrier to integrating edgeR's robust statistical methods into single-cell workflows.
• Statistical Gaps in Single-Cell Analysis: While edgeR handles quantification uncertainty (e.g., via sleuth-like adaptations), it lacks a native framework for multi-subject single-cell differential expression. Standard edgeR treats cells as independent replicates, ignoring the hierarchical structure of data (cells nested within subjects). This leads to inflated false positive rates. Existing mixed-model solutions (like NEBULA) do not yet incorporate edgeR's powerful empirical Bayes shrinkage for cell-level dispersion.

2. Methodology

The authors present edgePython, a comprehensive Python port of edgeR (version 4.8.2) that extends its capabilities for single-cell analysis.

A. Software Porting Strategy

• Translation: The core R source code (including dependencies in limma, statmod, and locfit) was translated into Python.
• AI-Assisted Development: The porting was executed using Large Language Models (Claude Opus 4.5 and 4.6). The process involved feeding R source code, documentation, and test cases to the LLM, followed by iterative debugging to resolve numerical discrepancies (e.g., optimization defaults, floating-point edge cases).
• Architecture:
  • Data Structures: Python dictionaries mirror edgeR's S3 list-based objects.
  • Computation: Operations utilize NumPy arrays and SciPy sparse matrices.
  • Performance: Performance-critical inner loops (specifically likelihood and gradient evaluations for the mixed model) are compiled using Numba (@njit) to achieve speeds comparable to or faster than the R implementation.
  • Coverage: The port includes 86 public functions across 24 modules, covering normalization (TMM, RLE), dispersion estimation, GLM fitting, all four hypothesis testing frameworks (Exact, LRT, QL F-test, TREAT), and gene set testing (camera, fry, roast, etc.).

B. Statistical Extension: Negative Binomial–Gamma Mixed Model

To address the multi-subject single-cell challenge, the authors implemented a Negative Binomial–Gamma mixed model (following the NEBULA-LN approach):

• Model Structure:
  $$Y_{gij} | b_{gj} \sim NB(\mu_{gij} e^{b_{gj}}, \phi_g)$$
  $$b_{gj} \sim N(0, \sigma^2_g)$$
  Where $Y_{gij}$ is the count for gene $g$, cell $i$, subject $j$; $b_{gj}$ is a subject-level random effect; and $\phi_g$ is the cell-level overdispersion.
• Empirical Bayes Shrinkage (Novel Contribution):
  • Standard NEBULA estimates cell-level dispersion ($\phi_g$) via Maximum Likelihood Estimation (MLE), which can be noisy with modest cell counts.
  • edgePython introduces empirical Bayes shrinkage of these cell-level dispersions toward an abundance-dependent prior trend.
  • This utilizes the squeezeVar function (originally from edgeR's quasi-likelihood workflow) to "borrow strength" across genes, stabilizing inference in small-sample regimes. Neither the original edgeR nor NEBULA currently performs this specific shrinkage step in a mixed-model setting.

C. Integration and Interoperability

• AnnData/Seurat Support: The tool provides bidirectional conversion between its internal DGEList and the AnnData (Scanpy) and Seurat formats.
• Input Handling: Direct import of kallisto HDF5 files (including bootstrap samples for quantification uncertainty), a feature previously missing in edgeR.
• AI Agents: A Model Context Protocol (MCP) server is included, allowing AI agents to execute the full analysis pipeline via natural language.

3. Key Results

The authors validated edgePython through extensive benchmarking against the original R edgeR and NEBULA packages.

• Numerical Concordance:
  • Bulk RNA-seq: On the HOXA1 (knockdown) and GSE60450 (mouse mammary gland) datasets, edgePython outputs (normalization factors, BCV, log-fold-changes, p-values) matched edgeR R outputs with a relative tolerance of $<10^{-3}$.
  • Gene Set Testing: Methods like camera and fry reproduced results across various configurations.
  • Specialized Analyses: Differential Transcript Usage (DTU) and scaled analysis workflows showed perfect agreement.
• Single-Cell Performance:
  • Validation: The NEBULA-LN port matched R's NEBULA package in coefficient estimates and p-values.
  • Shrinkage Efficacy: In a subsampled dataset (30 cells), the empirical Bayes shrinkage significantly stabilized noisy MLE dispersions, pulling them toward the prior trend and reducing false positives.
  • Biological Discovery: Applied to Clytia hemisphaerica (jellyfish) data, the method identified 689 differentially expressed genes (FDR < 0.05), revealing biological signals (e.g., upregulation of metabolic genes in fed states) that differed from previous analyses using other frameworks.
• Runtime:
  • For bulk RNA-seq, performance is comparable to R.
  • For single-cell mixed models, the Numba-compiled Python implementation is substantially faster than R's NEBULA, with the speed advantage increasing with dataset size.

4. Significance and Impact

• Bridging the Ecosystem Gap: edgePython removes the barrier between the statistical rigor of edgeR and the single-cell dominance of the Python ecosystem, enabling seamless integration with Scanpy and AnnData.
• Methodological Advancement: It is the first tool to combine mixed-models for multi-subject single-cell data with empirical Bayes shrinkage of cell-level dispersion, offering a more robust statistical framework for complex experimental designs.
• Demonstration of AI in Software Engineering: The project serves as a proof-of-concept that LLMs can successfully port complex, C-heavy statistical packages (edgeR is ~39% C code) across languages. The author notes the project was completed in one week, suggesting a paradigm shift in how scientific software is developed and maintained.
• Future-Proofing: The "edge*" concept implies that the barrier to porting to other languages (e.g., Rust, CUDA) is now low, potentially accelerating the development of high-performance genomics tools.

5. Limitations

• Excluded Functions: processAmplicons (CRISPR barcode processing) was not ported to maintain a clean separation of concerns (starting from count matrices).
• Annotation: Pathway analysis (goana, kegga) delegates to the external g:Profiler API rather than using local Bioconductor databases.
• Approximations: Only the NEBULA-LN (first-order Laplace) approximation is implemented; the higher-order NEBULA-HL is not yet available (though the authors argue the new shrinkage method mitigates the need for it in small samples).
• Single-Matrix Constraint: Like standard edgeR, it operates on a single count matrix, whereas emerging methods for spliced/unspliced RNA (e.g., Monod) utilize joint distributions of multimodal counts.

In summary, edgePython represents a major step forward in computational genomics, democratizing access to state-of-the-art differential expression statistics for Python users while introducing novel statistical refinements for single-cell data analysis.