# How to Perform Differential Expression Analysis in Single-Cell RNA-Seq: A Step-by-Step Workflow with Seurat and DESeq2


## Key Takeaways

- Pseudobulk analysis aggregates raw counts from individual cells within a biological sample and cell type, creating a matrix suitable for bulk RNA-seq differential expression (DE) tools like DESeq2. This approach preserves biological replication and accounts for inter-sample variability, mitigating the inflated significance and false positives inherent in applying bulk methods directly to single-cell data.
- The pseudobulk workflow necessitates robust upstream single-cell processing, including stringent quality control to remove low-quality cells and doublets, accurate cell type annotation based on marker genes or reference mapping, and complete sample metadata detailing biological replicates and covariates (e.g., batch, sex).
- DESeq2 is employed on the pseudobulk count matrix, utilizing a design formula that incorporates relevant covariates (e.g., `~ batch + condition`) to adjust for confounding factors and accurately estimate gene expression differences between experimental conditions. Benjamini-Hochberg correction is applied to control the false discovery rate for multiple testing.
- Insufficient biological replication (minimum of 3 samples per condition recommended) is a critical failure point, leading to unstable dispersion estimates and reduced statistical power in DESeq2. Pseudoreplication, arising from treating individual cells as independent replicates, is another common pitfall that inflates significance and yields unreliable DE calls.
- While pseudobulk analysis is effective for detecting DE, it inherently loses single-cell resolution and is dependent on the accuracy of cell type annotations. It does not detect compositional changes in cell type abundance between conditions, necessitating complementary differential abundance testing methods if such changes are of interest.

---

Differential expression (DE) analysis in single-cell RNA sequencing (scRNA-seq) asks whether gene expression levels differ between defined groups of cells, such as cell types, treatment conditions, or disease states. The core problem researchers face is that applying bulk RNA-seq DE methods directly to single-cell count matrices produces inflated significance and false positives because individual cells within a sample are not independent replicates. The practical solution is pseudobulk analysis, where counts are aggregated across all cells of a given cell type within each biological sample, followed by DE testing with tools designed for bulk data such as DESeq2. This workflow preserves biological replication, accounts for inter-sample variability, and produces results that are more reproducible across studies. This article provides a complete code walkthrough for performing pseudobulk DE analysis using Seurat for single-cell processing and DESeq2 for statistical testing, including covariate handling, multiple testing correction, and interpretation of results.

## Understanding Why Single-Cell DE Analysis Requires a Distinct Workflow

Single-cell RNA sequencing measures gene expression in thousands of individual cells from a single experiment, enabling researchers to resolve cellular heterogeneity that is invisible in bulk tissue measurements. The technology has transformed genomics by allowing cell-type-specific transcriptome profiling at unprecedented resolution. However, the analytical complexity of scRNA-seq data demands computational tools that differ substantially from those used for population-based RNA sequencing. A universal standardization of analytical methods does not yet exist, which can make it difficult for newcomers to navigate the field. Despite this lack of standardization, several canonical workflow steps have emerged as common practice, including read mapping, quality control, gene expression quantification, normalization, feature selection, dimensionality reduction, and cell clustering, with differential expression analysis serving as a downstream goal for many studies.

The statistical challenge in scRNA-seq DE analysis stems from the data structure itself. In a typical experiment, you may have thousands of cells sequenced from each of several biological samples, such as individual patients or animals. If you treat each cell as an independent observation, you are effectively claiming that all cells within a sample are independent replicates. This assumption is violated because cells from the same sample share technical batch effects, environmental influences, and genetic background. The result is that p-values become artificially small and the number of genes called as differentially expressed becomes inflated. Pseudobulk aggregation solves this problem by summing counts across all cells of a given cell type within each biological sample, creating one expression value per gene per sample per cell type. These aggregated values are then analyzed with DE methods designed for bulk RNA-seq data, which model biological variability between samples instead of technical variability between cells.

## At a Glance: Key Decisions in the Pseudobulk DE Workflow

The table below summarizes the major decision points in a pseudobulk DE analysis workflow, the recommended approach, and the rationale for each choice.

| Workflow Step | Recommended Approach | Rationale |
| --- | --- | --- |
| Data input | Processed Seurat object with cell type annotations and sample metadata | Enables aggregation by cell type and biological sample |
| Aggregation strategy | Sum raw counts across cells within each sample-cell type combination | Preserves count nature of data and biological replication structure |
| Minimum sample size | At least 3 biological samples per condition | Provides sufficient replication for variance estimation in DESeq2 |
| DE testing tool | DESeq2 on pseudobulk count matrix | Models negative binomial distribution and handles small sample sizes |
| Covariate handling | Include batch, sex, age, or other known confounders in the design formula | Reduces false positives from technical and biological variation |
| Multiple testing correction | Benjamini-Hochberg false discovery rate control | Balances discovery power with false positive control |
| Filtering thresholds | Adjusted p-value below 0.05 and log2 fold change above a defined cutoff | Focuses interpretation on biologically meaningful changes |

## Preparing Your Seurat Object for Pseudobulk Aggregation

Before you can perform pseudobulk DE analysis, your Seurat object must contain the necessary metadata and be properly processed. The quality of your downstream DE results depends directly on the quality of the upstream processing steps, including quality control, normalization, and cell type annotation.

### Quality Control and Cell Type Annotation

