18. Differential gene expression analysis

   Key takeaways

In scRNA-seq, there are two main approaches to identify differentially expressed genes: the sample-level view and the cell-level view. The sample-level approach has been shown to outperform the cell-level approach for scRNA-seq data.

Motivation

In the sample-level view, a pseudobulk sample is created for each patient and cell type by aggregating counts (e.g., summing across cells). These pseudobulk profiles can then be analyzed using tools such as edgeR or DESeq2, which were originally developed for bulk RNA-seq data.

Pseudobulking

Explore your data before modeling. Performing PCA on pseudobulk samples can reveal major sources of variation that should be accounted for in your design matrix.

Variability Exploration

Remove lowly expressed genes separately for each cell type, as different cell populations express distinct gene sets.

Feature selection

If your dataset contains different subgroups, construct a design matrix that includes all relevant covariates. You can then define contrasts to test the specific comparisons of interest.

Multiple cell types or groups

   Environment setup

Steps

1. Install conda:

   • Before creating the environment, ensure that conda is installed on your system.

2. Save the yml content:

   • Copy the content from the yml tab into a file named environment.yml.

3. Create the environment:

   • Open a terminal or command prompt.

   • Run the following command:

     conda env create -f environment.yml
     

4. Activate the environment:

   • After the environment is created, activate it using:

     conda activate <environment_name>
     

   • Replace <environment_name> with the name specified in the environment.yml file. In the yml file it will look like this:

     name: <environment_name>
     

5. Verify the installation:

   • Check that the environment was created successfully by running:

     conda env list
     

yml

name: differential-gene-expression
channels:
  - conda-forge
  - bioconda
dependencies:
  - conda-forge::decoupler-py=2.1.6
  - conda-forge::pertpy=1.0.6
  - bioconda::pydeseq2=0.5.4
  - conda-forge::python=3.13.13
  - conda-forge::scanpy=1.12.1

   Get data and notebooks

This book uses lamindb to store, share, and load datasets and notebooks using the theislab/sc-best-practices instance. We acknowledge free hosting from Lamin Labs.

1. Install lamindb

   • Install the lamindb Python package:

   pip install lamindb
   

2. Optionally create a lamin account

   • Sign up and log in following the instructions

3. Verify your setup

   • Run the lamin connect command:

   import lamindb as ln
   
   ln.Artifact.connect("theislab/sc-best-practices").df()
   

   You should now see up to 100 of the stored datasets.

4. Accessing datasets (Artifacts)

   • Search for the datasets on the Artifacts page

   • Load an Artifact and the corresponding object:

   import lamindb as ln
   af = ln.Artifact.connect("theislab/sc-best-practices").get(key="key_of_dataset", is_latest=True)
   obj = af.load()
   

   The object is now accessible in memory and is ready for analysis. Adapt the lamindb.Artifact.connect("theislab/sc-best-practices").get("SOMEIDXXXX") suffix to get respective versions.

5. Accessing notebooks (Transforms)

   • Search for the notebook on the Transforms page

   • Load the notebook:

   lamin load <notebook url>
   

   which will download the notebook to the current working directory. Analogously to Artifacts, you can adapt the suffix ID to get older versions.

18.1. Motivation

This chapter is a more detailed continuation of the annotation subchapter which already introduced differential gene expression (DGE) as a tool to annotate clusters with cell types. Here, we focus on more advanced use-cases of differential gene expression testing in complex experimental designs, which involve one or more conditions such as diseases, genetic knockouts or drugs. In such cases we are commonly interested in the magnitude and significance of differences in gene expression patterns between the condition of interest and a reference. This reference can be everything but is commonly a healthy sample. This statistical test can be applied to arbitrary groups, but in the case of single-cell RNA-Seq is commonly applied on the cell type level.

Fig. 18.1 DGE analysis attempts to infer genes that are statistically significantly over- or underexpressed between any compared groups (commonly between healthy and condition per cell type).

The outcome of such an analysis could be genesets which effect and potentially explain any observed phenotypes. These genesets can then be examined more closely with respect to, for example, affected pathways or induced cell-cell communication changes.

A differential gene expression test usually returns the log2 fold-change and the adjusted p-value per compared genes per compared conditions. This list can then be sorted by p-value and investigated in more detail.

The popular student’s t-test is one way of conducting such a test. However, it fails to take several single-cell RNA-seq peculiarities into account such as the excess number of zeros originating from dropouts or the need for complex experimental designs. More specifically, very rarely does one have sufficient sample numbers to accurately estimate the variance without pooling information across genes. Moreover, raw counts are never an absolute measurement of expression for a specific gene within a given sample. The actual read number per gene depends on the efficiency of the library preparation, the amount of contamination from non-coding transcripts and the sequencing depth. It does therefore lack in both, sensitivity and specificity for single-cell RNA-seq, let alone experimental design flexibility.

As a result, DGE testing is a classic bioinformatics problem which has been tackled by many tools already. Generally, the problem is currently being approached from two views, the sample-level view where expression is aggregated to create “pseudobulks” and then analysed with methods originally designed for bulk expression samples such as edgeR [Robinson et al., 2010] or DEseq2 [Love et al., 2014] and the cell-level view where cells are modeled individually using generalized mixed effect models such as MAST [Finak et al., 2015] or glmmTMB [Brooks et al., 2017]. The consensus and robustness across datasets for DGE tools is low [Das et al., 2021, Wang et al., 2019]. As previously described, although single-cell data contains technical noise artifacts such as dropout, zero-inflation and high cell-to-cell variability [Hicks et al., 2017, Luecken and Theis, 2019, Vallejos et al., 2017], methods designed for bulk RNA-seq data performed favorably compared to methods explicitly designed for scRNA-seq data [Das et al., 2021, Jaakkola et al., 2016, Soneson and Robinson, 2018, Squair et al., 2021]. Single-cell specific methods were found to be especially prone to wrongly labeling highly expressed genes as differentially expressed.

A recent study highlighted the issue of pseudoreplication where inferential statistics is applied to biological replicates which are not statistically independent. Failing to account for the inherent correlation of replicates (cells from the same individual) inflates the false discovery rate (FDR) [Junttila et al., 2022, Squair et al., 2021, Zimmerman et al., 2021]. Therefore, batch effect correction or the aggregation of cell-type-specific expression values within an individual through either a sum, mean or random effect per individual, that is pseudobulk generation, should be applied prior to DGE analysis to account for within-sample correlations [Zimmerman et al., 2021]. Generally, both, pseudobulk methods with sum aggregation such as edgeR, DESeq2, or Limma [Ritchie et al., 2015] and mixed models such as MAST with random effect setting were found to be superior compared to naive methods, such as the popular Wilcoxon rank-sum test or Seurat’s [Hao et al., 2021] latent models, which do not account for them [Junttila et al., 2022].

In matters arising from the Zimmerman paper, Murphy et al. critically examined the Zimmerman benchmarking strategy and improved it [Murphy and Skene, 2022]. They came to the conclusion that pseudobulk methods perform best but whether sum or mean aggregation works better requires further investigation.

