Troubleshooting Differential Expression in Single-Cell Data: Why Are My Results Not Reproducible?

By Dr. Zubair Khalid, DVM, MS, PhD ·

Troubleshooting Differential Expression in Single-Cell Data: Why Are My Results Not Reproducible?

Key Takeaways

  • Non-reproducible differential expression (DE) results in single-cell RNA sequencing (scRNA-seq) are primarily driven by identifiable sources of variability including batch effects, cell-level filtering thresholds, pseudobulk construction choices, normalization methods, and the statistical model used for hypothesis testing.
  • Batch effects are a critical confounder; visualizing data with PCA or UMAP colored by batch before integration is the first diagnostic check, as DE genes clustering by batch instead of biological condition indicates a significant issue.
  • Cell-level filtering thresholds directly impact cell population composition; re-running filtering with a range of thresholds and comparing gene counts and cell retention is essential to assess sensitivity and potential divergence in DE gene lists.
  • The choice of statistical model is paramount, with methods that ignore biological replicate variation being prone to hundreds of false positives; comparing DE results from models that explicitly account for donor or sample variation against those treating each cell as independent is a crucial diagnostic step.
  • Data sparsity and sequencing depth influence method selection; for low-depth data, methods robust to sparsity like limmatrend or Wilcoxon tests are preferred over zero-inflation models, while higher depth allows for more complex modeling approaches.
  • Pseudobulk construction versus cell-level analysis presents a trade-off between statistical power and reproducibility; comparing DE gene lists from both approaches can reveal discrepancies stemming from the statistical treatment of biological replicates.

Differential expression (DE) analysis in single-cell RNA sequencing (scRNA-seq) and single-nucleus RNA sequencing (snRNA-seq) frequently produces inconsistent results when the same biological question is analyzed across runs, datasets, or software versions. The primary causes are identifiable sources of variability: batch effects, filtering thresholds, pseudobulk construction choices, normalization methods, and the statistical model used for hypothesis testing. This article provides a diagnostic framework for researchers who observe non-reproducible DE results and need concrete checks to locate the source of divergence. The scope covers data inputs, workflow decisions, quality controls, and interpretation limits relevant to biology students, laboratory professionals, and life-science practitioners working with single-cell transcriptomic data.

At a Glance: Common Sources of Non-Reproducible Differential Expression

The table below summarizes the most frequent sources of variability in single-cell DE analysis, the typical symptom observed by the researcher, and the diagnostic check that should be applied first.

Source of VariabilityTypical Symptom in ResultsFirst Diagnostic Check
Batch effects between samples or sequencing runsDE genes cluster by experimental batch instead of biological conditionVisualize data with PCA or UMAP colored by batch before integration
Cell-level filtering thresholdsNumber of DE genes changes dramatically when minimum gene count or mitochondrial percentage thresholds are adjustedRe-run the filtering step with a range of thresholds and compare gene counts and cell retention
Pseudobulk construction methodDE results differ between pseudobulk and cell-level analysisCompare DE gene lists from both approaches and examine whether biological replicates are preserved
Normalization choiceRank order of DE genes shifts across normalization methodsRun two normalization approaches and compare the overlap of top-ranked genes
Statistical model selectionMethods that ignore biological replicate variation produce hundreds of false positivesUse methods that model donor or sample variation and compare against methods that treat each cell as independent
Data integration strategyBatch-corrected data changes DE outcomes compared to uncorrected dataTest DE on both corrected and uncorrected data and compare results for sparse versus deep datasets

Understanding Why Single-Cell Differential Expression Results Vary

Single-cell transcriptomic data contain technical noise and intrinsic biological variability that complicate the detection of differential expression signatures. A probabilistic model of expression-magnitude distortions typical of single-cell RNA-sequencing measurements was developed to enable detection of differential expression signatures and identification of subpopulations of cells in a way that is more tolerant of noise [<a href="#ref-1">1</a>]. This early recognition of the noise problem established that DE analysis in single-cell data requires statistical approaches distinct from bulk RNA-seq.

The core challenge is that each cell provides a sparse measurement of the transcriptome. Most genes are not detected in most cells, and the depth of sequencing varies substantially across cells within the same experiment. When a researcher observes non-reproducible DE results, the first question should be whether the variability originates from the biological material, the library preparation, the computational pipeline, or the statistical test.

A benchmark of 46 workflows for DE analysis of single-cell data with multiple batches showed that batch effects, sequencing depth, and data sparsity substantially impact performance [<a href="#ref-2">2</a>]. The same study found that the use of batch-corrected data rarely improves the analysis for sparse data, whereas batch covariate modeling improves the analysis for substantial batch effects [<a href="#ref-2">2</a>]. This finding has direct practical implications: the choice of whether to correct for batch effects before or during DE analysis depends on the depth of the data and the severity of the batch structure.

The relative performance of DE methods is contingent on their ability to account for variation between biological replicates. Methods that ignore this inevitable variation are biased and prone to false discoveries, and the most widely used methods can discover hundreds of differentially expressed genes in the absence of biological differences [<a href="#ref-3">3</a>]. This observation was demonstrated by exposing true and false discoveries of differentially expressed genes in the injured mouse spinal cord [<a href="#ref-3">3</a>]. For a researcher seeing inconsistent results, the first suspicion should be that the statistical method treats each cell as an independent observation instead of accounting for the fact that cells from the same biological sample are correlated.