Quality control is the first critical step in any scRNA-seq analysis. The goal is to remove low-quality cells, doublets, and empty droplets that would otherwise introduce noise into downstream analyses. Standard quality control metrics include the number of genes detected per cell, the total number of unique molecular identifiers (UMIs), and the percentage of mitochondrial reads. Cells with very low gene counts may represent empty droplets or dying cells, while cells with very high gene counts may represent doublets where two cells were captured together. High mitochondrial read fractions typically indicate stressed or dying cells. The specific thresholds for these metrics depend on your tissue type and experimental protocol, and you should examine the distributions in your own data instead of applying arbitrary cutoffs.

After quality control, the next steps are normalization, feature selection, dimensionality reduction, and clustering. These steps allow you to identify distinct cell populations in your data. Cell type annotation can be performed manually based on known marker genes, or through reference-based label transfer approaches that leverage previously annotated datasets. The accuracy of your cell type annotations directly affects the validity of your DE results, since pseudobulk aggregation assumes that all cells within a cluster represent the same cell type. Misannotated clusters will produce pseudobulk profiles that mix distinct cell types, diluting true DE signals and potentially creating spurious ones.

### Ensuring Sample Metadata Is Complete

For pseudobulk analysis, your Seurat object must contain metadata identifying which biological sample each cell came from. This is often stored in a column such as "sample_id" or "orig.ident". You also need metadata for any covariates you plan to include in your DE model, such as treatment condition, sex, age, or batch. If this information is not already present in your Seurat object, you need to add it before proceeding. The sample-level metadata should be stored in a separate data frame that maps each sample ID to its condition and covariate values. This data frame will be used when constructing the DESeq2 model.

The importance of complete sample metadata cannot be overstated. If you forget to record which cells came from which biological sample, you cannot perform pseudobulk analysis at all. If you omit important covariates such as batch or sex, your DE results may be confounded. Before starting any scRNA-seq experiment, you should plan your metadata collection carefully, ensuring that every sample has a unique identifier and that all relevant experimental variables are recorded.

## Creating Pseudobulk Count Matrices from Seurat

The pseudobulk approach aggregates raw counts across all cells of a given cell type within each biological sample. This creates a matrix where rows are genes, columns are sample-cell type combinations, and values are the summed counts across all cells in that group. The aggregation must be performed on raw counts, not normalized or scaled data, because DESeq2 expects count data and performs its own normalization internally.

### The Aggregation Function

The following R code demonstrates how to create a pseudobulk count matrix from a Seurat object. This function uses the `AggregateExpression` function from Seurat, which sums counts across cells based on grouping variables.

```r
## Load required libraries
library(Seurat)
library(dplyr)
library(tidyr)

## Function to create pseudobulk matrix
CreatePseudobulk <- function(seurat_obj, sample_col, celltype_col) {
  # Get raw counts
  counts <- GetAssayData(seurat_obj, assay = "RNA", slot = "counts")

  # Create a combined identifier for sample and cell type
  seurat_obj@meta.data$pseudobulk_id <- paste(
    seurat_obj@meta.data[[sample_col]],
    seurat_obj@meta.data[[celltype_col]],
    sep = "_"
  )

  # Aggregate counts by pseudobulk ID
  pseudobulk <- AggregateExpression(
    seurat_obj,
    assays = "RNA",
    group.by = "pseudobulk_id",
    return.seurat = FALSE
  )

  # Extract the count matrix
  pseudobulk_matrix <- as.matrix(pseudobulk$RNA)

  return(pseudobulk_matrix)
}

## Example usage
## pseudobulk_matrix <- CreatePseudobulk(seurat_obj, "sample_id", "cell_type")
```

The resulting matrix has genes as rows and sample-cell type combinations as columns. Each column name encodes both the sample and cell type, which you will need to parse when constructing the DESeq2 dataset.

### Handling Single-Cell and Single-Nucleus Data

The pseudobulk workflow applies equally to single-cell RNA-seq and single-nucleus RNA-seq (snRNA-seq) data. Single-nucleus sequencing is particularly useful for tissues that are difficult to dissociate into intact cells, such as brain tissue or archived frozen samples. The computational workflow for snRNA-seq data follows the same steps as scRNA-seq, including quality control, normalization, clustering, and pseudobulk aggregation. However, the quality control metrics may differ slightly, since nuclear transcripts have different characteristics than whole-cell transcripts. For example, the percentage of mitochondrial reads is typically lower in single-nucleus data because mitochondria are largely excluded from nuclear preparations.

The choice between single-cell and single-nucleus approaches depends on your biological question and tissue type. Fresh tissues that dissociate well are suitable for single-cell sequencing, while frozen tissues or tissues with complex morphology may require single-nucleus approaches. The DE analysis workflow described here works for both data types, provided you have performed appropriate quality control and cell type annotation for your specific data modality.

## Running DESeq2 on Pseudobulk Count Data

Once you have created your pseudobulk count matrix, the next step is to run DESeq2 for differential expression testing. DESeq2 is a widely used Bioconductor package for analyzing count-based transcriptomic data, and it is well suited for pseudobulk scRNA-seq analysis because it models the negative binomial distribution and handles small sample sizes appropriately.

### Constructing the DESeq2 Dataset