Hence, in this notebook, we demonstrate how to perform DGE analysis using the pseudobulk approach. We choose DESeq2 due to its Python implementation, PyDESeq2. However, the code can be easily adapted for use with edgeR by following the pertpy tutorial on differential gene expression. For this chapter, we combined parts of the tutorials from decoupler and pertpy. Please refer to their documentation for more details.

18.2. Environment setup

import decoupler as dc
import lamindb as ln
import numpy as np
import pandas as pd
import pertpy as pt
import scanpy as sc

ln.track()

Show code cell output

Hide code cell output

→ loaded Transform('zt7x1ebMcoIA0000', key='differential_gene_expression.ipynb'), re-started Run('Oh5MpUtPYV0UPunU') at 2026-05-12 06:49:26 UTC
→ notebook imports: decoupler==2.1.6 lamindb-core==2.4.2 numpy==2.4.3 pandas==2.3.3 pertpy==1.0.6 scanpy==1.12.1
• recommendation: to identify the notebook across renames, pass the uid: ln.track("zt7x1ebMcoIA")

18.3. Preparing the dataset

We will use the Kang dataset, which is a 10x droplet-based scRNA-seq peripheral blood mononuclear cell (PBMC) data from 8 Lupus patients before and after 6h-treatment with INF-β (16 samples in total) [Kang et al., 2018]. Interferon beta is used in the form of natural fibroblast or recombinant preparations (interferon beta-1a and interferon beta-1b) and exerts antiviral and antiproliferative properties similar to those of interferon alpha. Interferon beta has been approved for the treatment of relapsing–remitting multiple sclerosis and secondary progressive multiple sclerosis.

First, we load the full dataset.

adata = ln.Artifact.get(
    key="conditions/differential_gene_expression.h5ad",
).load()
adata

AnnData object with n_obs × n_vars = 24673 × 15706
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters'
    var: 'name'
    obsm: 'X_pca', 'X_umap'

We will need label (which contains the condition label), replicate (patient id) and cell_type columns of the .obs.

adata.obs[:5]

	nCount_RNA	nFeature_RNA	tsne1	tsne2	label	cluster	cell_type	replicate	nCount_SCT	nFeature_SCT	integrated_snn_res.0.4	seurat_clusters
index
AAACATACATTTCC-1	3017.0	877	-27.640373	14.966629	ctrl	9	CD14+ Monocytes	patient_1016	1704.0	711	1	1
AAACATACCAGAAA-1	2481.0	713	-27.493646	28.924885	ctrl	9	CD14+ Monocytes	patient_1256	1614.0	662	1	1
AAACATACCATGCA-1	703.0	337	-10.468194	-5.984389	ctrl	3	CD4 T cells	patient_1488	908.0	337	6	6
AAACATACCTCGCT-1	3420.0	850	-24.367997	20.429285	ctrl	9	CD14+ Monocytes	patient_1256	1738.0	653	1	1
AAACATACCTGGTA-1	3158.0	1111	27.952170	24.159738	ctrl	4	Dendritic cells	patient_1039	1857.0	928	12	12

We will need to work with raw counts so we check that .X indeed contains raw counts and put them into the counts layer of our AnnData object.

X = adata.X.data
np.array_equal(X, np.round(X))

True

adata.layers["counts"] = adata.X.copy()

We have 8 control and 8 disease patients.

print(len(adata[adata.obs["label"] == "ctrl"].obs["replicate"].cat.categories))
print(len(adata[adata.obs["label"] == "stim"].obs["replicate"].cat.categories))

8
8

We filter cells which have less than 200 genes and genes which were found in less than 3 cells for a rudimentary quality control.

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata.shape

(24562, 15701)

18.4. Pseudobulking

Since PyDESeq2 was introduced as a method for DGE analysis for bulk data, we first need to create pseudobulk samples from our single-cell dataset. For each patient we create one pseudobulk sample per cell type by aggregating the cells from each subpopulation and taking the sum of gene expression counts. Regardless of whether we want to run the analysis only on a few cell subpopulations and fit a model for each one of them separately or fit one model for all of them, we first need to prepare the data.

Since we need to create pseudobulks for each patient-condition combination, we first need to create such a column by concatenating replicate and label.

adata.obs["sample"] = pd.Categorical(
    f"{rep}_{l}" for rep, l in zip(
        adata.obs["replicate"],
        adata.obs["label"],
        strict=False
    )
)

Next, we generate pseudobulk samples using decoupler. With 8 patients, 8 cell types, and 2 conditions (ctrl and stim), this yields a total of 128 pseudobulks (8 x 8 x 2).

adata_pb = dc.pp.pseudobulk(
    adata, sample_col="sample",
    groups_col="cell_type",
    layer="counts",
    mode="sum"
)
adata_pb

AnnData object with n_obs × n_vars = 128 × 15701
    obs: 'sample', 'cell_type', 'label', 'replicate', 'psbulk_cells', 'psbulk_counts'
    var: 'name', 'n_cells'
    layers: 'psbulk_props'

Low-quality samples can be filtered using two criteria: the number of cells (psbulk_cells) and the total count sum (psbulk_counts). These metrics can be visualized in the following plot to guide filtering decisions.

dc.pl.filter_samples(
    adata=adata_pb,
    groupby=["label", "cell_type", "replicate"],
    min_cells=10,
    min_counts=1000,
    figsize=(5, 8),
)

While thresholds are dataset-specific and arbitrary, a commonly used guideline is to retain samples with a minimum of 10 cells and 1,000 total counts. Applying these thresholds (indicated by dashed lines) restricts the dataset to pseudobulks in the upper-right quadrant. The filtering step can be carried out with decoupler.pp.filter_samples().

dc.pp.filter_samples(adata_pb, min_cells=10,min_counts=1000)

After filtering out low-quality samples, we can visualize the remaining profiles. All megakaryocyte pseudobulk samples failed quality control, and two CD8 T cell samples were removed.

dc.pl.obsbar(adata=adata_pb, y="cell_type", hue="label", figsize=(6, 3))

18.5. Variability Exploration

The validity of DGE results highly depends on the capture of the major axis of variations in the statistical model. Intermediate data exploration steps such as principal component analysis (PCA) or multidimensional scaling (MDS) on pseudobulk samples allow for the identification of the sources of variation and thus can guide the construction of corresponding design and contrast matrices that model the data [Law et al., 2020].

Failing to account for multiple sources of biological variability for experiments which include biological replicates will inflate the FDR [Lähnemann et al., 2020, Thurman et al., 2021]. While increasing the number of cells per individual increases the precision, it has a limited effect on the power for the detection of differences across individuals. Therefore, the best way to increase statistical power is to increase the number of independent experimental samples [Zimmerman et al., 2021].

Since our data has already been generated, we cannot further increase the number of independent experimental samples. Nevertheless, we will now explore our data to determine the major axes of variation to properly generate our design matrices.