Core Principles of Reproducible Single-Cell Differential Expression

Biological Replicates Are the Foundation of Reproducibility

The single most important principle for reproducible DE analysis is that biological replicates must be preserved in the statistical model. Single-cell data provide a means to dissect the composition of complex tissues and specialized cellular environments, but the analysis of such measurements is complicated by high levels of technical noise and intrinsic biological variability [<a href="#ref-1">1</a>]. When cells from the same donor or the same animal are treated as independent observations, the effective sample size is inflated and the statistical test becomes anti-conservative.

The practical consequence is that DE genes identified without accounting for biological replicates are often not reproducible when the experiment is repeated with new biological samples. A reproducible Seurat-based protocol for analyzing PBMC CD4+ T-cell single-cell RNA sequencing data across malaria reinfection timepoints demonstrates the importance of a unified computational framework that includes standardized preprocessing, integration, clustering, and downstream transcriptomic analyses [<a href="#ref-4">4</a>]. The protocol enables systematic computation of module scores for predefined immune and CD4+ T-cell programs, identification of cluster-specific marker genes, timepoint-resolved differential expression analysis, and downstream Gene Ontology and KEGG pathway enrichment [<a href="#ref-4">4</a>]. The emphasis on a standardized workflow reflects the reality that small changes in any step can propagate through the analysis and alter DE results.

Data Sparsity and Zero Inflation Affect Method Choice

Single-cell data contain excessive zeros, and the handling of these zeros is a major challenge in DE analysis. A statistical framework that leverages UMI counts and zero proportions within a generalized Poisson or binomial mixed-effects model was proposed to account for batch effects and within-sample variation [<a href="#ref-5">5</a>]. The authors of that framework dissected four major challenges in single-cell DE analysis: excessive zeros, normalization, donor effects, and cumulative biases [<a href="#ref-5">5</a>]. These challenges underscore the limitations and conceptual pitfalls in existing workflows.

For low-depth data, single-cell techniques based on zero-inflation models deteriorate performance, whereas the analysis of uncorrected data using limmatrend, Wilcoxon test, and fixed effects models performs well [<a href="#ref-2">2</a>]. This finding suggests that the choice of statistical method should be guided by the depth of the sequencing data. A researcher working with shallow sequencing should avoid methods that model zero inflation explicitly and instead use methods that are robust to sparsity.

Normalization Choices Change Gene Rankings

Normalization is a source of variability that is often underestimated. The framework proposed in the Genome Biology study preserves biologically meaningful signals and offers improved performance in detecting differentially expressed genes by using absolute RNA expression instead of relative abundance [<a href="#ref-5">5</a>]. This paradigm shift challenges existing workflows and highlights the need for careful consideration of normalization procedures [<a href="#ref-5">5</a>].

When a researcher observes different DE results across runs, one of the first checks should be whether the normalization method changed between runs. Library size normalization, scran normalization, and methods that model absolute expression will produce different gene rankings, particularly for genes with moderate expression levels.

Practical Workflow for Diagnosing Non-Reproducible Results

Step 1: Document the Exact Pipeline Version and Parameters

Reproducibility begins with documentation. The Galaxy Training Network provides accessible workflow training and analysis tutorials with reproducibility context [<a href="#ref-6">6</a>]. The nf-core documentation describes community pipeline standards, usage, configuration, and reproducible workflow context [<a href="#ref-7">7</a>]. Both resources emphasize that pipeline version and parameter changes are common sources of apparent non-reproducibility.

Before investigating biological or statistical causes, confirm that the same software versions were used across runs. Changes in alignment tools, quantification methods, or DE packages can alter results even when the input data are identical. The Bioconductor project provides official package, workflow, installation, and reproducible genomic-analysis documentation [<a href="#ref-8">8</a>]. Version tracking for all packages used in the analysis is a prerequisite for diagnosing reproducibility problems.

Step 2: Compare Cell-Level Quality Control Metrics Across Runs

Cell-level quality control thresholds directly affect which cells enter the DE analysis. The isolation and preparation of cells for RNA analysis can introduce variability that is unrelated to biology. Protocols for isolating viable cells from complex tissues emphasize the importance of troubleshooting guidelines for efficient isolation [<a href="#ref-9">9</a>]. For example, a protocol for isolating viable adipocytes and stromal vascular fraction from human visceral adipose tissue describes steps for washing, mincing, and digesting tissue to generate a single-cell suspension [<a href="#ref-9">9</a>]. The RNA extraction protocol ensures a high yield of total RNA for downstream expression assays [<a href="#ref-9">9</a>].

For single-cell data, the key quality metrics are the number of genes detected per cell, the total number of UMI counts per cell, and the percentage of mitochondrial reads. When DE results differ across runs, compare the distributions of these metrics. If the filtering thresholds removed different proportions of cells, the cell populations being compared may differ.

Step 3: Examine Batch Structure Before Integration

Before applying any data integration method, visualize the data colored by batch, donor, or sequencing run. The NCBI provides official descriptions of databases, search systems, sequence resources, and analysis services that can be used to access public datasets for comparison [<a href="#ref-10">10</a>]. If cells cluster by batch instead of by biological condition, batch effects are present and must be addressed.

