# Retrieving Co expression data from a Seurat object

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
Drd1a1   13.753313 100.00000         Drd1a        drd1     1.04216182


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(striatum@meta.data)*100


And here again I reach a completely different proportion...

• Have you tried/considered the use of FeatureScatter? You can set a threshold and select the cells on the diagonal, above the said threshold. To get the list of cells fitting your criteria, you can either use the ggplot object or use CellSelector.