We perform very basic exploratory data analysis on the generated pseudo-replicates to identify potential outliers among patients or pseudobulks. These outliers can then be excluded to avoid biasing the differential expression results. We store the raw counts in the 'counts' layer, then normalize and scale them before computing PCA. Afterwards, we revert to the raw counts for downstream analysis.

adata_pb.layers['counts'] = adata_pb.X.copy()

sc.pp.normalize_total(adata_pb, target_sum=1e6)
sc.pp.log1p(adata_pb)
sc.pp.scale(adata_pb, max_value=10)
sc.tl.pca(adata_pb)

dc.pp.swap_layer(adata=adata_pb, key="counts", inplace=True)

We also include psbulk_counts_log and psbulk_cells_log because taking the log makes library sizes (psbulk_counts and psbulk_cells) less skewed and allows more reliable detection of correlations with principal components.

adata_pb.obs["psbulk_counts_log"] = np.log(adata_pb.obs["psbulk_counts"])
adata_pb.obs["psbulk_cells_log"] = np.log(adata_pb.obs["psbulk_cells"])

dc.tl.rankby_obsm(adata_pb, key="X_pca")
sc.pl.pca_variance_ratio(adata_pb)
dc.pl.obsm(
    adata=adata_pb,
    return_fig=True,
    nvar=5,
    titles=["PC scores", "Adjusted p-values"],
    figsize=(10, 5)
)

In this dataset, PC1 appears to explain the largest proportion of variance because it shows the highest variance ratio of around 0.13. It is associated with the metadata variables cell type and pseudobulk cells. Metadata variables associated with PCs that capture a substantial amount of variance are important and should be accounted for as relevant covariates (e.g., include in the design matrix) in downstream differential expression analysis when possible. The principal components can also be directly visualized, colored by these metadata variables.

adata_pb.obs = adata_pb.obs.sort_index(axis=1)
sc.pl.pca(adata_pb, color=adata_pb.obs, ncols=2, size=300)

We observe separation of cell types on the PCA plots as well as the separation into stimulated and unstimulated cells. For pseudobulk cells and counts, we can also observe some clustering, with high values appearing in the top left and top right regions. Looking more closely, we see that these clusters of high values partially overlap with specific cell type clusters: the top right cluster overlaps with CD4 T cells, while the top left cluster overlaps with CD14+ monocytes. This suggests that the number of cells and counts assigned to a pseudospot may depend on the pseudospot’s cell type. The covariates replicate and sample do not seem to be clearly correlated with the PCA components so we do not include any of them in our design matrix.

18.6. One cell type or group

If you already know which cell types to focus on, subset them at this step to reduce confounding factors. Otherwise, you can first analyze the full dataset to identify the most affected cell types, then return to this step.

If variability exploration showed that cell types are associated with PCs, it may help to rerun it on your cell type subsets to see if covariate effects disappear. In our case, the association with pseudobulk cells disappears in CD14+ monocytes. That’s why we won’t include it into our first design matrix.

18.6.1. Feature selection

In addition to filtering low-quality samples, lowly or noisily expressed genes can also be filtered prior to DGE analysis. This step should be performed at the cell type level, as different cell types may express distinct sets of genes. We run the pipeline on CD14+ monocytes subset of the data, as it was shown in the paper that the highest number of differentially expressed genes was identified in this subpopulation.

adata_mono = adata_pb[adata_pb.obs["cell_type"] == "CD14+ Monocytes"].copy()
adata_mono.shape

(16, 15701)

Two strategies are used to filter genes:

1. decoupler.pp.filter_by_expr(): Retains genes with a minimum total number of reads across all samples (min_total_count) and a minimum number of counts in a given number of samples (min_count). This approach was introduced in edgeR [Robinson et al., 2010].

2. decoupler.pp.filter_by_prop(): Retains genes that are expressed in at least a specified proportion of cells (min_prop) across a minimum number of samples (min_smpls). The number of retained genes can be visualized, and the filtering parameters can be adjusted interactively.

dc.pl.filter_by_expr(
    adata=adata_mono,
    group="label",
    min_count=10,
    min_total_count=15,
    large_n=10,
    min_prop=0.7,
)
dc.pl.filter_by_prop(
    adata=adata_mono,
    min_prop=0.1,
    min_smpls=2,
)

The top plot displays gene frequencies based on the filter_by_expr metrics, while the bottom plot corresponds to filter_by_prop. Dashed lines indicate the current threshold values. In the top plot, only genes in the upper-right quadrant are retained, in the bottom plot, only those to the right of the vertical line are kept.

Although filtering thresholds are arbitrary, a common heuristic is to look for bimodal distributions and set thresholds that separate low-quality genes from the rest. In this example, the default parameters retain a substantial number of genes while removing potentially noisy ones.

Once the threshold parameters are set, the actual gene filtering can be performed by simply changing pl to pp.

dc.pp.filter_by_expr(
    adata=adata_mono,
    group="label",
    min_count=10,
    min_total_count=15,
    large_n=10,
    min_prop=0.7,
)
dc.pp.filter_by_prop(
    adata=adata_mono,
    min_prop=0.1,
    min_smpls=2,
)
adata_mono.shape

(16, 2345)

18.6.2. Differential expression testing with PyDESeq2

We now switch to the pertpy package, as it makes it easier to use its built-in plotting functions. However, DGE analysis can also be performed manually using PyDESeq2. The interface of PyDESeq2 in pertpy is very similar.

At the next step, it is important to define an appropriate design matrix for the model. For this purpose, we strongly recommend reading this guide, which provides a helpful introduction to design matrices.

Understand the design syntax

We combine all the variables we would like to include in our model with a +. The variable names correspond to the column names ins .obs. Interactions between variables can be added to the model with :. For example, the model y ~ a + b + a:b includes the effect of a and b on y as well as how the effect of a on y changes depending on b and vice versa. In our context, y is the observed gene expression and a and b could be label and cell_type. Now it also becomes clear where the ~ comes from: It separates the two sides of a formula. So just add ~ at the beginning of your design. Besides that, it might be useful to know that a * b is shorthand for a + b + a:b. In most cases, this is all you need. For a deeper dive, we recommend reading the formulaic and patsy documentation.

pds2 = pt.tl.PyDESeq2(adata=adata_mono,design='~ label')

pds2.fit()

Show code cell output

Hide code cell output

Using None as control genes, passed at DeseqDataSet initialization

Fitting size factors...
... done in 0.00 seconds.

Fitting dispersions...
... done in 0.32 seconds.

Fitting dispersion trend curve...
... done in 0.03 seconds.

Fitting MAP dispersions...
... done in 0.36 seconds.

Fitting LFCs...
... done in 0.22 seconds.

Calculating cook's distance...
... done in 0.00 seconds.

Replacing 2 outlier genes.

Fitting dispersions...
... done in 0.01 seconds.

