Featureplot by co-expression of some genes

Question

I want to identify MDSCs cells which don't not express HLA-DRs (HLA-DR negative) but express CD14, CD33, and ITGAM

For this, I firstly plotted a feature plot to locate HLA-DR by

cd.features <- list(grep(pattern = '^HLA-DR', x = rownames(x = pbmc_small), value = TRUE))

pbmc_small <- AddModuleScore(object = pbmc_small, features = cd.features, ctrl = 5, name = 'CDFeatures')


FeaturePlot(object = pbmc_small, features = 'CDFeatures1')& 
+     scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

Which returned

For CD14, CD33, and ITGAM I done

> cd.features=list(c("CD14", "CD33",  "ITGAM"))
> pbmc_small <- AddModuleScore(object = pbmc_small, features = cd.features, ctrl = 5, name = 'CDFeatures')
> FeaturePlot(object = pbmc_small, features = 'CDFeatures1')& 
+     scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

> 

Which returned

By looking at these two plots I really can not say which part of plot is HLA-DR negative but expressing CD14, CD33, and ITGAM

Do you know a way to do this?

Plotting options don't help because looking at the below plots you can not say which cells are HLA negative and CD14 positive

Answer 1

Very old question, but recently modified, which means that someone might still be looking into this.

I usually prefer to use Ucell as the signatures are calculated within cells, which makes comparisons a bit easier. In addition is allows the formulation of signatures with both positive and negative genes.

hladr_genes <- grep(pattern = '^HLA-DR', x = rownames(x = pbmc_small), value = TRUE)
signatures_list <- list(MDSCs = c("CD14+","CD33+","ITGAM+",paste0(hladr_genes,"-")))

> signatures_list
$MDSCs
[1] "CD14+"     "CD33+"     "ITGAM+"    "HLA-DRA-"  "HLA-DRB5-" "HLA-DRB1-"

# calculate signature
pbmc_small <- AddModuleScore_UCell(pbmc_small,signatures_list,ncores = 8,name = "")
# perform smoothing
pbmc_small <- SmoothKNN(pbmc_small, signature.names = names(signatures_list), reduction="pca")

Answer 2

You need to make a transformed variable to the gene sets of interest. I would rescale the data on a range of 0-1. Then you can make a log2 ratio of the sum(set1)/sum(set2). You could explore other transformations, but this will give contrasting values for high set1 and high set2, but neutral if the expression is similar between the sets.