The decision to integrate data should be guided by the downstream analysis goal. A benchmark of integration strategies for DE analysis found that the use of batch-corrected data rarely improves the analysis for sparse data, whereas batch covariate modeling improves the analysis for substantial batch effects [<a href="#ref-2">2</a>]. This means that integration is not always beneficial for DE analysis. For sparse data, running DE on uncorrected data with batch included as a covariate in the model may produce more reproducible results.

Step 4: Test the Sensitivity of Results to Filtering Thresholds

Filtering thresholds are a common source of non-reproducibility because small changes in thresholds can change the cell population composition. A critical assessment of differential expression for bulk RNA-Seq projects found that finding the right balance of quality and quantity can be important, and it is essential that project quality does not drop below the level where important main conclusions are missed or misstated [<a href="#ref-11">11</a>]. The same study used knock-out and over-expression studies as a simplification to test recovery of a known causal gene in RNA-Seq cell line experiments [<a href="#ref-11">11</a>].

For single-cell data, run the filtering step with a range of thresholds. For example, vary the minimum number of genes detected per cell from 200 to 1000 and the maximum mitochondrial percentage from 5% to 20%. Record the number of cells retained and the number of DE genes identified at each threshold combination. If the DE gene list changes substantially across reasonable threshold ranges, the results are sensitive to filtering choices and should be interpreted with caution.

Step 5: Compare Pseudobulk and Cell-Level Analysis

Pseudobulk construction aggregates cells from the same biological sample into a single expression profile, which preserves the biological replicate structure. A reproducible protocol for single-cell analysis emphasizes the importance of standardized preprocessing and downstream transcriptomic analyses within a unified computational framework [<a href="#ref-4">4</a>]. The choice between pseudobulk and cell-level analysis is one of the most consequential decisions in single-cell DE analysis.

When results are not reproducible, compare the DE gene lists obtained from pseudobulk analysis and cell-level analysis. If the two approaches produce substantially different gene lists, the discrepancy likely stems from the statistical treatment of biological replicates. Pseudobulk analysis is generally more conservative and produces results that are more reproducible across independent experiments.

Step 6: Evaluate the Statistical Model Choice

The statistical model used for DE testing has a profound impact on results. Methods that ignore variation between biological replicates are biased and prone to false discoveries [<a href="#ref-3">3</a>]. The most widely used methods can discover hundreds of differentially expressed genes in the absence of biological differences [<a href="#ref-3">3</a>].

A benchmark of DE methods found that for single-end RNA-Seq reads aligned with STAR and quantified with htseq-count, there is potential value in testing the use of the generalized linear model implementation of edgeR with robust dispersion estimation more frequently for either single-variate or multi-variate two-group comparisons [<a href="#ref-11">11</a>]. The same study noted that when considering a limited number of patient sample comparisons with larger sample size, there might be some decreased variability between methods [<a href="#ref-11">11</a>].

For single-cell data, the choice of method should account for the experimental design. Methods that model donor effects and within-sample variation are more adaptable to diverse experimental designs [<a href="#ref-5">5</a>]. A method of moments framework for DE analysis of scRNA-seq data was introduced for robust and efficient differential analysis of mean expression, variability, and gene correlation from scRNA-seq data scalable to millions of cells and thousands of samples [<a href="#ref-12">12</a>]. This framework identified more significant and reproducible differences in mean expression compared to existing methods [<a href="#ref-12">12</a>].

Options and Tradeoffs in Single-Cell Differential Expression Workflows

Cell-Level Analysis Versus Pseudobulk Analysis

Cell-level analysis treats each cell as an observation and uses methods such as the Wilcoxon test, limma-trend, or zero-inflated models. Pseudobulk analysis aggregates cells within each biological sample and uses bulk RNA-seq DE methods such as edgeR or DESeq2. The tradeoff is between statistical power and reproducibility.

Cell-level analysis has higher power to detect small expression differences because the sample size is the number of cells, which can be thousands. However, this power is illusory when cells from the same biological sample are correlated. The effective sample size is the number of biological replicates, not the number of cells. Pseudobulk analysis correctly treats the biological replicate as the unit of analysis and produces results that are more reproducible across experiments.

A benchmark of 46 workflows found that DE analysis for a specific cell type outperforms that of large-scale bulk sample data in prioritizing disease-related genes [<a href="#ref-2">2</a>]. This finding supports the use of cell-type-specific DE analysis but does not resolve the question of whether cell-level or pseudobulk methods should be preferred.

Data Integration Before DE Versus Batch Covariate Modeling

Data integration methods such as Harmony, Seurat CCA, and scVI correct for batch effects by aligning cells across batches. An alternative approach is to include batch as a covariate in the DE model without correcting the data. The tradeoff depends on the severity of the batch effect and the sparsity of the data.

The benchmark of integration strategies found that the use of batch-corrected data rarely improves the analysis for sparse data, whereas batch covariate modeling improves the analysis for substantial batch effects [<a href="#ref-2">2</a>]. For low-depth data, single-cell techniques based on zero-inflation models deteriorate performance, whereas the analysis of uncorrected data using limmatrend, Wilcoxon test, and fixed effects models performs well [<a href="#ref-2">2</a>].

The practical implication is that researchers should test both approaches and compare the results. If the DE gene lists are substantially different between corrected and uncorrected data, the batch effect is influencing the results and the choice of approach should be documented and justified.