To run DESeq2, you need to create a DESeqDataSet object that contains the count matrix, sample metadata, and a design formula. The design formula specifies which variables explain the observed expression differences. At minimum, the design formula should include the condition of interest. If you have known covariates such as batch, sex, or age, these should be included in the design formula to account for their effects.

```r
## Load DESeq2
library(DESeq2)

## Parse the column names of the pseudobulk matrix to extract sample and cell type
## Assuming column names are in the format "sample_celltype"
col_data <- data.frame(
  pseudobulk_id = colnames(pseudobulk_matrix),
  stringsAsFactors = FALSE
)

## Split the pseudobulk ID into sample and cell type
col_data <- col_data %>%
  separate(pseudobulk_id, into = c("sample", "cell_type"), sep = "_", remove = FALSE)

## Merge with sample metadata
## sample_metadata should be a data frame with columns for sample ID, condition, and covariates
col_data <- col_data %>%
  left_join(sample_metadata, by = "sample")

## Set row names
rownames(col_data) <- col_data$pseudobulk_id

## Create DESeq2 dataset
dds <- DESeqDataSetFromMatrix(
  countData = pseudobulk_matrix,
  colData = col_data,
  design = ~ batch + condition
)
```

The design formula `~ batch + condition` indicates that you want to test for differences due to condition while accounting for batch effects. The order of terms matters: the variable of interest should come last in the formula. If you have multiple covariates, list them before the condition variable.

### Running the DE Analysis

After creating the DESeqDataSet, you run the DESeq function to perform the analysis. This function performs several steps internally, including estimation of size factors, estimation of dispersion, and fitting of the negative binomial model.

```r
## Run DESeq2 analysis
dds <- DESeq(dds)

## Get results for the condition comparison
## Specify the contrast you are interested in
results <- results(dds, contrast = c("condition", "treated", "control"))

## View summary of results
summary(results)

## Convert to data frame for further analysis
results_df <- as.data.frame(results)
results_df$gene <- rownames(results_df)
```

The results object contains columns for baseMean, log2FoldChange, lfcSE, stat, pvalue, and padj. The padj column contains the Benjamini-Hochberg adjusted p-values, which control the false discovery rate across all tested genes.

### Filtering Significant Genes

After obtaining the results, you typically filter for genes that meet your significance thresholds. A common approach is to use an adjusted p-value cutoff of 0.05 and a log2 fold change cutoff of 1 (corresponding to a 2-fold change). However, the appropriate thresholds depend on your biological context and the power of your experiment.

```r
## Filter for significant genes
significant_genes <- results_df %>%
  filter(padj < 0.05, abs(log2FoldChange) > 1)

## Order by adjusted p-value
significant_genes <- significant_genes %>%
  arrange(padj)

## View top genes
head(significant_genes)
```

The number of significant genes you obtain depends on several factors, including the biological difference between conditions, the number of biological replicates, the depth of sequencing, and the variability between samples. Experiments with more biological replicates and larger effect sizes will yield more significant genes.

## Handling Covariates and Batch Effects

Covariates are variables that are not your primary interest but may influence gene expression and confound your DE analysis if not accounted for. Common covariates in scRNA-seq studies include sequencing batch, sex, age, and genetic background. Ignoring these variables can lead to false positives or false negatives in your DE results.

### Including Covariates in the Design Formula

The design formula in DESeq2 allows you to include multiple covariates. For example, if your samples come from two sequencing batches and include both male and female individuals, your design formula might be `~ batch + sex + condition`. This model estimates the effect of condition while adjusting for batch and sex differences.

The choice of which covariates to include requires careful thought. Including too many covariates can reduce your statistical power, especially with small sample sizes. Including too few can leave confounding variables unaccounted for. A reasonable approach is to include covariates that you know or suspect affect gene expression based on your experimental design or prior knowledge. You should also check for associations between your covariates and your condition of interest. If a covariate is perfectly confounded with your condition, you cannot separate their effects, and your DE analysis will be unreliable.

### Technical Batch Effects in Single-Cell Data

Technical batch effects are a major concern in scRNA-seq experiments. Differences in sample preparation, sequencing runs, and reagent lots can introduce systematic variation that is unrelated to biology. The multiplexed droplet single-cell RNA-sequencing approach using natural genetic variation, implemented in the tool demuxlet, allows cells from multiple individuals to be pooled and captured in a single experiment, reducing batch effects and increasing throughput. This approach uses natural genetic variation to determine the sample identity of each droplet and to detect doublets. In pools of up to 64 individuals, as few as 50 single-nucleotide polymorphisms per cell are sufficient to assign the correct sample identity for the vast majority of singlets and to identify doublets at rates consistent with previous estimates.

If you have used multiplexing approaches such as demuxlet, your sample metadata should include the genetic identity of each cell, which you can then use for pseudobulk aggregation. The batch effects that would have arisen from processing samples separately are reduced because all samples were processed together. However, you should still check for any remaining batch structure in your data and include batch as a covariate if necessary.

### Integration Approaches and Their Limitations

Data integration methods such as Harmony, Seurat's integration functions, or scVI can be used to align single-cell datasets across batches or conditions. These methods are valuable for visualization and clustering, as they remove technical variation and allow cells from different batches to be compared directly. However, for DE analysis, you should be cautious about using integrated data. Integration methods alter the expression values, and the corrected values are not suitable for count-based DE analysis with DESeq2. The recommended approach is to perform clustering and cell type annotation on integrated data, then return to the raw counts for pseudobulk aggregation and DE testing. This separation of concerns ensures that your DE analysis uses unmodified count data while still benefiting from the improved cell type identification that integration provides.