Fitting MAP dispersions...
... done in 0.01 seconds.

Fitting LFCs...
... done in 0.01 seconds.

res_df = pds2.test_contrasts(pds2.contrast(
    column="label",
    baseline="ctrl",
    group_to_compare="stim")
)

Show code cell output

Hide code cell output

Running Wald tests...
... done in 0.10 seconds.

Log2 fold change & Wald test p-value, contrast vector: [0. 1.]
              baseMean  log2FoldChange     lfcSE       stat        pvalue  \
index                                                                       
HES4        179.827743        6.600825  0.379008  17.416074  6.230913e-68   
ISG15     19252.031373        7.076825  0.420073  16.846656  1.110180e-63   
SDF4         28.927758       -0.757192  0.178380  -4.244825  2.187642e-05   
UBE2J2       25.529327       -0.889011  0.187880  -4.731802  2.225355e-06   
AURKAIP1    115.547324       -0.572923  0.097563  -5.872361  4.296310e-09   
...                ...             ...       ...        ...           ...   
CSTB       1201.426144        0.864196  0.189358   4.563812  5.023307e-06   
SUMO3        72.227191       -0.429374  0.109306  -3.928191  8.558717e-05   
PTTG1IP      54.543949        0.213316  0.135542   1.573794  1.155351e-01   
ITGB2        41.979324        1.341235  0.240249   5.582690  2.368272e-08   
PRMT2        70.413437        0.566704  0.149073   3.801520  1.438112e-04   

                  padj  
index                   
HES4      2.981937e-66  
ISG15     4.412494e-62  
SDF4      4.560019e-05  
UBE2J2    5.234160e-06  
AURKAIP1  1.303344e-08  
...                ...  
CSTB      1.133749e-05  
SUMO3     1.671123e-04  
PTTG1IP   1.438820e-01  
ITGB2     6.699154e-08  
PRMT2     2.748470e-04  

[2345 rows x 6 columns]

res_df.head(10)

	variable	baseMean	log_fc	lfcSE	stat	p_value	adj_p_value	contrast
0	IFITM2	491.867844	4.670296	0.151406	30.846095	6.319564e-209	1.481938e-205	None
1	IL1RN	900.090770	6.466407	0.223553	28.925558	5.697320e-184	6.680107e-181	None
2	NT5C3A	221.226625	5.583271	0.193526	28.850299	5.023241e-183	3.926500e-180	None
3	SSB	365.880933	3.336528	0.116120	28.733517	1.455497e-181	8.532854e-179	None
4	RABGAP1L	211.055152	5.994200	0.230397	26.016826	3.194924e-149	1.498419e-146	None
5	RTCB	152.556535	3.431247	0.134175	25.572931	3.052671e-144	1.193086e-141	None
6	DYNLT1	948.183449	2.609925	0.102510	25.460223	5.439750e-143	1.822316e-140	None
7	CCL8	7259.221485	9.389468	0.369940	25.381043	4.083958e-142	1.197110e-139	None
8	CXCL11	1433.690007	9.021039	0.357845	25.209353	3.163259e-140	8.242047e-138	None
9	DDX58	201.134679	4.946996	0.200043	24.729696	5.127619e-135	1.202427e-132	None

Let’s take a look at the columns in the results table:

• variable. gene name

• baseMean: mean of normalized counts for all samples

• log_fc: log2 fold change

• lfcSE: standard error

• stat: Wald statistic

• p_value: Wald test p-value

• adj_p_value: Benjamini-Hochberg adjusted p-values

The rows are sorted by adj_p_value. The log_fc always depends on the defined baseline. A positive log_fc means gene expression is higher in the group of interest compared to the baseline group, while negative values indicate lower expression. In this case, IFITM2 increased by around 4.7 in stimulated compared to control CD14+ monocytes. In contrast, VCAN decreased by around 4.4 (see the plot below).

Next, we visualize the results.

pds2.plot_volcano(res_df, log2fc_thresh=0)

The higher a point is on the y-axis, the smaller the adj_p_value. The further a point is from 0 on the x-axis, the larger the change in expression. Overall, this means that genes in the upper left and upper right regions of the plot show the strongest changes in gene expression and, at the same time, are also the least likely to have changes that occurred by chance.

We can also visualize our results for different subgroups. In this case, we plot the results for individual patients.

pds2.plot_paired(
    adata_mono,
    results_df=res_df,
    n_top_vars=4,
    groupby="label",
    pairedby="replicate"
)

We can see that the CD14+ monocytes from patient 1015 appear to be highly responsive to the treatment.

18.7. Multiple cell types or groups

If you’re unsure which cell type will be most affected by a treatment, a good starting point is to fit a model to the entire dataset. Afterwards, you can return to the analysis steps from the first part of this chapter and focus on the cell type of interest.

In addition, the following steps demonstrate how to handle analyses involving multiple groups within the data. You may already know which cell type you want to examine more closely, but your dataset might include multiple groups (e.g., sex, responder vs. non-responder status or age groups). Then, this section provides an initial glimpse of the types of questions you can explore for a single cell type across multiple groups.

Since we do not have particularly interesting subgroups within our CD14+ monocyte subset (adata_mono), we will instead analyze multiple cell types together. Because this analysis includes several cell types simultaneously, we will skip the feature selection step at this stage. Feature selection is better performed at the level of individual cell types, as different cell types can express distinct gene sets.

We will begin by constructing a simple design matrix.

pds2 = pt.tl.PyDESeq2(adata=adata_pb, design="~ label")

pds2.fit()

Show code cell output

Hide code cell output

Fitting size factors...
... done in 0.03 seconds.

Using None as control genes, passed at DeseqDataSet initialization

Fitting dispersions...
... done in 2.05 seconds.

Fitting dispersion trend curve...
... done in 0.18 seconds.

Fitting MAP dispersions...
... done in 2.69 seconds.

Fitting LFCs...
... done in 2.39 seconds.

Calculating cook's distance...
... done in 0.08 seconds.

Replacing 36 outlier genes.

Fitting dispersions...
... done in 0.02 seconds.

Fitting MAP dispersions...
... done in 0.02 seconds.

Fitting LFCs...
... done in 0.05 seconds.

res_df = pds2.compare_groups(
    adata_pb,
    column="label",
    baseline="ctrl",
    groups_to_compare=["stim"]
)

Show code cell output

Hide code cell output

Fitting size factors...
... done in 0.03 seconds.

Using None as control genes, passed at DeseqDataSet initialization

Fitting dispersions...
... done in 2.24 seconds.

Fitting dispersion trend curve...
... done in 0.19 seconds.

Fitting MAP dispersions...
... done in 2.59 seconds.

Fitting LFCs...
... done in 1.70 seconds.

Calculating cook's distance...
... done in 0.08 seconds.

Replacing 36 outlier genes.

Fitting dispersions...
... done in 0.02 seconds.

Fitting MAP dispersions...
... done in 0.01 seconds.

Fitting LFCs...
... done in 0.01 seconds.