Normalization Methods

Normalization methods differ in how they scale expression values across cells. Library size normalization divides each cell by its total UMI count. More sophisticated methods account for composition differences or use spike-in controls. The choice of normalization affects DE results because it changes the relative expression values.

A statistical framework that uses absolute RNA expression instead of relative abundance improves sensitivity, reduces false discoveries, and enhances biological interpretability [<a href="#ref-5">5</a>]. This approach challenges existing workflows and highlights the need for careful consideration of normalization procedures [<a href="#ref-5">5</a>].

When troubleshooting non-reproducible results, run the analysis with two different normalization methods and compare the overlap of DE genes. If the overlap is low, the results are sensitive to normalization and should be interpreted with caution.

Records and Measurements for Reproducible Analysis

What to Record for Each Analysis Run

Reproducible analysis requires detailed records of every step. The Carpentries provides foundational computing, data, shell, Git, and programming training context [<a href="#ref-13">13</a>]. Version control is a core component of reproducible analysis. The following records should be maintained for each DE analysis run:

  1. Software versions for all packages used in the pipeline, including the operating system and programming language version
  2. Exact parameters for each step, including filtering thresholds, normalization method, integration method, and DE model formula
  3. Number of cells retained at each filtering step and the number of cells per biological sample
  4. Number of genes detected per cell and the distribution of UMI counts
  5. Batch structure and the assignment of cells to batches
  6. DE method and model formula, including covariates included in the model
  7. Random seed if any stochastic steps are used in the pipeline

The nf-core documentation describes community pipeline standards for reproducible workflow context [<a href="#ref-7">7</a>]. These standards include version tracking and parameter documentation. The Galaxy Training Network provides accessible workflow training and analysis tutorials with reproducibility context [<a href="#ref-6">6</a>].

Quality Metrics to Track Across Runs

The following quality metrics should be tracked across runs to diagnose non-reproducibility:

  1. Median genes detected per cell
  2. Median UMI counts per cell
  3. Percentage of cells passing filtering thresholds
  4. Number of clusters identified and the cell type composition of each cluster
  5. Proportion of cells in each cluster that come from each batch or donor
  6. Number of DE genes identified at a fixed significance threshold
  7. Overlap of DE gene lists between runs

If any of these metrics change substantially between runs, the source of the change should be identified before the DE results are compared.

Common Failure Patterns in Single-Cell Differential Expression

Failure Pattern 1: Treating Cells as Independent Observations

The most common cause of non-reproducible DE results is the use of statistical methods that treat each cell as an independent observation. This approach ignores the correlation between cells from the same biological sample. Methods that ignore variation between biological replicates are biased and prone to false discoveries [<a href="#ref-3">3</a>]. The most widely used methods can discover hundreds of differentially expressed genes in the absence of biological differences [<a href="#ref-3">3</a>].

The diagnostic check is to compare the number of DE genes identified by a cell-level method with the number identified by a pseudobulk method. If the cell-level method identifies substantially more DE genes, the difference is likely due to false discoveries from ignoring biological replicate variation.

Failure Pattern 2: Inconsistent Filtering Thresholds

Filtering thresholds that differ across runs change the cell population composition. If one run uses a minimum of 500 genes per cell and another uses 1000 genes per cell, the cell populations being compared are different. The DE results will differ because the cell populations are not the same.

The diagnostic check is to run the filtering step with a range of thresholds and compare the number of cells retained and the number of DE genes identified. If the results are sensitive to threshold choices, the thresholds should be specified in advance and documented.

Failure Pattern 3: Batch Effects Confounded with Biological Condition

When batch effects are confounded with the biological condition of interest, DE results are not reproducible because the batch effect is not separable from the biological effect. For example, if all control samples are processed in one batch and all treated samples in another batch, any DE genes may reflect batch differences instead of treatment effects.

The diagnostic check is to examine the experimental design. If batch and condition are confounded, the experiment cannot produce reproducible DE results without additional data. The EMBL-EBI Training provides bioinformatics learning pathways and data-resource training that cover experimental design considerations [<a href="#ref-14">14</a>].

Failure Pattern 4: Integration Method Alters the DE Signal

Data integration methods can remove biological signal along with batch effects. The benchmark of integration strategies found that the use of batch-corrected data rarely improves the analysis for sparse data [<a href="#ref-2">2</a>]. If integration removes the biological signal of interest, DE results will not be reproducible.

The diagnostic check is to compare DE results from corrected and uncorrected data. If the DE gene lists are substantially different, the integration method is altering the signal and the choice of approach should be reconsidered.

Failure Pattern 5: Normalization Method Changes Gene Rankings

Normalization methods that produce different gene rankings will produce different DE results. The use of absolute RNA expression instead of relative abundance improves sensitivity and reduces false discoveries [<a href="#ref-5">5</a>]. If the normalization method changed between runs, the DE results will differ.

The diagnostic check is to run the analysis with two normalization methods and compare the overlap of DE genes. If the overlap is low, the results are sensitive to normalization.

Failure Pattern 6: Low Sequencing Depth Produces Unstable Results

Low sequencing depth increases data sparsity and reduces the reliability of DE analysis. For low-depth data, single-cell techniques based on zero-inflation models deteriorate performance [<a href="#ref-2">2</a>]. If the sequencing depth differs between runs, the DE results will differ.