## Multiple Testing Correction and Interpretation

When you test thousands of genes for differential expression, you face a multiple testing problem. If you use a nominal p-value threshold of 0.05, you would expect 5% of all tested genes to appear significant by chance alone. With 20,000 genes tested, this means approximately 1,000 false positives. Multiple testing correction methods address this issue by adjusting p-values to control the rate of false discoveries.

### Benjamini-Hochberg False Discovery Rate Control

The Benjamini-Hochberg procedure is the most commonly used multiple testing correction method in genomics. It controls the false discovery rate, which is the expected proportion of false positives among all genes called significant. DESeq2 applies this correction by default and reports the adjusted p-values in the padj column. An adjusted p-value of 0.05 means that approximately 5% of the genes called significant are expected to be false positives.

The false discovery rate approach is preferred over more stringent methods such as Bonferroni correction because it provides a better balance between discovering true positives and controlling false positives. Bonferroni correction controls the family-wise error rate, which is the probability of making even one false positive, but it is very conservative and may miss many true DE genes, especially in experiments with limited power.

### Interpreting Fold Changes and Effect Sizes

While adjusted p-values tell you whether a gene's expression change is statistically significant, the log2 fold change tells you the magnitude of the change. A gene with a log2 fold change of 1 has doubled its expression, while a gene with a log2 fold change of -1 has halved its expression. The biological significance of a fold change depends on the gene and the context. Some genes have large expression changes that are biologically meaningful, while others have small changes that may still be important if the gene is highly regulated.

DESeq2 applies shrinkage to log2 fold changes, which moderates the fold change estimates for genes with low counts or high variability. This shrinkage reduces the tendency for genes with low expression to show artificially large fold changes. The shrunken fold changes are more reliable for ranking genes and for downstream analyses such as gene set enrichment.

### Complementing DE with Differential Detection Analysis

Standard DE analysis tests for differences in the average expression of genes between conditions. However, single-cell data also allows you to test for differences in the fraction of cells or samples in which a gene is detected, an approach called differential detection (DD) analysis. DE and DD analyses provide complementary information, both in terms of the individual genes they report and in the functional interpretation of those genes. A gene may show similar average expression between conditions but be detected in a much larger fraction of cells in one condition, indicating a change in the proportion of cells expressing that gene instead of a change in expression level per cell. Incorporating DD analysis into your workflow can provide a more complete picture of transcriptional changes between conditions.

## Practical Workflow: From Seurat to DESeq2 Results

This section provides a complete, step-by-step workflow that you can adapt to your own data. The workflow assumes you have a Seurat object with cell type annotations and sample metadata.

### Step 1: Verify Your Seurat Object

Before creating pseudobulk data, verify that your Seurat object contains the necessary information. Check that the metadata contains sample identifiers and cell type annotations, and that the raw counts are accessible.

```r
## Check metadata columns
head(seurat_obj@meta.data)

## Verify that raw counts are available
GetAssayData(seurat_obj, assay = "RNA", slot = "counts")[1:5, 1:5]

## Check the number of cells per sample and cell type
table(seurat_obj@meta.data$sample_id, seurat_obj@meta.data$cell_type)
```

This verification step helps you catch problems early. If sample identifiers are missing or cell type annotations are incomplete, you need to address these issues before proceeding.

### Step 2: Create the Pseudobulk Matrix

Use the aggregation function described earlier to create the pseudobulk count matrix. After creating the matrix, verify that the dimensions are sensible. The number of columns should equal the number of unique sample-cell type combinations, and the number of rows should equal the number of genes in your dataset.

```r
## Create pseudobulk matrix
pseudobulk_matrix <- CreatePseudobulk(seurat_obj, "sample_id", "cell_type")

## Check dimensions
dim(pseudobulk_matrix)

## Check column names
colnames(pseudobulk_matrix)
```

### Step 3: Filter for Cell Types of Interest

You may not want to test DE for every cell type in your dataset. Focus on cell types that are relevant to your biological question and that have sufficient cell numbers and sample representation. Cell types with very few cells or represented in only a few samples will produce unreliable pseudobulk estimates.

```r
## Define cell types of interest
celltypes_of_interest <- c("T_cells", "B_cells", "Macrophages")

## Filter the pseudobulk matrix
pseudobulk_filtered <- pseudobulk_matrix[, grepl(paste(celltypes_of_interest, collapse = "|"), colnames(pseudobulk_matrix))]
```

### Step 4: Construct the DESeq2 Dataset

Create the column data for the DESeq2 dataset, ensuring that sample metadata is correctly matched to each pseudobulk column.

```r
## Create column data
col_data <- data.frame(
  pseudobulk_id = colnames(pseudobulk_filtered),
  stringsAsFactors = FALSE
)

## Extract sample and cell type
col_data <- col_data %>%
  separate(pseudobulk_id, into = c("sample", "cell_type"), sep = "_", remove = FALSE)

## Merge with sample metadata
col_data <- col_data %>%
  left_join(sample_metadata, by = "sample")

## Ensure row names match column names
rownames(col_data) <- col_data$pseudobulk_id

## Create DESeq2 dataset
dds <- DESeqDataSetFromMatrix(
  countData = pseudobulk_filtered,
  colData = col_data,
  design = ~ batch + condition
)
```

