Nullstrap-DE: Synthetic Calibration for RNA-seq

• The paper introduces Nullstrap-DE as a calibration layer that replaces fixed significance cutoffs with a threshold derived from synthetic null data to control the FDR.
• It preserves the original DESeq2/edgeR framework, maintaining negative binomial modeling, log fold-change estimation, and gene ranking while adjusting significance.
• Simulation studies and real-data applications demonstrate that Nullstrap-DE achieves accurate FDR control with high power and improved biological specificity.

Nullstrap-DE is a synthetic–null-based calibration layer for RNA-seq differential expression analysis that is designed to augment parametric pipelines such as DESeq2 and edgeR without modifying their internal implementations. It retains the same negative binomial generalized linear models, the same log fold-changes and dispersions, and the same gene ranking, but replaces the usual significance cutoff with a threshold calibrated from synthetic null data so as to provide reliable false discovery rate control while preserving power (Jiang et al., 28 Jul 2025).

1. Statistical setting and motivation

Differential expression analysis in bulk RNA-seq starts from a count matrix

$\mathbf{Y} = [Y_{ij}] \in \mathbb{R}^{n \times m},$

where $n$ denotes samples and $m$ denotes genes. For each gene $j$, the inferential target is the null hypothesis of no condition effect on mean expression,

$$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$

Because thousands of genes are tested simultaneously, false discovery rate control is central: $$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$

The framework is motivated by the familiar FDR–power trade-off. DESeq2 and edgeR are described as the dominant tools because they model RNA-seq counts directly and borrow strength across genes. In the general formulation used by Nullstrap-DE, counts are modeled as

$Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$

usually with $F$ taken to be negative binomial, and the mean is written on the log scale as

$$\log(\mu_{ij}) = \log(s_i) + \alpha_j + \mathbf{x}_i^\top \boldsymbol\beta_j + \mathbf{z}_i^\top \boldsymbol\gamma_j. \tag{1}$$

Here $s_i$ is the size factor, $n$0 is the intercept, $n$1 is the treatment design, $n$2 is the treatment effect, $n$3 denotes additional covariates, and $n$4 their coefficients.

The paper characterizes DESeq2 and edgeR as highly powered because they combine negative binomial modeling with shrinkage of gene-wise dispersion estimates. At the same time, it argues that mild deviations from the negative binomial model can lead to FDR inflation, especially in large-sample settings, where small misspecifications become statistically visible across many genes. Non-parametric procedures such as Wilcoxon are presented as more robust for FDR control in some settings, but they often have lower power and do not support covariate adjustment. Nullstrap-DE is introduced to bridge this gap by preserving the modeling flexibility and ranking behavior of DESeq2 and edgeR while calibrating the final decision rule (Jiang et al., 28 Jul 2025).

2. Core construction: synthetic null parallelism

The central idea is to treat the parent differential expression method as a black box that maps counts and design matrices to gene-level statistics, then run that same procedure on both observed data and synthetic null data. The paper describes this as “synthetic null parallelism.”

Let

$n$5

where $n$6 denotes the estimation-and-testing routine supplied by DESeq2 or edgeR. For multi-parameter treatment effects,

$n$7

and in the two-condition case this simplifies to a scalar Wald form,

$n$8

or equivalently

$n$9

In practice, the paper notes that $m$0 is often used because it is monotonically related to the internal test statistic and preserves gene ranking.

Synthetic null generation proceeds gene by gene under the null constraint $m$1. Using fitted size factors, intercepts, covariate effects, and dispersions from the real data, the null mean is

$m$2

or, in Algorithm 1, with resampled size factors $m$3,

$m$4

Synthetic counts are then sampled from

$m$5

typically

$m$6

Once the synthetic matrix $m$7 has been generated, the same DE method is rerun to obtain null statistics $m$8. Nullstrap-DE then estimates the false discovery proportion at threshold $m$9 by comparing exceedances in real and synthetic data: $j$0 and

