Seurat: FeaturePlot for average expression level of a set of genes

[YangShi-PKU]

HI
Thank you for developing such a powerful and user-friendly software. I am analyzing some drop-seq data by Seurat. In your vignette, you show how to visualize a feature (usually the expression level of a gene) on the tSNE plot. But as you know, some cell types cannot be well defined by only one marker gene; using a set of genes may be a better choice. So I wonder if I could plot the average expression level of a gene set as a feature on the tSNE plot. Is this possible using FeaturePlot? If not, is there some tricks to do it with the Seurat object? Any advice will be appreciated.
Thank you very much!

Yang

[leonfodoulian]

Hi Yang,

I am not sure if Seurat already has an implemented function for this, but here is my workaround for it:

# Select genes of interest (using sample() here for demonstration purposes)
gene.set <- sample(x = rownames(x = object@data), size = 100, replace = FALSE)

# Get mean expression of genes of interest per cell
mean.exp <- colMeans(x = object@data[gene.set, ], na.rm = TRUE)

# Add mean expression values in '[email protected]$gene.set.score'
if (all(names(x = mean.exp) == rownames(x = [email protected]))) {
  cat("Cell names order match in 'mean.exp' and '[email protected]':\n", 
      "adding gene set mean expression values in '[email protected]$gene.set.score'")
  [email protected]$gene.set.score <- mean.exp
}

# Plot mean expression using Seurat::FeaturePlot()
FeaturePlot(object = object, features.plot = "gene.set.score")

Steps are commented, but if you need further clarification as to what each code line is doing, please let me know.

Best,
Leon

[leonfodoulian]

Also, I don't know if averaging scaled values might be a better choice here compared to simply averaging raw or normalised expression values. You can simply exchange object@data[gene.set, ] with [email protected][gene.set, ] when computing colMeans().

Best,
Leon

[YangShi-PKU]

Thanks for your quick reply! Your advice perfectly solve my problem!

[YiweiNiu]

Hi,

I would like to provide another solution using AverageExpression

# Select genes of interest (using sample() here for demonstration purposes)
gene.set <- sample(x = rownames(x = object@data), size = 100, replace = FALSE)

# stash original cell identity
obj <- StashIdent(obj, save.name = 'cluster')

# add decoy cell identity
Idents(obj) <- '1'

# compute average gene expression
gene.set.exp <- geneAverageExpression(obj, features = gene.set)

# add it to meta.data
[email protected]$gene.set.exp <- gene.set.exp[,1]

# plot
FeaturePlot(object = obj, features.plot = "gene.set.exp")

# restore original cell identity
Idents(obj) <- 'cluster'

Through changing the cell identity, average gene expression can also be easily computed across other groups (e.g. control/case, tumor/normal).