### Step 5: Run DESeq2 and Extract Results

Run the DESeq2 analysis and extract results for your comparison of interest.

```r
## Run DESeq2
dds <- DESeq(dds)

## Extract results
results <- results(dds, contrast = c("condition", "treated", "control"))

## Convert to data frame
results_df <- as.data.frame(results)
results_df$gene <- rownames(results_df)

## Add cell type information
results_df$cell_type <- col_data$cell_type[match(results_df$gene, rownames(results_df))]

## Filter significant genes
sig_genes <- results_df %>%
  filter(padj < 0.05, abs(log2FoldChange) > 1)
```

### Step 6: Visualize and Interpret Results

Visualize your results to check for consistency and to identify patterns. Common visualizations include MA plots, volcano plots, and heatmaps of top DE genes.

```r
## MA plot
plotMA(results, main = "MA Plot", ylim = c(-5, 5))

## Volcano plot
library(ggplot2)
volcano_plot <- ggplot(results_df, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = padj < 0.05 & abs(log2FoldChange) > 1), size = 0.5) +
  scale_color_manual(values = c("grey", "red")) +
  theme_minimal() +
  labs(title = "Volcano Plot", x = "Log2 Fold Change", y = "-Log10 Adjusted P-value")
print(volcano_plot)
```

## Common Failure Patterns and How to Avoid Them

Several recurring problems can compromise pseudobulk DE analysis. Recognizing these failure patterns early can save time and prevent incorrect biological conclusions.

### Insufficient Biological Replication

The most common failure in scRNA-seq DE analysis is having too few biological samples per condition. If you have only two samples per condition, DESeq2 cannot reliably estimate biological variability, and the results will be unstable. With one sample per condition, DE analysis is essentially impossible because there is no replication at the biological level. The minimum recommended number of biological replicates is three per condition, and more replicates provide greater statistical power and more reliable dispersion estimates.

This limitation is particularly relevant for studies using multiplexed approaches. While multiplexing can increase the number of biological samples processed in a single experiment, the number of samples per condition is still determined by your experimental design. Plan your experiment with sufficient biological replication from the start, because you cannot add replicates after sequencing.

### Pseudoreplication from Treating Cells as Independent Samples

The pseudoreplication problem occurs when researchers perform DE analysis directly on single-cell data without aggregation, treating each cell as an independent biological replicate. This approach produces severely inflated significance because the effective sample size is the number of cells, which can be thousands, instead of the number of biological samples, which may be only three or four. The result is that nearly every gene appears differentially expressed, and the top hits are often driven by technical artifacts instead of biology.

The pseudobulk approach avoids this problem by aggregating cells within each biological sample before DE testing. The effective sample size becomes the number of biological samples, which is the correct unit of replication. If you see results where thousands of genes are called significant with very small effect sizes, you may be falling victim to pseudoreplication.

### Incomplete or Incorrect Sample Metadata

Another common failure is incomplete or incorrect sample metadata. If sample identifiers are missing, duplicated, or incorrectly assigned, the pseudobulk aggregation will produce incorrect results. Similarly, if covariate information such as batch or sex is missing, you cannot include these variables in your DE model, and your results may be confounded.

To avoid this problem, maintain a sample metadata file that is separate from your Seurat object and that records all relevant experimental variables. Verify this metadata before performing pseudobulk aggregation, and cross-check that the number of cells per sample matches your expectations.

### Cell Type Misannotation

If your cell type annotations are incorrect, your pseudobulk profiles will mix distinct cell populations, diluting true DE signals and potentially creating spurious ones. Cell type misannotation can occur when clusters are poorly separated, when marker genes are not specific, or when reference-based annotation tools are applied to data from a different tissue or species than the reference.

To reduce the risk of misannotation, validate your cell type labels using multiple marker genes, examine the expression of known cell type markers in your clusters, and consider using reference-based label transfer as a complement to manual annotation. If you are unsure about the identity of a cluster, it is better to exclude it from DE analysis than to include it with an incorrect label.

### Ignoring Ambient RNA Contamination

Ambient RNA contamination occurs when free-floating RNA from lysed cells is captured in droplets along with intact cells, adding spurious transcripts to the measured expression profiles. This contamination can affect DE analysis by introducing background expression that is not specific to the cell type being measured. Tools such as FastCAR have been developed to correct for ambient RNA in scRNA-seq datasets, and using such corrections can improve the accuracy of downstream DE analysis. If you suspect ambient contamination in your data, consider applying an appropriate correction method before pseudobulk aggregation.

## Limitations of Pseudobulk DE Analysis

While pseudobulk analysis is the recommended approach for DE testing in scRNA-seq data, it has limitations that you should understand when interpreting your results.

### Loss of Single-Cell Resolution

Pseudobulk aggregation averages expression across all cells of a given cell type within a sample. This averaging loses information about cell-to-cell variability within the cell type. If your biological question concerns heterogeneity within a cell type, such as the presence of distinct cell states or subpopulations, pseudobulk analysis will not capture this heterogeneity. In such cases, you may need to use alternative approaches that model cell-level variability while properly accounting for biological replication.

### Dependence on Cell Type Annotations

