Differential gene expression analysis

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.

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 Wang et al., 2019Das et al., 2021. 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., 2017Vallejos et al., 2017Luecken & Theis, 2019, methods designed for bulk RNA-seq data performed favorably compared to methods explicitly designed for scRNA-seq data Das et al., 2021Soneson & Robinson, 2018Jaakkola et al., 2016Squair 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) Squair et al., 2021Zimmerman et al., 2021Junttila et al., 2022. 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 & 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.

Environment setup¶

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.

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

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.

We have 8 control and 8 disease patients.

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

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.

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).

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.

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().

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.

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 Thurman et al., 2021Lähnemann et al., 2020. 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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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.

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 design matrix and then use contrasts to specify the conditions of interest.

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.

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).

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.

Questions¶

Flipcards¶

Multiple choice questions¶

Contributors¶

We gratefully acknowledge the contributions of:

Authors¶

• Lukas Heumos

• Anastasia Litinetskaya

• Soroor Hediyeh-Zadeh

• Luis Heinzlmeier

Reviewers¶