Discordance in gene signature behavior between bulk and single-cell RNASeq

Question

The objective of the following analysis is to identify an activation signature of a specific phenotype on bulk RNASeq and to apply it to single-cell RNA-Seq, in order to identify the population of activated cells in scRNA-Seq samples.

In the first step I used DESeq2 to identify genes that increase during different timepoints (I have 4 timepoints):

• Bulk samples are normalized using VST (Variance Stabilizing Transformation - DESeq2)
• For each timepoint, differentially expressed genes are computed with respect to the previous timepoint (i.e. 2 vs 1, 3 vs 2, 4 vs 3)
• Differentially expressed genes that increase with timepoint (log-FC > 1, padj < .05) are selected (n=24)

As a control, I plotted the activation score (computed with scanpy) based on the 24-gene signature applied to the bulk samples from which it has been extracted. As you can see, the score always increases from the first timepoint (blue) to the last (red).

However, when I use the same scoring procedure on a single-cell RNASeq sample containing the same populations, I think something may be wrong.

In the following the violin plot of the extracted activation score on an unsupervised clustering on single-cell data. Excluding cluster 4, which we identified as contamination, while cluster 1 and cluster 0 are supposed to have the correct degree of activation, clusters 2 and 3 have a lower degree of activation with respect to cluster 0 (which we identified as the first timepoint - lowest activation on bulk).

What may be the causes of this behavior? Right now the possible explanation I came out with is that:

• Clusters 2 and 3 have global shifts on gene expression that make TPM normalization over-normalize activation genes

Are there other possible explanations for this behavior, and how to check for them?

Answer 1

One explanation could be that your mapping of clusters to timepoints is not accurate. There are other methods you could look at for doing this, for example scMap, scPred, or Seurat v3 (disclosure: I am one of the Seurat developers).

Another possibility is that there is some problem in the calculation of the activation signature. You could instead look at the expression of the 24 genes themselves in each cluster (for example by drawing a heatmap).

The difference between bulk and single-cell here could also be biological. This could happen if the increase in expression of those genes in bulk was due to a change in cell type frequency over time rather than increased expression.