单细胞测序—标准分析流程(4)—GSEA与GSVA

原创

sheldor没耳朵

发布于 2024-09-05 18:04:04

1710

发布于 2024-09-05 18:04:04

文章被收录于专栏：单细胞测序

单细胞测序—标准分析流程(4)—GSEA与GSVA

这部分代码是我综合了好几篇帖子手打的代码

主要参考的是单细胞绘图之GSEA & GSVA

再调用GSVA函数出问题时主要参考：

https://github.com/rcastelo/GSVA/issues/172

GSEA和GSVA，再也不用去下载gmt文件咯

1 GSEA

接着单细胞测序—标准流程代码（3）—marker 基因富集分析_差异基因,继续分析

gesa_gsva_bymyself.R

rm(list=ls())
options(stringsAsFactors = F) 
source('scRNA_scripts/lib.R')
library(clusterProfiler)
library(enrichplot)

getwd()
d='gesa_gsva_bymyself/'
dir.create(d)
setwd(d) 
# 数据导入
sce = readRDS('../2-harmony/sce.all_int.rds')
load('../phe.Rdata')
sce$celltype=phe$celltype
sce$group=sce$stim
table(sce$celltype,sce$group )

######## GSEA #####
table(sce@meta.data$celltype)
cells_sub <- subset(sce@meta.data,celltype %in% c("CD4 Naive T")) 
scRNA_sub <- subset(sce,cells = row.names(cells_sub))
#logfc.threshold设置很小，因为GSEA尽量要所有的基因
sub.markers <- FindMarkers(scRNA_sub,group.by = 'group',
                           ident.1 = 'CTRL',ident.2='STIM',logfc.threshold=0.01)

sub.markers.sig <- subset(sub.markers,p_val_adj<0.05 & avg_log2FC > 0.5)

# GSEA

mydata <- data.frame(Gene = rownames(sub.markers),
                     logFC = sub.markers$avg_log2FC) %>% arrange(desc(logFC))
genelist <- bitr(mydata$Gene,fromType = "SYMBOL",toType = "ENTREZID",OrgDb = "org.Hs.eg.db") %>% 
  inner_join(mydata,by=c("SYMBOL"="Gene")) %>%arrange(desc(logFC))

gsea_input <- genelist$logFC 
names(gsea_input) <- genelist$ENTREZID
        
head(gsea_input)


geneset <- read.gmt("../c5.go.bp.v2024.1.Hs.entrez.gmt")
gg <- GSEA(gsea_input,TERM2GENE = geneset,verbose = F,
           pvalueCutoff = 0.1,pAdjustMethod = "BH")
sortgg <- gg[order(gg@result$NES,decreasing = T),]
sortgg <- sortgg[sortgg$p.adjust < 0.05,]


#这里取前3个通路进行绘图
head(sortgg)
geneset_plot <- c("GOBP_UROGENITAL_SYSTEM_DEVELOPMENT",
                  "GOBP_CELL_CELL_JUNCTION_ASSEMBLY",
                  "GOBP_CELL_CELL_JUNCTION_ORGANIZATION")
mycol <- pal_nejm()(8)
gseaplot2(gg,
          geneSetID =geneset_plot,
          color = mycol[c(1:3)],
          title='Specific GO_BP',
          rel_heights = c(1.3,0.3,0.6),
          pvalue_table = F)
ggsave("gsea.pdf",width = 10,height = 10)
save(gg,sortgg,file = "gsea.Rdata")

代码解释

同样的在数据导入部分，我认为可以直接直接load，降维聚类分群注释后的seurat对象（sce.all.int），这样的效果和上述三行的效果一致。

接下来，如本次以对照组(CTRL)与刺激组(STIM)中的CD4 Naive T细胞亚群为例，进行GSEA分析，后续重要关注哪两个分组或哪些细胞亚群，修改对应的代码即可。

提取CD4 Naive T细胞亚群的差异表达基因：

使用subset提取细胞类型为CD4 Naive T的亚群。
使用FindMarkers函数找到对照组(CTRL)与刺激组(STIM)之间的差异表达基因。为了保证GSEA能使用所有基因，logfc.threshold设置为0.01（阈值非常低）。
筛选出显著的差异基因：调整后的p值小于0.05，且log2倍数变化大于0.5。

准备GSEA输入数据：

将差异基因的基因名与logFC（log2倍数变化）值整理为一个数据框mydata。
使用bitr函数将基因符号（SYMBOL）转换为基因的Entrez ID，并与logFC数据合并。
创建GSEA输入向量gsea_input，其中基因的Entrez ID为向量的名字，logFC为向量的值。