The validity of pseudobulk results depends entirely on the accuracy of your cell type annotations. If your annotations are incorrect or too coarse, your pseudobulk profiles will not represent homogeneous cell populations. For example, if you annotate all T cells as a single group when your data actually contains distinct CD4 and CD8 T cell subtypes with different expression profiles, your pseudobulk analysis will average across these subtypes and may miss subtype-specific DE.

### Sensitivity to Sample Size

Pseudobulk analysis requires a sufficient number of biological samples per condition to estimate biological variability reliably. With small sample sizes, DESeq2 may have low statistical power, meaning that true DE genes may not reach significance. This is particularly problematic for genes with low expression or high variability between samples. Increasing the number of biological replicates is the most effective way to improve power, but this increases the cost and complexity of the experiment.

### Inability to Detect Compositional Changes

Pseudobulk DE analysis tests for changes in gene expression within a cell type, but it does not test for changes in the proportions of cell types between conditions. If a treatment causes an expansion of one cell type and a contraction of another, this compositional change will not be detected by DE analysis. Differential abundance testing methods, such as Milo, which assigns cells to partially overlapping neighborhoods on a k-nearest neighbor graph, can identify perturbations in cell abundance that are obscured by discretizing cells into clusters. If you are interested in both expression changes and abundance changes, you should run both types of analysis.

## Records and Measurements for Reproducible DE Analysis

Reproducibility is a central concern in bioinformatics, and DE analysis is no exception. Keeping detailed records of your analysis steps, parameters, and software versions allows others to reproduce your results and allows you to revisit your analysis with confidence.

### Documenting Software Versions and Parameters

Record the versions of all software packages used in your analysis, including R, Seurat, DESeq2, and any other tools. The Bioconductor project provides official documentation for its packages, and recording the session information ensures that your analysis can be reproduced with the same software environment. The `sessionInfo()` function in R prints the versions of all loaded packages, and you should save this output with your analysis results.

Also document the parameters you used at each step, including quality control thresholds, normalization methods, clustering resolution, and DE cutoffs. These parameters can substantially affect your results, and different choices may lead to different conclusions. A reproducible analysis records these choices explicitly.

### Using Workflow Management Tools

Workflow management tools can help ensure reproducibility by automating the execution of analysis steps and recording the parameters and software versions used. The nf-core community provides documentation for standardized bioinformatics pipelines that follow best practices for reproducibility. While nf-core pipelines are primarily designed for bulk RNA-seq and other genomics data, the principles of pipeline standardization and documentation apply equally to scRNA-seq analysis. Similarly, the Galaxy Training Network provides accessible workflow training that emphasizes reproducible analysis practices.

For scRNA-seq analysis, you can create your own workflow scripts that document each step of the analysis. These scripts should be version-controlled using Git, and you should record the commit hash associated with each analysis run. The Carpentries provides lessons on foundational computing skills, including shell, Git, and programming, that are valuable for managing bioinformatics workflows.

### Storing and Sharing Data

Store your raw sequencing data, processed count matrices, and analysis results in organized directories with clear naming conventions. Raw data should be deposited in appropriate repositories, such as those maintained by the National Center for Biotechnology Information (NCBI), to ensure long-term accessibility. Processed data and analysis scripts should be shared alongside publications to allow others to reproduce your results.

The NCBI provides a range of data resources for genomics research, including databases for raw sequencing data, processed expression data, and metadata. Depositing your data in these repositories ensures that it remains accessible to the research community and satisfies the data sharing requirements of most journals and funding agencies.

## Quality Control Checks Throughout the Workflow

Quality control is not a single step at the beginning of the analysis but an ongoing process that should be applied at multiple points throughout the workflow.

### Quality Control at the Single-Cell Level

Before pseudobulk aggregation, verify that your single-cell data passed appropriate quality control. Check the distributions of gene counts, UMI counts, and mitochondrial read fractions to ensure that low-quality cells were removed. Examine the clustering results to confirm that clusters are well separated and that marker genes are expressed in the expected cell types.

### Quality Control at the Pseudobulk Level

After creating pseudobulk matrices, check that the aggregated counts are sensible. Verify that the total counts per sample-cell type combination are reasonable and that there are no samples with extremely low counts that might indicate problems with cell capture or sequencing. Examine the correlation between biological replicates to ensure that replicates are more similar to each other than to samples from different conditions.

### Quality Control After DE Analysis

After running DESeq2, examine the diagnostic plots to check that the model fits are reasonable. The dispersion plot should show decreasing dispersion with increasing mean expression, and the MA plot should show a relatively symmetric distribution of fold changes around zero for non-significant genes. If the diagnostic plots look unusual, investigate potential problems with your data or model specification.

## Safety and Regulatory Context for Research Use

Differential expression analysis in scRNA-seq is a research tool, and the results should be interpreted within the context of the study design and limitations. The computational workflow described here does not involve direct manipulation of human subjects or animals, but the data may come from studies that involved such manipulation. Researchers should ensure that all data used in their analyses were collected under appropriate ethical approvals and that patient or animal identifiers are handled according to relevant regulations.