Running Wald tests...

Log2 fold change & Wald test p-value, contrast vector: [0. 1.]
                baseMean  log2FoldChange     lfcSE      stat    pvalue  \
index                                                                    
AL627309.1      0.010607        0.043299  2.310776  0.018738  0.985050   
RP11-206L10.2   0.030188       -0.014207  2.115790 -0.006715  0.994642   
RP11-206L10.9   0.047391        0.120755  1.838275  0.065689  0.947625   
FAM87B          0.053617       -0.028965  2.603665 -0.011125  0.991124   
LINC00115       0.389827       -0.221986  0.420872 -0.527442  0.597886   
...                  ...             ...       ...       ...       ...   
C21orf58        0.121437       -0.262831  0.920014 -0.285682  0.775122   
PCNT            0.676612       -0.112711  0.333573 -0.337892  0.735445   
DIP2A           1.532521       -0.289122  0.273588 -1.056781  0.290612   
S100B           0.638731       -0.509243  0.795387 -0.640245  0.522013   
PRMT2          24.013222        0.241891  0.132709  1.822721  0.068346   

                   padj  
index                    
AL627309.1          NaN  
RP11-206L10.2       NaN  
RP11-206L10.9       NaN  
FAM87B              NaN  
LINC00115           NaN  
...                 ...  
C21orf58            NaN  
PCNT           0.841589  
DIP2A          0.487455  
S100B          0.694474  
PRMT2          0.182654  

[15701 rows x 6 columns]

... done in 0.49 seconds.

pds2.plot_paired(
    adata_pb,
    results_df=res_df,
    n_top_vars=4,
    groupby="label",
    pairedby="cell_type"
)

• Performing pseudobulk for paired samples

We can see that the top four differentially expressed genes show their strongest expression changes not only in CD14+ monocytes but also in CD4 T cells.

Now we could ask: Is the gene expression difference between ctrl and stim different between CD14+ Monocytes and CD4 T cells? In other words, could the gene expression changes depend on the cell type?

To do so, we create a more complex desgin matrix and then use contrasts to specify the conditions of interest.

pds2 = pt.tl.PyDESeq2(adata=adata_pb, design="~ cell_type * label")

pds2.fit()

Show code cell output

Hide code cell output

Fitting size factors...
... done in 0.03 seconds.

Using None as control genes, passed at DeseqDataSet initialization

Fitting dispersions...
... done in 2.20 seconds.

Fitting dispersion trend curve...
... done in 0.19 seconds.

Fitting MAP dispersions...
... done in 2.85 seconds.

Fitting LFCs...
... done in 3.00 seconds.

Calculating cook's distance...
... done in 0.05 seconds.

Replacing 3 outlier genes.

Fitting dispersions...
... done in 0.01 seconds.

Fitting MAP dispersions...
... done in 0.01 seconds.

Fitting LFCs...
... done in 0.01 seconds.

interaction_contrast = (
    pds2.cond(cell_type="CD14+ Monocytes", label="stim") -
    pds2.cond(cell_type="CD14+ Monocytes", label="ctrl")
) - (
    pds2.cond(cell_type="CD4 T cells", label="stim") -
    pds2.cond(cell_type="CD4 T cells", label="ctrl")
)

res_df = pds2.test_contrasts(interaction_contrast)

Show code cell output

Hide code cell output

Running Wald tests...

Log2 fold change & Wald test p-value, contrast vector: [ 0.  0.  0.  0.  0.  0.  0.  0.  1. -1.  0.  0.  0.  0.]
                baseMean  log2FoldChange      lfcSE      stat    pvalue  \
index                                                                     
AL627309.1      0.010607        1.740385  10.956334  0.158847  0.873789   
RP11-206L10.2   0.030188       -0.208432  10.939235 -0.019054  0.984798   
RP11-206L10.9   0.047391        0.917190  10.951037  0.083754  0.933252   
FAM87B          0.053617        0.882214  10.986500  0.080300  0.935999   
LINC00115       0.389827        0.455551   1.374820  0.331353  0.740378   
...                  ...             ...        ...       ...       ...   
C21orf58        0.121437        0.744942   4.379467  0.170099  0.864932   
PCNT            0.676612       -0.612964   1.089164 -0.562784  0.573582   
DIP2A           1.532521       -0.203075   0.684943 -0.296484  0.766861   
S100B           0.638731        0.746339   2.759453  0.270466  0.786802   
PRMT2          24.013222        0.705998   0.293153  2.408295  0.016027   

                   padj  
index                    
AL627309.1          NaN  
RP11-206L10.2       NaN  
RP11-206L10.9       NaN  
FAM87B              NaN  
LINC00115           NaN  
...                 ...  
C21orf58            NaN  
PCNT                NaN  
DIP2A          0.898603  
S100B               NaN  
PRMT2          0.079216  

[15701 rows x 6 columns]

... done in 0.66 seconds.

pds2.plot_volcano(res_df, log2fc_thresh=0)

NaNs encountered, dropping rows with NaNs

As we can see, there are some genes that stand out in their differences in gene expression changes. For example, the change in expression of RABGAP1L is higher in CD14+ monocytes compared to CD4 T cells, while ENO1 shows a lower change.

We can also plot the fold changes of the top differentially expressed genes. However, we only include results that pass a certain alpha threshold, for example 0.01.

pds2.plot_fold_change(res_df[res_df["adj_p_value"] < 0.01].copy(), n_top_vars=15)

Finally, we can also plot multiple comparisons. In this case, we are comparing the gene expression of CD14+ monocytes to all other cell types. Comparing cell types with each other should roughly reflect the marker genes used during annotation. This therefore mainly shows what is possible and should ideally be replaced by your groups of interest (e.g., sex, responder vs. non-responder status or age groups).

res_df = pds2.compare_groups(adata_pb,
    column="cell_type",
    baseline="CD14+ Monocytes",
    groups_to_compare=list(set(adata_pb.obs["cell_type"].unique()) - {"CD14+ Monocytes"})
)

Show code cell output

Hide code cell output

Fitting size factors...
... done in 0.03 seconds.

Using None as control genes, passed at DeseqDataSet initialization

Fitting dispersions...
... done in 2.75 seconds.

Fitting dispersion trend curve...
... done in 0.25 seconds.

Fitting MAP dispersions...
... done in 3.39 seconds.

Fitting LFCs...
... done in 2.42 seconds.

Calculating cook's distance...
... done in 0.05 seconds.

Replacing 6 outlier genes.

Fitting dispersions...
... done in 0.01 seconds.

Fitting MAP dispersions...
... done in 0.01 seconds.

Fitting LFCs...
... done in 0.01 seconds.

Running Wald tests...
... done in 0.57 seconds.

Running Wald tests...

Log2 fold change & Wald test p-value, contrast vector: [ 0. -1.  0.  0.  0.  0.  1.]
                baseMean  log2FoldChange     lfcSE      stat    pvalue  \