The diagnostic check is to compare the median UMI counts per cell between runs. If the depth is substantially different, the comparison may not be valid.

Limitations of Single-Cell Differential Expression Analysis

Statistical Power Depends on the Number of Biological Replicates

The number of biological replicates is the primary determinant of statistical power for DE analysis. Single-cell data provide many cells per sample, but the effective sample size is the number of biological replicates. A study with three biological replicates per condition has limited power to detect small expression differences, regardless of the number of cells sequenced.

The critical assessment of DE for bulk RNA-Seq projects noted that any estimate from a simplification would be an underestimation of the true need for methods testing for every project [<a href="#ref-11">11</a>]. The same principle applies to single-cell data: the complexity of the analysis means that methods testing is needed for every project.

Data Sparsity Limits Detection of Low-Expression Genes

Genes with low expression levels are detected in few cells, which limits the statistical power to detect DE for these genes. The excessive zeros in single-cell data are a major challenge in DE analysis [<a href="#ref-5">5</a>]. Methods that model zero inflation may perform poorly for low-depth data [<a href="#ref-2">2</a>].

Batch Effects Cannot Be Fully Removed

Data integration methods can reduce batch effects but cannot fully remove them. The benchmark of integration strategies found that batch effects, sequencing depth, and data sparsity substantially impact the performance of DE workflows [<a href="#ref-2">2</a>]. Researchers should interpret DE results from integrated data with caution.

Cell Type Composition Affects Results

The proportion of each cell type in a sample affects the DE results. If the cell type composition differs between conditions, DE genes may reflect composition differences instead of expression changes within cell types. The protocol for identifying recirculating thymic regulatory T cells and characterizing the role of Eos in their function using scRNA-seq and TCR-seq describes a refined classification based on analysis of gene expression from single-cell RNA sequencing and TCR sequencing data [<a href="#ref-15">15</a>]. This type of refined classification is necessary to distinguish cell type composition effects from expression changes.

Safety and Regulatory Context for Single-Cell Data Analysis

Data Management and Privacy

Single-cell data from human subjects are subject to privacy regulations. The NCBI provides official descriptions of databases and search systems for accessing public datasets [<a href="#ref-10">10</a>]. Researchers should ensure that their data management practices comply with applicable regulations and institutional policies.

Reproducibility Requirements for Publication

Many journals require that analysis code and parameters be made available for reproducibility. The nf-core documentation describes community pipeline standards for reproducible workflow context [<a href="#ref-7">7</a>]. The Galaxy Training Network provides accessible workflow training and analysis tutorials with reproducibility context [<a href="#ref-6">6</a>]. Researchers should prepare their analysis code and documentation for public release.

Quality Control Standards

Quality control standards for single-cell data are still evolving. The Bioconductor project provides official package, workflow, installation, and reproducible genomic-analysis documentation [<a href="#ref-8">8</a>]. Researchers should follow established quality control practices and document their decisions.

Professional Escalation Criteria

When to Seek Expert Assistance

The following situations warrant consultation with a bioinformatics expert or statistician:

  1. DE results are not reproducible across runs despite following the diagnostic checks described in this article
  2. Batch effects are confounded with the biological condition of interest
  3. The number of biological replicates is very small and the DE results are unstable
  4. The analysis requires advanced statistical methods that are not familiar to the research team
  5. The DE results are critical for a regulatory submission or clinical decision

When to Reconsider the Experimental Design

If the diagnostic checks reveal that the experimental design is inadequate for the research question, the experiment may need to be redesigned. The EMBL-EBI Training provides bioinformatics learning pathways and data-resource training that cover experimental design considerations [<a href="#ref-14">14</a>]. Common design problems include:

  1. Too few biological replicates
  2. Batch effects confounded with condition
  3. Inadequate sequencing depth
  4. Cell type composition differences between conditions that are not accounted for

A Decision Framework for Selecting Differential Expression Methods Based on Data Characteristics

When DE results are not reproducible, the cause is often a mismatch between the statistical method and the structural features of the dataset. instead of testing every available method, researchers can apply a structured decision framework that matches method choice to three measurable data characteristics: biological replicate count, sequencing depth, and batch effect severity. This framework converts the abstract problem of method selection into concrete diagnostic steps that can be completed in a single working session.

Step 1: Measure the Effective Sample Size

The first decision point is the number of biological replicates per condition. This number, not the number of cells, determines the statistical foundation of the analysis. Methods that ignore variation between biological replicates are biased and prone to false discoveries, and the most widely used methods can discover hundreds of differentially expressed genes in the absence of biological differences [<a href="#ref-3">3</a>]. The effective sample size determines which class of methods is appropriate.

For datasets with three or fewer biological replicates per condition, the analysis operates in a low-replicate regime. In this regime, the researcher should prioritize methods that borrow information across genes to stabilize dispersion estimates. A critical assessment of differential expression for bulk RNA-Seq projects found potential value in testing the use of the generalized linear model implementation of edgeR with robust dispersion estimation more frequently for either single-variate or multi-variate two-group comparisons [<a href="#ref-11">11</a>]. This recommendation extends to pseudobulk analysis of single-cell data, where the same dispersion challenges apply.

