I know that this topic has been discussed in other threads, but I never found a definitive answers (and maybe there isn't one) about assessing coexpression in single cell data.
I want to assess the expression of specific genes in a subpopulation of cells and report the percentage of co expression in each class. My focus is in the striatum and I have been using previously published data to do it so. This brain region is mainly constituted of two non overlapping population of cells, Drd1a positive and drd2 (or Adora2a) positive neurons and I want to investigate how my genes of interest are expressed in these subset of neurons. Based on previous threads about this subject, I was able to come to a solution, however, I am reaching different results depending on the strategy applied.
My input data is a counting matrix (GSE139265). After trimming, scaling and normalizing the data (no clustering), I use two different approaches to reach my goal:
Strategy 1, Dotplot
Idents(striatum, WhichCells(object = striatum, expression = Drd1a > 0, slot = 'data')) <- 'drd1_positive' Idents(striatum, WhichCells(object = striatum, expression = Adora2a > 0, slot = 'data')) <- 'a2a_positive'
And before assessing the overlap with my gene of interest, I checked what is the overlap between my newly created classes, and then my brain stated cracking.
percentage=DotPlot(striatum, features= c("Drd1a", "Adora2a"), cols=c("green", "red")) percentage$data
And this is what I get:
avg.exp pct.exp features.plot id avg.exp.scaled Drd1a 5.940562 57.19298 Drd1a a2a -0.09047198 Adora2a 8.210048 100.00000 Adora2a a2a 1.15470054 Drd1a1 13.753313 100.00000 Drd1a drd1 1.04216182 Adora2a1 0.000000 0.00000 Adora2a drd1 -0.57735027
So based on the results (zero overlap of Adora2a within Drd1a neurons but 57% overlap of Drd1a within Adora2a neurons), I am clearly I am doing something really wrong. But I can't really tell what it is.
My second approach is jus do the calculation of the percentage of overlap without inputting the new identities into the Seurat object, like this:
length(WhichCells(object = striatum, expression = Drd1a > 0 & Adora2a > 0))/nrow(email@example.com)*100
And here again I reach a completely different proportion...
I would really appreciate if somebody could clear my mind about this issue. Thanks a lot