When DE analysis is used in translational research, such as identifying biomarkers for disease, the results should be validated in independent cohorts before any clinical application is considered. The machine learning framework for keloid biomarker discovery demonstrates how transcriptomic analysis can identify candidate biomarkers, but also highlights the importance of validation across independent datasets. Transcriptomic meta-analysis provides a framework for integrating gene expression studies across independent datasets to identify expression patterns that are reproducible, but heterogeneity in experimental design, sequencing platforms, and sample composition must be explicitly considered to avoid misleading conclusions.

## Professional Escalation Criteria

Knowing when to seek additional expertise can prevent wasted effort and incorrect conclusions. Consider consulting with a bioinformatics specialist or statistician in the following situations.

### When to Consult a Bioinformatics Specialist

If you are unfamiliar with R programming or the Seurat and DESeq2 packages, consider taking a training course or consulting with a colleague who has experience with these tools. The EMBL-EBI Training program provides learning pathways for bioinformatics data resources and practical analysis education, and the Galaxy Training Network offers accessible workflow training for researchers at all levels. Investing time in training can prevent many common errors and save time in the long run.

If your data has complex structure, such as multiple batches, nested experimental designs, or repeated measures from the same individuals, consult a statistician to ensure that your DE model is correctly specified. Incorrect model specification can lead to biased results and incorrect conclusions.

### When to Seek Help with Data Integration

If you are combining data from multiple experiments or sequencing platforms, data integration becomes a major challenge. Differences in experimental design, sequencing platforms, and sample composition introduce substantial heterogeneity that limits direct comparability between studies. If you are attempting a meta-analysis across independent datasets, consult with researchers who have experience with cross-study integration to ensure that your approach is appropriate.

### When Results Are Unexpected

If your DE analysis produces results that contradict known biology or that are difficult to interpret, investigate potential causes before accepting the results. Check for sample mislabeling, batch effects, cell type misannotation, and other technical issues. If the results remain puzzling after these checks, consult with colleagues who have domain expertise in your biological system.

## Frequently Asked Questions

### Why can I not run DESeq2 directly on single-cell counts?

DESeq2 assumes that each column in the count matrix represents an independent biological sample. When you provide single-cell counts, each cell is treated as an independent sample, which violates the assumption of independence because cells from the same biological sample share technical and biological variation. This pseudoreplication leads to inflated significance and excessive false positives. Pseudobulk aggregation combines cells from the same biological sample and cell type into a single count value, creating a count matrix where each column represents an independent biological replicate.

### How many biological replicates do I need for pseudobulk DE analysis?

The minimum recommended number of biological replicates per condition is three. With fewer than three replicates, DESeq2 cannot reliably estimate biological variability, and the results will be unstable. More replicates provide greater statistical power and more reliable dispersion estimates. The optimal number depends on the variability of your system, the magnitude of the effect you are trying to detect, and the depth of sequencing.

### What is the difference between differential expression and differential abundance?

Differential expression tests whether the average expression level of a gene differs between conditions within a given cell type. Differential abundance tests whether the proportion of cells belonging to a particular cell type or state differs between conditions. These are complementary analyses that answer different biological questions. Pseudobulk DE analysis addresses differential expression, while methods such as Milo address differential abundance.

### Should I use normalized or raw counts for pseudobulk aggregation?

You should use raw counts for pseudobulk aggregation. DESeq2 performs its own normalization internally through the estimation of size factors, and providing normalized data can distort the count distribution and lead to incorrect results. The raw counts should be summed across cells within each sample-cell type combination, and the resulting matrix should be provided to DESeq2 without additional normalization.

### How do I choose the thresholds for significant genes?

The choice of thresholds depends on your biological question and the power of your experiment. A common approach is to use an adjusted p-value cutoff of 0.05 and a log2 fold change cutoff of 1, corresponding to a 2-fold change. However, you may choose more stringent thresholds if you are prioritizing a small number of genes for validation, or more lenient thresholds if you are exploring the data and plan to validate results with additional experiments.

### Can I use this workflow for single-nucleus RNA-seq data?

Yes, the pseudobulk DE workflow applies equally to single-nucleus RNA-seq data. The upstream processing steps, including quality control, normalization, and clustering, may need to be adjusted for the characteristics of nuclear transcripts, but the pseudobulk aggregation and DESeq2 analysis steps are the same. Ensure that your quality control thresholds are appropriate for your data modality.

### What should I do if my cell types are not well separated in clustering?

If your clusters are poorly separated, your cell type annotations may be unreliable, and your pseudobulk DE results may be compromised. Consider using a higher clustering resolution to identify more distinct subpopulations, or use reference-based label transfer to annotate your clusters based on known cell type markers. Validate your annotations using multiple marker genes and examine the expression of these markers in your clusters.

### How do I account for multiple cell types in my DE analysis?

You should perform DE analysis separately for each cell type. Create a pseudobulk matrix for each cell type, run DESeq2 on each matrix, and interpret the results in the context of each cell type. This approach allows you to identify cell-type-specific DE genes, which is one of the main advantages of single-cell over bulk RNA-seq analysis.

## Related Bioinformatics Guides