For datasets with four to ten biological replicates per condition, the analysis operates in a moderate-replicate regime. The same assessment noted that when considering a limited number of patient sample comparisons with larger sample size, there might be some decreased variability between methods [<a href="#ref-11">11</a>]. In this regime, the researcher has more flexibility in method choice, and the decision should shift to the other two data characteristics.

For datasets with more than ten biological replicates per condition, the analysis operates in a high-replicate regime. A method of moments framework for differential expression analysis of scRNA-seq data was introduced for robust and efficient differential analysis of mean expression, variability, and gene correlation from scRNA-seq data scalable to millions of cells and thousands of samples [<a href="#ref-12">12</a>]. This framework identified more significant and reproducible differences in mean expression compared to existing methods [<a href="#ref-12">12</a>]. The scalability of this approach makes it suitable for large cohort studies where the number of samples is the limiting factor.

Step 2: Assess Sequencing Depth and Sparsity

The second decision point is sequencing depth, which directly determines data sparsity. A benchmark of 46 workflows for differential expression analysis of single-cell data with multiple batches showed that batch effects, sequencing depth, and data sparsity substantially impact performance [<a href="#ref-2">2</a>]. The same study found that for low-depth data, single-cell techniques based on zero-inflation models deteriorate performance, whereas the analysis of uncorrected data using limmatrend, Wilcoxon test, and fixed effects models performs well [<a href="#ref-2">2</a>].

To measure sequencing depth, calculate the median UMI counts per cell for each biological sample. If the median is below approximately 1000 UMI counts per cell, the data are in the low-depth regime. In this regime, avoid methods that explicitly model zero inflation. The benchmark evidence indicates that these methods perform poorly when the data are already sparse [<a href="#ref-2">2</a>]. Instead, use methods that are robust to sparsity, such as limmatrend, the Wilcoxon test, or fixed effects models applied to uncorrected data [<a href="#ref-2">2</a>].

If the median is above approximately 5000 UMI counts per cell, the data are in the high-depth regime. In this regime, the researcher has more method options because the sparsity problem is less severe. A statistical framework that leverages UMI counts and zero proportions within a generalized Poisson or binomial mixed-effects model was proposed to account for batch effects and within-sample variation [<a href="#ref-5">5</a>]. This framework is more adaptable to diverse experimental designs and mitigates key shortcomings of current approaches, particularly those related to normalization procedures [<a href="#ref-5">5</a>].

Step 3: Quantify Batch Effect Severity

The third decision point is the severity of batch effects. The benchmark of integration strategies found that the use of batch-corrected data rarely improves the analysis for sparse data, whereas batch covariate modeling improves the analysis for substantial batch effects [<a href="#ref-2">2</a>]. This finding provides a clear decision rule: the choice between data integration and batch covariate modeling depends on the severity of the batch structure.

To quantify batch effect severity, visualize the data with PCA or UMAP colored by batch before any integration step. If cells from different batches form distinct clusters that do not overlap with biological condition, the batch effect is substantial. In this case, include batch as a covariate in the DE model instead of correcting the data before analysis. The benchmark evidence supports batch covariate modeling for substantial batch effects [<a href="#ref-2">2</a>].

If cells from different batches intermingle and the clustering is driven primarily by biological condition, the batch effect is mild. In this case, the researcher can proceed with uncorrected data or use batch correction without substantial risk of altering the DE signal. The benchmark found that the use of batch-corrected data rarely improves the analysis for sparse data [<a href="#ref-2">2</a>], so for mild batch effects in sparse data, the simpler approach of uncorrected data with batch as a covariate is preferred.

Step 4: Apply the Decision Matrix

The three measured characteristics combine into a decision matrix that guides method selection. The matrix below summarizes the recommended approach for each combination of replicate count, sequencing depth, and batch effect severity.

Biological ReplicatesSequencing DepthBatch Effect SeverityRecommended Approach
Low (3 or fewer)LowMildPseudobulk with edgeR and robust dispersion estimation
Low (3 or fewer)LowSubstantialPseudobulk with batch covariate in the model
Low (3 or fewer)HighMildPseudobulk with edgeR or DESeq2
Low (3 or fewer)HighSubstantialPseudobulk with batch covariate and robust dispersion
Moderate (4 to 10)LowMildCell-level analysis with limmatrend or Wilcoxon test
Moderate (4 to 10)LowSubstantialCell-level analysis with batch covariate modeling
Moderate (4 to 10)HighMildCell-level analysis with mixed-effects model
Moderate (4 to 10)HighSubstantialCell-level analysis with batch covariate and mixed-effects model
High (more than 10)LowMildMethod of moments framework or limmatrend
High (more than 10)LowSubstantialMethod of moments framework with batch covariate
High (more than 10)HighMildMethod of moments framework or mixed-effects model
High (more than 10)HighSubstantialMethod of moments framework with batch covariate

The decision matrix is a starting point, not a final answer. The critical assessment of DE for bulk RNA-Seq projects noted that any estimate from a simplification would be an underestimation of the true need for methods testing for every project [<a href="#ref-11">11</a>]. The same principle applies to single-cell data: the complexity of the analysis means that methods testing is needed for every project [<a href="#ref-11">11</a>]. The matrix narrows the method space but does not eliminate the need for empirical validation.

