5_Gene set enrichment

Gene set enrichment is a widely used analysis to get a sense of what cells are doing. There are many packages out there, so we listed the packages/tutorials we have tried and example code of what we frequently use. Most, if not all, terms used in gene set enrichment come from MSigDb, so it's best to browse this page prior to analysis.

Tutorials

clusterProfiler

piano

scPubR

fgsea

Visualizing individual MSigDB gene sets

library(ArchR)
library(msigbr)
library(Seurat)
library(scCustomize)

rna_cols<-paletteContinuous(n=8)[c(1:3, 6:8)]

#get gene ontology biological process collection
msig<-msigdbr(subcategory = "GO:BP")

#find term for extracellular matrix
msig$gs_name[grep("EXTRACELLULAR_MATRIX",msig$gs_name)] %>% unique()

#get ecm gene list
ecm<-msig[msig$gs_name == "GOBP_REGULATION_OF_EXTRACELLULAR_MATRIX_ORGANIZATION",]$gene_symbol

seurat_object<-AddModuleScore(seurat_object, features = list(ecm), name = "ecm")
FeaturePlot_scCustom(seurat_object, features = "ecm1")+scale_color_gradientn(colors = rna_cols)

Using GSEA to compare terms across seurat clusters

library(Seurat)
library(ArchR)
library(clusterProfiler)
library(msigdbr)
library(ggplot2)
library(dplyr)


# Initialize a results table
results <- data.frame(cluster_id = character(),
                      geneID = character(),
                      p_value = numeric(),
                      ID = character(),
                      Description = character(),
                      BgRatio = character(),
                      p.adjust = numeric(),
                      Count = numeric(),
                      qvalue = numeric(),
                      stringsAsFactors = FALSE)

DefaultAssay(seurat_object)<-"RNA"
markers<-FindAllMarkers(seurat_object)

# Loop through each cluster and gene set, and perform enrichment analysis

enrich_markers<-function(cluster_id){
  cluster_id <- cluster_id
  marker_sub<-markers[markers$cluster == cluster_id,]
  genes<-marker_sub$gene
  
  # you will need to change species, category, and subcategory to which MSigDB gene set you want
  
  #this selected homo spaiens, gene ontology: biological process
  m_t2g <- msigdbr(species = "Homo sapiens", category = "C5", subcategory = "GO:BP") %>% 
  dplyr::select(gs_name, gene_symbol) 

  result <- enricher(gene =genes,
                     TERM2GENE = m_t2g,
                     pvalueCutoff = 0.05,
                     universe = rownames(seurat_object),
                     qvalueCutoff = 0.05,
                     pAdjustMethod = "BH")
  if (length(result@result) > 0) {
    results <- result@result
    results$cluster_id<-cluster_id
    return(results)
  }else{
    return(NULL)
  }
  
}

# Use lapply to perform enrichment analysis for all cluster-gene set combinations
results_list <- lapply(levels(factor(seurat_object$seurat_clusters)),function(i){
  enrich_markers(i)})

# Flatten the results list and remove NULL elements
results <- do.call(rbind, results_list)

#this gets the top ~100 terms per cluster based on pvalue & count, feel free to change the number as you see fit
top_terms<-lapply(levels(factor(results$cluster_id)), function(x){
  
  top_gse_sub<-results  %>% dplyr::filter(cluster_id == x & pvalue < 0.05) 
  top_gse_filt<-top_gse_sub %>% slice_min(pvalue, n =100) %>% slice_max(Count, n =100)
  top_gse_filt$Description
  
})

#get unique top terms
top_terms<-top_terms %>% unlist() %>% unique()

#get values for top terms from each cluster
results_sub<-results[results$Description %in% top_terms,]

#make sure to display significant values
results_sub<-results_sub[results_sub$pvalue< 0.05,]

#organize data frame by cluster
results_sub<-results_sub %>% group_by(cluster_id)

#order levels to make plotting nice :)
results_sub$Description<-factor(results_sub$Description, levels =rev(unique(results_sub$Description)))
results_sub$cluster_id<-factor(results_sub$cluster_id, levels =unique(results_sub$cluster_id))

ggplot(results_sub, aes(x = Description, y = cluster_id, size = Count, fill = pvalue))+
geom_point(shape = 21, color="black")+
theme_minimal()+
RotatedAxis()+
coord_flip()+
scale_fill_gradientn(colors = paletteContinuous(set = "whiteBlue", reverse = T))