$j$1

For a target FDR level $j$2, the calibrated threshold is

$j$3

Genes with $j$4 are then called differentially expressed (Jiang et al., 28 Jul 2025).

3. Add-on workflow for DESeq2 and edgeR

Nullstrap-DE is explicitly formulated as an add-on rather than a replacement. The parent pipeline is run on the observed data first, yielding the standard fitted quantities: size factors $j$5, dispersions $j$6, fitted coefficients $j$7, $j$8, $j$9, and gene-level statistics or p-values. The framework then generates a single synthetic null dataset by erasing the treatment effect while preserving gene-specific intercepts and dispersions, covariate effects, and the empirical distribution of size factors.

The implementation described in the paper fixes the synthetic analysis as closely as possible to the original one. DESeq2 or edgeR is rerun on $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$0, with the size factors fixed at $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$1 and dispersions fixed at $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$2. This makes the real and synthetic statistics directly comparable, because the only intended difference is that the synthetic data arise under the null.

This workflow preserves several features of the parent method. It preserves gene ranking, because the observed statistics $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$3 are unchanged. It preserves effect-size estimation, because log fold-changes are still estimated by DESeq2 or edgeR on the real data. It also inherits the ability of those methods to handle multiple conditions, batch effects, and other covariates through $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$4 and $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$5. In that sense, Nullstrap-DE changes only the final calibration of significance, not the inferential model fitted to the observed counts (Jiang et al., 28 Jul 2025).

4. Asymptotic guarantees

The theoretical analysis is developed primarily for the two-condition setting, with scalar treatment effects $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$6. Let $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$7 be the set of truly null genes and $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$8 the set of truly differentially expressed genes, with $$H_{0,j}: \text{no effect of condition on mean expression}, \quad H_{1,j}: \text{nonzero effect}.$$9. The theory rests on three assumptions: uniform estimation accuracy for real and synthetic fitted effects, regularity of the null statistic distribution around the data-driven threshold, and a beta-min condition ensuring separation of non-null effects from estimation noise.

Under Assumptions 1 and 2, the paper gives an asymptotic FDR bound for the selected set $$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$0: $$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$1 with

$$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$2

If Assumption 3 also holds, then the power satisfies

$$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$3

The interpretation given in the paper is that FDR is controlled asymptotically at approximately the target level $$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$4, with corrections determined by the estimation error scale $$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$5 and vanishing terms in $$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$6. For negative binomial generalized linear models, the supplement shows that maximum likelihood estimation yields

$$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$7

under regularity conditions, implying that when $$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$8, Nullstrap-DE applied to NB-GLMs achieves asymptotic FDR control and power tending to 1.