Step 5: Validate the Selected Method with a Known Signal

After selecting a method from the decision matrix, validate the choice using a known biological signal. The critical assessment of DE used knock-out and over-expression studies as a simplification to test recovery of a known causal gene in RNA-Seq cell line experiments [<a href="#ref-11">11</a>]. The same validation strategy applies to single-cell data.

Identify a gene or set of genes with a well-established expression difference between the conditions being compared. This could be a marker gene for a cell type that is known to change in abundance or a gene with a documented response to the perturbation under study. Run the selected DE method and check whether the known signal is recovered. If the known gene is not detected as differentially expressed, the method may be too conservative for the data characteristics, and an alternative from the decision matrix should be tested.

A reproducible Seurat-based protocol for analyzing PBMC CD4+ T-cell single-cell RNA sequencing data across malaria reinfection timepoints demonstrates the value of a unified computational framework that includes standardized preprocessing, integration, clustering, and downstream transcriptomic analyses [<a href="#ref-4">4</a>]. The protocol enables systematic computation of module scores for predefined immune and CD4+ T-cell programs, identification of cluster-specific marker genes, timepoint-resolved differential expression analysis, and downstream Gene Ontology and KEGG pathway enrichment [<a href="#ref-4">4</a>]. The use of predefined gene programs provides a built-in validation set for the DE analysis.

Step 6: Document the Decision and Record the Outcome

The final step in the decision framework is documentation. Record the measured data characteristics, the selected method, the validation result, and the final DE gene list. The nf-core documentation describes community pipeline standards for reproducible workflow context [<a href="#ref-7">7</a>]. These standards include version tracking and parameter documentation. The Galaxy Training Network provides accessible workflow training and analysis tutorials with reproducibility context [<a href="#ref-6">6</a>].

The documentation should include the median UMI counts per cell, the number of biological replicates per condition, the batch effect severity assessment, the selected method and its version, the validation gene set and whether it was recovered, and the final DE gene list with the significance threshold. This record enables the researcher to reproduce the analysis at a later date and provides a basis for comparison when new data are added.

Common Failure Patterns in Method Selection

The decision framework addresses several common failure patterns that produce non-reproducible DE results. The first failure pattern is using a cell-level method that ignores biological replicate variation when the replicate count is low. This pattern produces hundreds of false discoveries [<a href="#ref-3">3</a>]. The framework prevents this by directing low-replicate datasets to pseudobulk methods.

The second failure pattern is using a zero-inflation model on low-depth data. The benchmark found that single-cell techniques based on zero-inflation models deteriorate performance for low-depth data [<a href="#ref-2">2</a>]. The framework prevents this by directing low-depth datasets to limmatrend, Wilcoxon test, or fixed effects models.

The third failure pattern is using batch-corrected data for DE analysis when the data are sparse. The benchmark found that the use of batch-corrected data rarely improves the analysis for sparse data [<a href="#ref-2">2</a>]. The framework prevents this by directing sparse datasets to uncorrected data with batch covariate modeling.

The fourth failure pattern is using a method that is not scalable to the dataset size. A method of moments framework was introduced for robust and efficient differential analysis scalable to millions of cells and thousands of samples [<a href="#ref-12">12</a>]. For large datasets, methods that cannot scale will produce unstable results or fail to complete. The framework directs high-replicate datasets to scalable methods.

Limitations of the Decision Framework

The decision framework has limitations that should be acknowledged. The thresholds for low, moderate, and high replicate counts are practical guidelines, not statistical guarantees. The boundary between low and moderate depth at approximately 1000 UMI counts per cell is similarly approximate. The framework does not replace the need for methods testing on each dataset.

The critical assessment of DE for bulk RNA-Seq projects noted that analysis of public data does not consider all experimental designs, and presentation of downstream analysis is limited [<a href="#ref-11">11</a>]. Any estimate from a simplification would be an underestimation of the true need for methods testing for every project [<a href="#ref-11">11</a>]. The decision framework is a simplification that narrows the method space but does not eliminate the need for empirical validation.

The framework also does not address all sources of non-reproducibility. Cell type composition differences between conditions can produce DE genes that reflect composition changes instead of expression changes within cell types. A protocol for identifying recirculating thymic regulatory T cells and characterizing the role of Eos in their function using scRNA-seq and TCR-seq describes a refined classification based on analysis of gene expression from single-cell RNA sequencing and TCR sequencing data [<a href="#ref-15">15</a>]. This type of refined classification is necessary to distinguish cell type composition effects from expression changes [<a href="#ref-15">15</a>].

Professional Escalation Criteria for Method Selection

The decision framework provides a structured approach, but some situations warrant consultation with a bioinformatics expert or statistician. If the DE results remain non-reproducible after applying the framework and validating with a known signal, expert assistance is needed. If the dataset has a complex experimental design with multiple covariates or nested batch structures, the framework may not capture the full complexity. If the DE results are critical for a regulatory submission or clinical decision, expert review is required.

The EMBL-EBI Training provides bioinformatics learning pathways and data-resource training that cover experimental design considerations [<a href="#ref-14">14</a>]. The Bioconductor project provides official package, workflow, installation, and reproducible genomic-analysis documentation [<a href="#ref-8">8</a>]. These resources can support researchers who need to deepen their understanding of method selection before consulting an expert.

