很久以前，我們在寫單細(xì)胞分析的時候，就寫過GSVA分析了，GSVA分析的全名叫Gene Set Variation Analysis，能夠分析基因集/通路的活性程度，具體的原理我們這里就不再贅述了（有了這個包，豬的GSEA和GSVA分析也不在話下（第一集））。這里我們對之前的GSVA分析進(jìn)行一個更新，看一下不一樣的差異分析思路和可視化！其他內(nèi)容請參考：https://bioconductor.org/packages/devel/bioc/vignettes/GSVA/inst/doc/GSVA.html
我們這里以單細(xì)胞數(shù)據(jù)為例：

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

關(guān)于基因集，可以在GSEA官網(wǎng)上下載，也可以利用R包，例如：跟著Cell學(xué)單細(xì)胞轉(zhuǎn)錄組分析(十三):單細(xì)胞GSVA分析|這個包涵蓋大多數(shù)物種。我們這里是直接從GSEA官網(wǎng)下載的，然后讀入即可！

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

準(zhǔn)備分析文件進(jìn)行分析，也很簡單，其實就是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))

分析完成之后，就可以進(jìn)行差異分析了。之前我們講的或者官網(wǎng)上寫的都是利用limma包進(jìn)行的分析。這里我們參考一篇NC文章，他將每個細(xì)胞的GSVA score構(gòu)建為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"))

今天內(nèi)容就分享到這里了，覺得有用的點(diǎn)個贊再走唄！