For finite samples, the paper recommends a conservative safety margin. Instead of enforcing $$\mathrm{FDR} = \mathbb{E}\left[\frac{\#\text{ false positives}}{\max(\#\text{ calls},1)}\right].$$9, one may define the operational cutoff by

$Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$0

This recommendation reflects the rate $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$1 and is intended to improve finite-sample FDR control (Jiang et al., 28 Jul 2025).

5. Simulation studies and real-data applications

The simulation study compares DESeq2, edgeR, Wilcoxon rank-sum on raw counts, Wilcoxon on normalized counts, Nullstrap-DESeq2, and Nullstrap-edgeR. In the first setting there are two conditions, no extra covariates, $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$2 genes, differential expression proportions in $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$3, balanced sample sizes $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$4, fold changes in $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$5, target FDR levels $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$6, and size factors sampled from $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$7. The reported pattern is that Nullstrap-DESeq2 and Nullstrap-edgeR track the nominal FDR well across settings while retaining power close to DESeq2 and edgeR; by contrast, DESeq2 and edgeR show clear FDR inflation as $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$8, fold change, or target $Y_{ij} \sim F(\cdot \mid \mu_{ij}, \phi_j),$9 increase, and Wilcoxon methods are substantially less powerful, especially at small $F$0.

The second simulation setting introduces confounding covariates: an additional binary covariate correlated with treatment, $F$1 of genes with nonzero covariate effects $F$2, and fold changes in $F$3. Here the paper reports that Nullstrap-DESeq2 and Nullstrap-edgeR again maintain reliable FDR control, whereas DESeq2 and edgeR exhibit substantial FDR inflation and Wilcoxon-based methods fail because they cannot adjust for the confounder. Across both settings, the main conclusion is that Nullstrap-DE achieves the combination of calibrated FDR and power close to that of the parent parametric methods (Jiang et al., 28 Jul 2025).

Three real-data applications are used to support the method. In the monocyte negative-control experiment, RNA-seq data from 34 human samples and 52,376 genes are permuted 1000 times so that no true differential expression is expected. DESeq2 and edgeR often detect many genes on these permuted datasets, Wilcoxon on normalized counts detects fewer but still non-negligible false discoveries, whereas Nullstrap-DESeq2 and Nullstrap-edgeR detect essentially zero differentially expressed genes in almost all permutations.

Using the true labels for classical and non-classical monocytes, the paper reports that Nullstrap-DESeq2 and Nullstrap-edgeR call fewer genes at FDR 0.05 than DESeq2, edgeR, or Wilcoxon. The retained genes show clearer expression separation in heatmaps and stronger enrichment for immune-related Gene Ontology terms such as “leukocyte migration” (GO:0050900) and “positive regulation of cytokine production” (GO:0001819). Genes called by edgeR but not by Nullstrap-edgeR are described as exhibiting weaker or noisier patterns and more generic GO enrichment.

A third application uses RNA-seq from airway smooth muscle cells, with 4 untreated and 4 dexamethasone-treated samples and 63,677 genes. At FDR 0.05, Nullstrap-DESeq2 and Nullstrap-edgeR call fewer genes than DESeq2 and edgeR, while Wilcoxon finds no genes at all. The retained genes are reported to produce strong separation between treated and untreated samples, and pathway enrichment emphasizes calcium signaling, cytokine–cytokine receptor interaction, and cAMP signaling, all described as directly related to glucocorticoid biology and airway smooth muscle function (Jiang et al., 28 Jul 2025).

6. Scope, limitations, and extensions

The main strengths claimed for Nullstrap-DE are robust FDR control, high power, compatibility with the full DESeq2 and edgeR workflow, support for covariates and complex designs, and improved biological specificity in downstream interpretation. Because the framework leaves the observed-data fit unchanged, effect sizes, dispersions, and the ranking of genes are preserved. The practical output differs only in the significance threshold and hence in the final set of declared differentially expressed genes.

The limitations are also explicit. Nullstrap-DE requires running DESeq2 or edgeR twice—once on the real data and once on the synthetic null data—so computational cost is roughly doubled, and would increase further if multiple synthetic null datasets were used. The theoretical guarantees depend on reasonably accurate estimation under the parent model; the framework is presented as robust to mild misspecification, not arbitrary failure of the underlying generalized linear model. It also cannot repair every weakness of the parent procedure: unstable dispersion estimates or mis-specified designs remain problematic. Finally, the finite-sample safety margin can make the procedure conservative.

The paper positions Nullstrap-DE as a general framework rather than a method tied only to DESeq2 and edgeR. A stated implication is that the same synthetic-null calibration strategy could be extended to other differential expression methods and to other high-dimensional omics settings, provided there is a clear null generative model and a well-defined gene-level test statistic. Proposed directions include limma-voom, generalized linear mixed models, count-based proteomics, ChIP-seq, ATAC-seq, and single-cell RNA-seq. In that broader view, Nullstrap-DE is presented as an instance of synthetic null parallelism for calibrating FDR in complex model-based inference while preserving the practical advantages of the underlying analysis pipeline (Jiang et al., 28 Jul 2025).