Frequently Asked Questions

Why do I get different DE genes when I run the same analysis twice?

If the same analysis is run twice with the same software versions and parameters, the results should be identical unless the pipeline includes stochastic steps. Check whether any step uses random initialization or subsampling. If the pipeline is deterministic and the results differ, the input data or parameters changed between runs. Verify that the input files are identical and that the parameters are the same.

What is the difference between pseudobulk and cell-level DE analysis?

Pseudobulk analysis aggregates cells from the same biological sample into a single expression profile and uses bulk RNA-seq DE methods. Cell-level analysis treats each cell as an observation and uses methods designed for single-cell data. Pseudobulk analysis preserves the biological replicate structure and produces more conservative results. Cell-level analysis has higher power but can produce false discoveries when biological replicate variation is ignored [<a href="#ref-3">3</a>].

Should I use batch-corrected data for DE analysis?

The use of batch-corrected data rarely improves the analysis for sparse data, whereas batch covariate modeling improves the analysis for substantial batch effects [<a href="#ref-2">2</a>]. Test both approaches and compare the results. If the DE gene lists are substantially different, the batch effect is influencing the results and the choice should be documented.

How many biological replicates do I need for reproducible DE results?

The number of biological replicates depends on the effect size and the variability between samples. More replicates produce more reproducible results. A study with three biological replicates per condition has limited power to detect small expression differences. The critical assessment of DE for bulk RNA-Seq projects noted that finding the right balance of quality and quantity can be important [<a href="#ref-11">11</a>].

Why do I get hundreds of DE genes when there should be no biological difference?

This result indicates that the statistical method is producing false discoveries. Methods that ignore variation between biological replicates are biased and prone to false discoveries, and the most widely used methods can discover hundreds of differentially expressed genes in the absence of biological differences [<a href="#ref-3">3</a>]. Switch to a method that accounts for biological replicate variation.

What should I do if my DE results change when I change the filtering thresholds?

Run the filtering step with a range of thresholds and compare the number of cells retained and the number of DE genes identified. If the results are sensitive to threshold choices, the thresholds should be specified in advance and documented. Consider using a pseudobulk approach, which is less sensitive to cell-level filtering choices.

How do I know if my normalization method is affecting my DE results?

Run the analysis with two different normalization methods and compare the overlap of DE genes. If the overlap is low, the results are sensitive to normalization. A statistical framework that uses absolute RNA expression instead of relative abundance improves sensitivity and reduces false discoveries [<a href="#ref-5">5</a>].

What is the best DE method for single-cell data?

There is no single best method for all situations. The choice of method depends on the experimental design, the depth of the sequencing data, and the severity of batch effects. A benchmark of 46 workflows suggested several high-performance methods under different conditions based on simulation and real data analyses [<a href="#ref-2">2</a>]. Test multiple methods and compare the results to identify the most robust approach for your data.

Related Bioinformatics Guides

Related Clinical & Scientific Guides

References and Further Reading

[1] [Bayesian approach to single-cell differential expression analysis](https://doi.org/10.1038/nmeth.2967). Nature Methods, 2014. [2] [Benchmarking integration of single-cell differential expression](https://doi.org/10.1038/s41467-023-37126-3). Nature Communications, 2023. [3] [Confronting false discoveries in single-cell differential expression](https://doi.org/10.1038/s41467-021-25960-2). Nature Communications, 2021. [4] [A Reproducible Seurat-Based Protocol for Single-Cell RNA Sequencing Analysis of Peripheral Blood Mononuclear Cell CD4+ T Cells During Malaria Reinfection.](https://pubmed.ncbi.nlm.nih.gov/42612107). Journal of visualized experiments : JoVE, 2026. [5] [Exploring and mitigating shortcomings in single-cell differential expression analysis with a new statistical paradigm](https://doi.org/10.1186/s13059-025-03525-6). Genome Biology, 2025. [6] [Galaxy Training Network](https://training.galaxyproject.org/). Galaxy Project. [7] [nf-core Documentation](https://nf-co.re/docs). nf-core. [8] [Bioconductor](https://bioconductor.org/). Bioconductor Project. [9] [Isolation of Viable Adipocytes and Stromal Vascular Fraction from Human Visceral Adipose Tissue Suitable for RNA Analysis and Macrophage Phenotyping.](https://pubmed.ncbi.nlm.nih.gov/33191941). Journal of visualized experiments : JoVE, 2020. [10] [NCBI Data Resources](https://www.ncbi.nlm.nih.gov/). National Center for Biotechnology Information. [11] [Critical Differential Expression Assessment for Individual Bulk RNA-Seq Projects.](https://pubmed.ncbi.nlm.nih.gov/38405814). bioRxiv : the preprint server for biology, 2024. [12] [Method of moments framework for differential expression analysis of single-cell RNA-sequencing data](https://doi.org/10.1016/j.cell.2024.09.044). Cell, 2024. [13] [The Carpentries Lessons](https://carpentries.org/lessons). The Carpentries. [14] [EMBL-EBI Training](https://www.ebi.ac.uk/training). European Bioinformatics Institute. [15] [Protocol for identifying recirculating thymic regulatory T cells and characterizing the role of Eos in their function using scRNA-seq and TCR-seq.](https://doi.org/10.1016/j.xpro.2026.104620). 2026.

This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.