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?