执行GSEA分析：

使用read.gmt加载GO生物过程（BP）基因集。

这里我是按照视频里去官网https://www.gsea-msigdb.org/gsea/msigdb/human/collections.jsp#C5下载了对应的GO BP gmt文件。

这段代码是有优化的空间的，应该可以用clusterProfiler这个包直接分析，这里还没有研究，之后有时间再优化吧

使用GSEA函数进行富集分析，设定p值阈值为0.1，p值校正方法为Benjamini-Hochberg。

对GSEA结果按照NES（标准化富集得分）进行排序，并筛选调整后的p值小于0.05的结果。

绘制GSEA结果：

选择前三个通路进行绘图
使用gseaplot2绘制GSEA结果图，并保存为PDF文件。
使用save函数保存GSEA分析的结果到文件gsea.Rdata。

gsea.pdf

2 GSVA

脚本接上

######## GSVA #####

expr <- as.data.frame(scRNA_sub@assays$RNA@data)
expr$ID <- rownames(expr)
s2e <- bitr(expr$ID,
                 fromType = "SYMBOL",
                 toType = "ENTREZID",
                 OrgDb = "org.Hs.eg.db")
expr <- inner_join(expr,s2e,by=c("ID"="SYMBOL"))
rownames(expr) <- expr$ENTREZID

expr[1:3,1:3]
expr[1:5,c((ncol(expr)-1),ncol(expr))]
#删除最后两列，即SYMBOL，ENTREZID两列
#这对表达矩阵来说，这两列并不需要
expr <- expr[,-c((ncol(expr)-1),ncol(expr))]
geneset <- read.gmt("../c2.cp.kegg_medicus.v2024.1.Hs.entrez.gmt")
geneset_list <- split(geneset$gene,geneset$term)


library(GSVA)
expr <- as.matrix(expr)
#gsva1 <- gsva(gsvaParam(xpr,geneset_list,kcdf="Gaussian",method = "gsva",parallel.sz = 12))
gsva1 <- gsva(gsvaParam(expr, geneset_list))

meta <- as.data.frame(scRNA_sub@meta.data[,c("orig.ident", "group","celltype")])
mydata <- t(as.data.frame(gsva1)) %>% as.data.frame()
mydata <- merge(meta,mydata,by.x = 0,by.y = 0)
rownames(mydata) <- mydata$Row.names
mydata <- mydata[,c(-1,-2,-4)]

#差异分析
library(limma)
group_list <- factor(mydata$group)
design <- model.matrix(~0+group_list)
colnames(design) <- levels(group_list)
rownames(design) <- rownames(mydata)
contrast.matrix <- makeContrasts('STIM-CTRL',levels = design)
fit <- lmFit(t(mydata[,-1]),design = design)
fit2 <- contrasts.fit(fit,contrast.matrix)
fit2 <- eBayes(fit2)
alldiff = topTable(fit2,coef =1,n = Inf)
sigdiff <- subset(alldiff,abs(logFC)>0.5 & adj.P.Val < 0.05)
save(alldiff,sigdiff,file = "GSVA.Rdata")

#Heatmap
alldiff <- alldiff[order(abs(alldiff$logFC),decreasing = T),]
plotdata <- t(mydata[,rownames(alldiff[1:20,])])
library(ComplexHeatmap)
library(circlize)

Group <- c('#D2691E','#DDA0DD')
names(Group) <- c('STIM','CTRL')

Top <- HeatmapAnnotation(Group=factor(mydata$group),
                         annotation_legend_param = list(labels_gp = gpar(fontsize = 10),
                                                        title_gp = gpar(fontsize=10,fontface = "bold"),ncol=1),
                         border = T,
                         col = list(Group=Group),
                         show_annotation_name = TRUE,
                         annotation_name_side = "left",
                         annotation_name_gp = gpar(fontsize=10))
pdf(file = 'GSVA_heatmap.pdf',width=8,height = 6) 
Heatmap(plotdata,
        name="Scores",
        top_annotation = Top,
        cluster_rows = FALSE,
        col = colorRamp2(c(-1,0,1),c('#008B8B','#F5F5F5','#DC143C')),
        color_space = "RGB",
        cluster_columns = FALSE,border =T,
        row_names_side ='left',
        column_order=NULL,
        show_column_names =FALSE,
        row_names_gp = gpar(fontsize = 8),
        column_split = c(rep(1,table(mydata$group)[1]),rep(2,table(mydata$group)[2])),
        gap = unit(4,"mm"),
        use_raster=FALSE,
        column_title = NULL,
        column_title_gp = gpar(fontsize = 10),
        show_heatmap_legend = TRUE,
        heatmap_legend_param=list(labels_gp = gpar(fontsize = 10),
                                  title_gp = gpar(fontsize=10,fontface = "bold")))