index                                                                    
AL627309.1      0.010607        2.474859  5.503200  0.449713  0.652918   
RP11-206L10.2   0.030188        2.282718  5.401389  0.422617  0.672575   
RP11-206L10.9   0.047391        2.677699  4.433057  0.604030  0.545824   
FAM87B          0.053617        2.282707  5.494428  0.415459  0.677806   
LINC00115       0.389827        1.199922  0.788086  1.522576  0.127865   
...                  ...             ...       ...       ...       ...   
C21orf58        0.121437        2.515586  1.834763  1.371069  0.170353   
PCNT            0.676612        2.175944  0.651016  3.342384  0.000831   
DIP2A           1.532521        1.044960  0.414023  2.523918  0.011606   
S100B           0.638731        3.882751  1.306370  2.972168  0.002957   
PRMT2          24.013222        0.373729  0.174184  2.145598  0.031905   

                   padj  
index                    
AL627309.1          NaN  
RP11-206L10.2       NaN  
RP11-206L10.9       NaN  
FAM87B              NaN  
LINC00115      0.190911  
...                 ...  
C21orf58            NaN  
PCNT           0.002286  
DIP2A          0.024516  
S100B          0.007275  
PRMT2          0.058900  

[15701 rows x 6 columns]

... done in 0.57 seconds.

Running Wald tests...

Log2 fold change & Wald test p-value, contrast vector: [ 0. -1.  1.  0.  0.  0.  0.]
                baseMean  log2FoldChange     lfcSE      stat    pvalue  \
index                                                                    
AL627309.1      0.010607        0.161952  5.470482  0.029605  0.976382   
RP11-206L10.2   0.030188       -0.246900  5.375897 -0.045927  0.963368   
RP11-206L10.9   0.047391       -0.389012  4.432479 -0.087764  0.930064   
FAM87B          0.053617       -0.768141  5.492885 -0.139843  0.888784   
LINC00115       0.389827        0.908217  0.612293  1.483305  0.137993   
...                  ...             ...       ...       ...       ...   
C21orf58        0.121437        1.199311  1.736388  0.690693  0.489759   
PCNT            0.676612        1.931798  0.502570  3.843835  0.000121   
DIP2A           1.532521        0.735263  0.325135  2.261409  0.023734   
S100B           0.638731        4.501091  1.235405  3.643412  0.000269   
PRMT2          24.013222        0.239543  0.153839  1.557103  0.119446   

                   padj  
index                    
AL627309.1          NaN  
RP11-206L10.2       NaN  
RP11-206L10.9       NaN  
FAM87B              NaN  
LINC00115      0.181667  
...                 ...  
C21orf58            NaN  
PCNT           0.000267  
DIP2A          0.036822  
S100B          0.000566  
PRMT2          0.159545  

[15701 rows x 6 columns]

... done in 0.58 seconds.

Running Wald tests...

Log2 fold change & Wald test p-value, contrast vector: [ 0. -1.  0.  0.  0.  0.  0.]
                baseMean  log2FoldChange     lfcSE      stat        pvalue  \
index                                                                        
AL627309.1      0.010607        1.472575  5.504369  0.267528  7.890624e-01   
RP11-206L10.2   0.030188        1.492242  5.391594  0.276772  7.819552e-01   
RP11-206L10.9   0.047391        1.491757  4.444019  0.335677  7.371141e-01   
FAM87B          0.053617        1.492273  5.484781  0.272075  7.855642e-01   
LINC00115       0.389827       -0.002214  0.804918 -0.002751  9.978054e-01   
...                  ...             ...       ...       ...           ...   
C21orf58        0.121437        1.376525  1.846643  0.745420  4.560177e-01   
PCNT            0.676612        0.990738  0.660558  1.499851  1.336531e-01   
DIP2A           1.532521       -0.665660  0.460780 -1.444638  1.485596e-01   
S100B           0.638731        1.090877  1.448696  0.753006  4.514461e-01   
PRMT2          24.013222       -1.228460  0.183585 -6.691519  2.208651e-11   

                       padj  
index                        
AL627309.1              NaN  
RP11-206L10.2           NaN  
RP11-206L10.9           NaN  
FAM87B                  NaN  
LINC00115      9.984005e-01  
...                     ...  
C21orf58                NaN  
PCNT           1.999094e-01  
DIP2A          2.184766e-01  
S100B          5.456290e-01  
PRMT2          1.420383e-10  

[15701 rows x 6 columns]

... done in 0.64 seconds.

Running Wald tests...

Log2 fold change & Wald test p-value, contrast vector: [ 0. -1.  0.  0.  1.  0.  0.]
                baseMean  log2FoldChange     lfcSE      stat    pvalue  \
index                                                                    
AL627309.1      0.010607        2.368673  5.495653  0.431009  0.666462   
RP11-206L10.2   0.030188        2.356868  5.385137  0.437662  0.661632   
RP11-206L10.9   0.047391        2.176679  4.447242  0.489445  0.624527   
FAM87B          0.053617        2.176522  5.486869  0.396678  0.691605   
LINC00115       0.389827        0.574520  0.876521  0.655454  0.512175   
...                  ...             ...       ...       ...       ...   
C21orf58        0.121437        2.037958  1.867656  1.091185  0.275192   
PCNT            0.676612        1.752980  0.702328  2.495958  0.012562   
DIP2A           1.532521       -0.969059  0.628585 -1.541651  0.123158   
S100B           0.638731        1.797302  1.482714  1.212170  0.225447   
PRMT2          24.013222       -0.333650  0.184535 -1.808052  0.070598   

                   padj  
index                    
AL627309.1          NaN  
RP11-206L10.2       NaN  
RP11-206L10.9       NaN  
FAM87B              NaN  
LINC00115      0.660674  
...                 ...  
C21orf58            NaN  
PCNT           0.046558  
DIP2A          0.255190  
S100B          0.385549  
PRMT2          0.172537  

[15701 rows x 6 columns]

... done in 0.62 seconds.

Running Wald tests...

Log2 fold change & Wald test p-value, contrast vector: [ 0. -1.  0.  1.  0.  0.  0.]
                baseMean  log2FoldChange     lfcSE      stat    pvalue  \
index                                                                    
AL627309.1      0.010607        3.396434  5.696649  0.596216  0.551031   
RP11-206L10.2   0.030188        3.074759  5.599790  0.549085  0.582947   
RP11-206L10.9   0.047391        3.140181  4.616491  0.680210  0.496372   
FAM87B          0.053617        3.074880  5.696030  0.539829  0.589315   
LINC00115       0.389827        1.314403  0.839342  1.565992  0.117351   
...                  ...             ...       ...       ...       ...   
C21orf58        0.121437        2.718599  1.954670  1.390823  0.164279   
PCNT            0.676612        2.018859  0.690688  2.922968  0.003467   
DIP2A           1.532521        1.157852  0.443806  2.608917  0.009083   
S100B           0.638731        3.976854  1.383560  2.874364  0.004048   
PRMT2          24.013222       -0.180823  0.202844 -0.891434  0.372696   

                   padj  