- [Proteomics Data Analysis in R: A Practical Workflow for Differential Expression and Visualization](/knowledge/bioinformatics/proteomics-data-analysis-in-r-a-practical-workflow-for-differential-expression-and-visualization)
- [RNA Sequencing Data Analysis: From Raw Reads to Differential Expression](/knowledge/bioinformatics/rna-sequencing-data-analysis-from-raw-reads-to-differential-expression)
- [RNA-Seq Data Analysis Workflow: From Raw Reads to Insights](/knowledge/bioinformatics/rna-seq-data-analysis-workflow-from-raw-reads-to-insights)
- [Single-Cell Sequencing Workflow: From Sample Preparation to Data Analysis](/knowledge/bioinformatics/single-cell-sequencing-workflow-from-sample-preparation-to-data-analysis)
- [Single-Cell RNA Sequencing Quality Control: A Practical Guide to Filtering and Metrics](/knowledge/bioinformatics/single-cell-rna-sequencing-quality-control-a-practical-guide-to-filtering-and-metrics)

## Related Clinical & Scientific Guides

* [A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data](/knowledge/bioinformatics/a-practical-guide-to-detecting-antimicrobial-resistance-genes-in-shotgun-metagenomic-data)
* [Computational Immunology: Modeling the Immune System](/knowledge/bioinformatics/computational-immunology-modeling-the-immune-system)
* [How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices](/knowledge/bioinformatics/how-to-set-hard-filters-for-germline-variant-calling-a-practical-guide-to-gatk-best-practices)


## References and Further Reading

- [NCBI Data Resources](https://www.ncbi.nlm.nih.gov/). National Center for Biotechnology Information.
- [EMBL-EBI Training](https://www.ebi.ac.uk/training). European Bioinformatics Institute.
- [Bioconductor](https://bioconductor.org/). Bioconductor Project.
- [Galaxy Training Network](https://training.galaxyproject.org/). Galaxy Project.
- [nf-core Documentation](https://nf-co.re/docs). nf-core.
- [The Carpentries Lessons](https://carpentries.org/lessons). The Carpentries.
- [Single-Cell RNA Sequencing Analysis: A Step-by-Step Overview.](https://pubmed.ncbi.nlm.nih.gov/33835452). Methods in molecular biology (Clifton, N.J.), 2021.
- [Best practices on the differential expression analysis of multi-species RNA-seq.](https://pubmed.ncbi.nlm.nih.gov/33926528). Genome biology, 2021.
- [Differential abundance testing on single-cell data using k-nearest neighbor graphs.](https://pubmed.ncbi.nlm.nih.gov/34594043). Nature biotechnology, 2022.
- [Plant Single-Cell/Nucleus RNA-seq Workflow.](https://pubmed.ncbi.nlm.nih.gov/36495448). Methods in molecular biology (Clifton, N.J.), 2023.
- [Multiplexed droplet single-cell RNA-sequencing using natural genetic variation.](https://pubmed.ncbi.nlm.nih.gov/29227470). Nature biotechnology, 2018.
- [Guidelines for bioinformatics of single-cell sequencing data analysis in Alzheimer's disease: review, recommendation, implementation and application.](https://pubmed.ncbi.nlm.nih.gov/35236372). Molecular neurodegeneration, 2022.
- [Computational Analysis of Single-Cell RNA-Seq Data.](https://pubmed.ncbi.nlm.nih.gov/36264495). Methods in molecular biology (Clifton, N.J.), 2023.
- [Analysis of a Single Cell RNA-seq Workflow by Random Matrix Theory Methods.](https://pubmed.ncbi.nlm.nih.gov/39585539). Bulletin of mathematical biology, 2024.
- [Multisite Assessment of Methods for Cell Preservation Upstream of Single-Cell RNA Sequencing.](https://doi.org/10.7171/001c.162768). 2026.
- [A Robust Machine Learning Framework for Keloid Biomarker Discovery Beyond Differential Expression](https://doi.org/10.64898/2026.06.24.734231). 2026.
- [Transcriptomic Meta-Analysis as a Framework for Robust Cross-Study Biological Inference.](https://doi.org/10.3390/ijms27114674). 2026.
- [Differential detection workflows for multi-sample single-cell RNA-seq data](https://doi.org/10.1101/2023.12.17.572043). bioRxiv, 2023.
- [Spatial Transcriptomic Modeling of Vascular Remodeling in Aortic Aneurysm Using Integrated Single-Cell RNA Sequencing Analysis.](https://doi.org/10.1016/j.slast.2025.100377). SLAS technology, 2025.
- [Single-cell RNA sequencing data analysis of the inner ear in gentamicin-treated mice via intraperitoneal injection](https://doi.org/10.1515/med-2025-1242). Open Medicine, 2025.
- [Automation of RNA-Seq Sample Preparation and Miniaturized Parallel Bioreactors Enable High-Throughput Differential Gene Expression Studies](https://doi.org/10.3390/microorganisms13040849). Microorganisms, 2025.
- [Single-Cell RNA Sequencing Reveals Cardiac Fibroblast-Specific Transcriptomic Changes in Dilated Cardiomyopathy](https://doi.org/10.3390/cells13090752). Cells, 2024.
- [Differential detection workflows for multi-sample single-cell RNA-seq data](https://doi.org/10.1186/s12864-025-12102-x). BMC Genomics, 2025.
- [Data analysis in single-cell RNA-Seq](https://doi.org/10.1016/B978-0-12-814919-5.00019-1). Single Cell Omics Volume 1 Technological Advances and Applications, 2019.
- [FastCAR: fast correction for ambient RNA to facilitate differential gene expression analysis in single-cell RNA-sequencing datasets](https://doi.org/10.1186/s12864-023-09822-3). BMC Genomics, 2023.

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