dev.off()

代码解释

这段代码主要进行GSVA（Gene Set Variation Analysis）分析，识别在不同条件下基因集活性的变化，并生成热图展示基因集的差异表达情况。以下是代码的逐步解释：

GSVA分析数据准备：

expr$ID <- rownames(expr)：将基因名作为expr数据框的一列。
s2e <- bitr(expr$ID, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db")：使用bitr函数将基因符号（SYMBOL）转换为Entrez ID，以便后续分析使用。
expr <- inner_join(expr, s2e, by = c("ID" = "SYMBOL"))：将Entrez ID与表达数据结合。
rownames(expr) <- expr$ENTREZID：将基因的Entrez ID作为表达矩阵的行名。
expr <- expr,-c((ncol(expr)-1), ncol(expr))：删除最后两列（SYMBOL和ENTREZID），因为表达矩阵中不需要这两列。

读取基因集并进行GSVA分析：

geneset <- read.gmt("../c2.cp.kegg_medicus.v2024.1.Hs.entrez.gmt")：读取KEGG基因集文件。
geneset_list <- split(geneset$gene, geneset$term)：将基因集按通路名称（term）进行分类，生成每个通路包含的基因列表。
expr <- as.matrix(expr)：将表达数据转换为矩阵形式。
gsva1 <- gsva(gsvaParam(expr, geneset_list))：使用gsva函数进行基因集变异分析，输入为表达矩阵和基因集列表，输出每个基因集在不同细胞中的活性评分。

这一步是参考了github的报错，gsva更新了调用方式，导致之前的许多代码报错。要加上gsvaParam参数。

合并元数据和GSVA结果：

meta <- as.data.frame(scRNA_sub@meta.data, c("orig.ident", "group", "celltype"))：提取子集的元数据。
mydata <- t(as.data.frame(gsva1)) %>% as.data.frame()：将GSVA结果转置并转换为数据框形式。
mydata <- merge(meta, mydata, by.x = 0, by.y = 0)：将元数据与GSVA结果合并。
rownames(mydata) <- mydata$Row.names：设置合并后的数据框行名。
mydata <- mydata,c(-1,-2,-4)：删除不需要的列，保留关键的表达数据和分组信息。

差异分析：

group_list <- factor(mydata$group)：将分组信息转换为因子。
design <- model.matrix(~0+group_list)：构建设计矩阵，描述分组情况。
colnames(design) <- levels(group_list)：将设计矩阵列名设为分组名。
contrast.matrix <- makeContrasts('STIM-CTRL', levels = design)：创建对比矩阵，用于比较STIM组和CTRL组之间的差异。
fit <- lmFit(t(mydata,-1), design)：使用线性模型拟合GSVA结果与分组信息。
fit2 <- contrasts.fit(fit, contrast.matrix)：应用对比矩阵进行组间差异比较。
fit2 <- eBayes(fit2)：使用eBayes方法对结果进行贝叶斯调整。
alldiff = topTable(fit2, coef = 1, n = Inf)：提取差异分析的结果，按logFC排序。

sigdiff <- subset(alldiff, abs(logFC) > 0.5 & adj.P.Val < 0.05)：筛选出具有显著差异的基因集（logFC绝对值大于0.5，且调整后的p值小于0.05）。
save(alldiff, sigdiff, file = "GSVA.Rdata")：保存差异分析结果。

热图绘制：

alldiff <- alldifforder(abs(alldiff$logFC), decreasing = T),：按logFC绝对值对差异分析结果排序。
plotdata <- t(mydata[, rownames(alldiff1:20,)])：提取前20个差异最大的基因集，生成用于绘图的数据。
定义Group颜色，用于表示分组信息。
使用HeatmapAnnotation创建热图的分组注释。
使用Heatmap函数绘制热图，显示基因集的活性评分，按分组对数据进行分列，颜色从青色（低表达）到红色（高表达）渐变。
最终将热图保存为GSVA_heatmap.pdf文件。

总结来说，这段代码进行了基因集变异分析（GSVA），找出了在不同组别（STIM vs. CTRL）之间显著差异的基因集，并将前20个基因集的活性评分通过热图可视化。

GSVA_heatmap.pdf