index                    
AL627309.1          NaN  
RP11-206L10.2       NaN  
RP11-206L10.9       NaN  
FAM87B              NaN  
LINC00115      0.179558  
...                 ...  
C21orf58       0.234673  
PCNT           0.009134  
DIP2A          0.021086  
S100B          0.010424  
PRMT2          0.456279  

[15701 rows x 6 columns]
Log2 fold change & Wald test p-value, contrast vector: [ 0. -1.  0.  0.  0.  1.  0.]
                baseMean  log2FoldChange     lfcSE      stat    pvalue  \
index                                                                    
AL627309.1      0.010607        2.151508  5.502361  0.391015  0.695786   
RP11-206L10.2   0.030188        1.959356  5.400545  0.362807  0.716749   
RP11-206L10.9   0.047391        1.959338  4.454816  0.439825  0.660064   
FAM87B          0.053617        2.295035  5.477388  0.419002  0.675215   
LINC00115       0.389827        1.067380  0.763056  1.398823  0.161866   
...                  ...             ...       ...       ...       ...   
C21orf58        0.121437        1.637420  1.891524  0.865662  0.386675   
PCNT            0.676612        2.224298  0.612325  3.632546  0.000281   
DIP2A           1.532521       -0.492096  0.514316 -0.956796  0.338670   
S100B           0.638731        1.730428  1.455929  1.188539  0.234621   
PRMT2          24.013222        0.148906  0.172725  0.862102  0.388632   

                   padj  
index                    
AL627309.1          NaN  
RP11-206L10.2       NaN  
RP11-206L10.9       NaN  
FAM87B              NaN  
LINC00115      0.343732  
...                 ...  
C21orf58            NaN  
PCNT           0.002758  
DIP2A          0.541026  
S100B          0.432238  
PRMT2          0.590751  

[15701 rows x 6 columns]

... done in 0.58 seconds.

pds2.plot_multicomparison_fc(res_df, n_top_vars=5, figsize=(12, 1.5))

The plot shows genes that are differentially expressed in all other cell types compared to CD14+ monocytes. For example, CLL2 and CCL7 are significantly less expressed in all other cell types compared to CD14+ monocytes.

18.8. Questions

18.8.1. Flipcards

What is differential gene expression and in which cases are we interested in testing for it?

DGE analysis identifies genes with significant expression differences between conditions, such as healthy versus diseased states, to elucidate underlying biological mechanisms.

What is the 'pseudoreplication' problem and how can it be circumvented?

In single-cell RNA sequencing, treating individual cells from the same subject as independent samples can inflate false discovery rates. This issue can be mitigated by aggregating data into 'pseudobulks' per subject or modeling subjects as random effects.

What is the 'multiple testing' problem and how can it be eluded?

Testing thousands of genes simultaneously increases the likelihood of false positives. To address this, corrections like the Benjamini-Hochberg procedure are applied to control the false discovery rate.

18.8.2. Multiple choice questions

What does a positive log2 fold change indicate in DESeq2 results?

The gene is expressed equally in both groups
The gene is more expressed in the baseline group
The gene is more expressed in the group of interest compared to baseline
The gene was filtered out during preprocessing

✔ Correct Answer: The gene is more expressed in the group of interest compared to baseline

What is the purpose of PCA on pseudobulk samples before DGE testing?

To normalize gene expression counts
To identify sources of variation that should be captured in the design matrix
To cluster cells into subtypes
To remove batch effects automatically

✔ Correct Answer: To identify sources of variation that should be captured in the design matrix

What does an interaction term like `cell_type:label` in a design matrix allow you to test?

Whether cell types differ in total read counts
Whether the effect of treatment differs across cell types
Whether replicates are correlated within patients
Whether batch effects are present in the data

✔ Correct Answer: Whether the effect of treatment differs across cell types

18.9. References

[deBKvB+17]

Mollie E. Brooks, Kasper Kristensen, Koen J. van Benthem, Arni Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin Mächler, and Benjamin M. Bolker. glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal, 9(2):378–400, 2017. URL: https://doi.org/10.32614/RJ-2017-066, doi:10.32614/RJ-2017-066.

[deDRM+21] (1,2)

Samarendra Das, Anil Rai, Michael L Merchant, Matthew C Cave, and Shesh N Rai. A comprehensive survey of statistical approaches for differential expression analysis in single-cell term`RNA` sequencing studies. Genes (Basel), 12(12):1947, December 2021.

[deFMY+15]

Greg Finak, Andrew McDavid, Masanao Yajima, Jingyuan Deng, Vivian Gersuk, Alex K. Shalek, Chloe K. Slichter, Hannah W. Miller, M. Juliana McElrath, Martin Prlic, Peter S. Linsley, and Raphael Gottardo. Mast: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell term`rna` sequencing data. Genome Biology, 16(1):278, Dec 2015. URL: https://doi.org/10.1186/s13059-015-0844-5, doi:10.1186/s13059-015-0844-5.

[deHHAN+21]

Yuhan Hao, Stephanie Hao, Erica Andersen-Nissen, William M. Mauck III, Shiwei Zheng, Andrew Butler, Maddie J. Lee, Aaron J. Wilk, Charlotte Darby, Michael Zagar, Paul Hoffman, Marlon Stoeckius, Efthymia Papalexi, Eleni P. Mimitou, Jaison Jain, Avi Srivastava, Tim Stuart, Lamar B. Fleming, Bertrand Yeung, Angela J. Rogers, Juliana M. McElrath, Catherine A. Blish, Raphael Gottardo, Peter Smibert, and Rahul Satija. Integrated analysis of multimodal single-cell data. Cell, 2021. URL: https://doi.org/10.1016/j.cell.2021.04.048, doi:10.1016/j.cell.2021.04.048.

[deHTTI17]

Stephanie C Hicks, F William Townes, Mingxiang Teng, and Rafael A Irizarry. Missing data and technical variability in single-cell term`RNA`-sequencing experiments. Biostatistics, 19(4):562–578, 11 2017. URL: https://doi.org/10.1093/biostatistics/kxx053, arXiv:https://academic.oup.com/biostatistics/article-pdf/19/4/562/26346801/kxx053.pdf, doi:10.1093/biostatistics/kxx053.

[deJSeyetermDNAsrollahME16]

Maria K Jaakkola, Fatemeh Seyeterm`DNA`srollah, Arfa Mehmood, and Laura L Elo. Comparison of methods to detect differentially expressed genes between single-cell populations. Briefings in Bioinformatics, 18(5):735–743, 07 2016. URL: https://doi.org/10.1093/bib/bbw057, arXiv:https://academic.oup.com/bib/article-pdf/18/5/735/25581122/bbw057.pdf, doi:10.1093/bib/bbw057.

[deJSE22] (1,2)

