单细胞转录组GSVA分析更新及可视化

很久以前，我们在写单细胞分析的时候，就写过GSVA分析了，GSVA分析的全名叫Gene Set Variation Analysis，能够分析基因集/通路的活性程度，具体的原理我们这里就不再赘述了（有了这个包，猪的GSEA和GSVA分析也不在话下（第一集））。这里我们对之前的GSVA分析进行一个更新，看一下不一样的差异分析思路和可视化！其他内容请参考：https://bioconductor.org/packages/devel/bioc/vignettes/GSVA/inst/doc/GSVA.html
我们这里以单细胞数据为例：

library(Seurat)
library(RColorBrewer)
library(GSVA)
library(GSEABase)
human_data <- readRDS("D:/KS项目/human_data.rds")

关于基因集，可以在GSEA官网上下载，也可以利用R包，例如：跟着Cell学单细胞转录组分析(十三):单细胞GSVA分析|这个包涵盖大多数物种。我们这里是直接从GSEA官网下载的，然后读入即可！

genesets <- getGmt('./c2.cp.v2023.2.Hs.symbols.gmt')

准备分析文件进行分析，也很简单，其实就是gene expression matrix。用counts或者normalized都行。在分析的时候，选择对method即可。

gene_matrix <- GetAssayData(human_data, layer = 'data')
GSVA <- gsva(as.matrix(gene_matrix), genesets, min.sz=10, max.sz=1000,verbose=TRUE)
head(rownames(GSVA))

分析完成之后，就可以进行差异分析了。之前我们讲的或者官网上写的都是利用limma包进行的分析。这里我们参考一篇NC文章，他将每个细胞的GSVA score构建为seurat obj，这样做差异分析和可视化都会很方便。

meta=human_data@meta.data[,c("orig.ident","celltype","group")]
row.names(meta)=colnames(GSVA)
GSVA_Seurat <- CreateSeuratObject(counts = GSVA, meta.data = meta, project = "GSVA_singleCell")
#celltype
Idents(GSVA_Seurat) <- "celltype"
GSVA_celltype=FindAllMarkers(GSVA_Seurat, only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.1)

#group
Macrophage <- subset(GSVA_Seurat, celltype=='Macrophage')
GSVA_DE = FindMarkers(Macrophage, ident.1 = "GM",ident.2 = "BM",
                      group.by = "group", min.pct = 0, logfc.threshold = 0)

火山图展示：如果需要展示通路也可以按照自己的需求展示！

seurat气泡图展示：

DotPlot(GSVA_Seurat, features = c("Wp circadian rhythm genes",
                                  "Reactome pi3k akt signaling in cancer",
                                  "Kegg jak stat signaling pathway",
                                  "Wp aerobic glycolysis",
                                  "Wp glycolysis in senescence",
                                  "Kegg medicus reference translation initiation"),group.by = "orig.ident")+coord_flip()+NoLegend()+
  theme(axis.title = element_blank())

然后就是热图展示·：

GSVA_exp <- AverageExpression(GSVA_Seurat,features = rownames(GSVA),group.by = 'celltype', assays = 'RNA', slot = "counts") 
GSVA_exp <- as.data.frame(GSVA_exp$RNA)

pathways <- c("Biocarta thelper pathway",
              "Biocarta tcytotoxic pathway",
              "Reactome digestion and absorption",
              "Reactome aspirin adme",
              "Biocarta il12 pathway",
              "Biocarta no2il12 pathway",
              "Reactome response to metal ions",
              "Kegg medicus reference translation initiation",
              "Reactome pd 1 signaling",
              "Kegg asthma",
              "Kegg arachidonic acid metabolism")

gsva_plot <- GSVA_exp[pathways,]

library(ComplexHeatmap)
pheatmap(gsva_plot,cluster_rows = F,cluster_cols = F,scale = 'none',
         colorRampPalette(c("#2166AC",'#478ABF','#90C0DC', "white",'#EF8C65','#CF4F45',"#B2182B"))(100),
         border=FALSE,cellwidth = 25, cellheight = 25,heatmap_legend_param = list(title="GSVA score"))

今天内容就分享到这里了，觉得有用的点个赞再走呗！