Sini Junttila, Johannes Smolander, and Laura L Elo. Benchmarking methods for detecting differential states between conditions from multi-subject single-cell term`rna`-seq data. bioRxiv, 2022. URL: https://www.biorxiv.org/content/early/2022/02/19/2022.02.16.480662, arXiv:https://www.biorxiv.org/content/early/2022/02/19/2022.02.16.480662.full.pdf, doi:10.1101/2022.02.16.480662.

[deKST+18]

Hyun Min Kang, Meena Subramaniam, Sasha Targ, Michelle Nguyen, Lenka Maliskova, Elizabeth McCarthy, Eunice Wan, Simon Wong, Lauren Byrnes, Cristina M Lanata, and others. Multiplexed droplet single-cell term`rna`-sequencing using natural genetic variation. Nature biotechnology, 36(1):89–94, 2018.

[deLZD+20]

Charity W. Law, Kathleen Zeglinski, Xueyi Dong, Monther Alhamdoosh, Gordon K. Smyth, and Matthew E. Ritchie. A guide to creating design matrices for gene expression experiments. F1000Research, 9:1444–1444, Dec 2020. 33604029[pmid]. URL: https://pubmed.ncbi.nlm.nih.gov/33604029.

[deLHA14]

Michael I. Love, Wolfgang Huber, and Simon Anders. Moderated estimation of fold change and dispersion for term`rna`-seq data with deseq2. Genome Biology, 15(12):550, Dec 2014. URL: https://doi.org/10.1186/s13059-014-0550-8, doi:10.1186/s13059-014-0550-8.

[deLT19]

Malte D Luecken and Fabian J Theis. Current best practices in single-cell term`rna`-seq analysis: a tutorial. Molecular Systems Biology, 15(6):e8746, 2019. URL: https://www.embopress.org/doi/abs/10.15252/msb.20188746, arXiv:https://www.embopress.org/doi/pdf/10.15252/msb.20188746, doi:https://doi.org/10.15252/msb.20188746.

[deLahnemannKosterS+20]

David Lähnemann, Johannes Köster, Ewa Szczurek, Davis J. McCarthy, Stephanie C. Hicks, Mark D. Robinson, Catalina A. Vallejos, Kieran R. Campbell, Niko Beerenwinkel, Ahmed Mahfouz, Luca Pinello, Pavel Skums, Alexandros Stamatakis, Camille Stephan-Otto Attolini, Samuel Aparicio, Jasmijn Baaijens, Marleen Balvert, Buys de Barbanson, Antonio Cappuccio, Giacomo Corleone, Bas E. Dutilh, Maria Florescu, Victor Guryev, Rens Holmer, Katharina Jahn, Thamar Jessurun Lobo, Emma M. Keizer, Indu Khatri, Szymon M. Kielbasa, Jan O. Korbel, Alexey M. Kozlov, Tzu-Hao Kuo, Boudewijn P.F. Lelieveldt, Ion I. Mandoiu, John C. Marioni, Tobias Marschall, Felix Mölder, Amir Niknejad, Lukasz Raczkowski, Marcel Reinders, Jeroen de Ridder, Antoine-Emmanuel Saliba, Antonios Somarakis, Oliver Stegle, Fabian J. Theis, Huan Yang, Alex Zelikovsky, Alice C. McHardy, Benjamin J. Raphael, Sohrab P. Shah, and Alexander Schönhuth. Eleven grand challenges in single-cell data science. Genome Biology, 21(1):31, Feb 2020. URL: https://doi.org/10.1186/s13059-020-1926-6, doi:10.1186/s13059-020-1926-6.

[deMS22]

Alan E. Murphy and Nathan G. Skene. A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis. Nat Commun, dec 2022. URL: https://doi.org/10.1038%2Fs41467-022-35519-4, doi:10.1038/s41467-022-35519-4.

[deRPW+15]

Matthew E. Ritchie, Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, and Gordon K. Smyth. Limma powers differential expression analyses for term`rna`-sequencing and microarray studies. Nucleic acids research, 43(7):e47–e47, Apr 2015. gkv007[PII]. URL: https://doi.org/10.1093/nar/gkv007, doi:10.1093/nar/gkv007.

[deRMS10] (1,2)

Mark D. Robinson, Davis J. McCarthy, and Gordon K. Smyth. Edger: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (Oxford, England), 26(1):139–140, Jan 2010. btp616[PII]. URL: https://doi.org/10.1093/bioinformatics/btp616, doi:10.1093/bioinformatics/btp616.

[deSR18]

Charlotte Soneson and Mark D. Robinson. Bias, robustness and scalability in single-cell differential expression analysis. Nature Methods, 15(4):255–261, Apr 2018. URL: https://doi.org/10.1038/nmeth.4612, doi:10.1038/nmeth.4612.

[deSGK+21] (1,2)

Jordan W. Squair, Matthieu Gautier, Claudia Kathe, Mark A. Anderson, Nicholas D. James, Thomas H. Hutson, Rémi Hudelle, Taha Qaiser, Kaya J. E. Matson, Quentin Barraud, Ariel J. Levine, Gioele La Manno, Michael A. Skinnider, and Grégoire Courtine. Confronting false discoveries in single-cell differential expression. Nature Communications, 12(1):5692, Sep 2021. URL: https://doi.org/10.1038/s41467-021-25960-2, doi:10.1038/s41467-021-25960-2.

[deTRCP21]

Andrew L Thurman, Jason A Ratcliff, Michael S Chimenti, and Alejandro A Pezzulo. Differential gene expression analysis for multi-subject single cell term`RNA` sequencing studies with aggregateBioVar. Bioinformatics, 37(19):3243–3251, May 2021.

[deVRS+17]

Catalina A. Vallejos, Davide Risso, Antonio Scialdone, Sandrine Dudoit, and John C. Marioni. Normalizing single-cell term`rna` sequencing data: challenges and opportunities. Nature Methods, 14(6):565–571, Jun 2017. URL: https://doi.org/10.1038/nmeth.4292, doi:10.1038/nmeth.4292.

[deWLNN19]

Tianyu Wang, Boyang Li, Craig E. Nelson, and Sheida Nabavi. Comparative analysis of differential gene expression analysis tools for single-cell term`rna` sequencing data. BMC Bioinformatics, 20(1):40, Jan 2019. URL: https://doi.org/10.1186/s12859-019-2599-6, doi:10.1186/s12859-019-2599-6.

[deZEL21] (1,2,3)

Kip D. Zimmerman, Mark A. Espeland, and Carl D. Langefeld. A practical solution to pseudoreplication bias in single-cell studies. Nature Communications, 12(1):738, Feb 2021. URL: https://doi.org/10.1038/s41467-021-21038-1, doi:10.1038/s41467-021-21038-1.

18.10. Contributors

We gratefully acknowledge the contributions of:

18.10.1. Authors

• Lukas Heumos

• Anastasia Litinetskaya

• Soroor Hediyeh-Zadeh

• Luis Heinzlmeier

18.10.2. Reviewers