连续两次求贤令:曾经我给你带来了十万用户,但现在祝你倒闭,以及 生信技能树知识整理实习生招募,让我走大运结识了几位优秀小伙伴!大家开始根据我的ngs组学视频进行一系列公共数据集分析实战,其中几个小伙伴让我非常惊喜,不需要怎么沟通和指导,就默默的完成了一个实战!
他前面的分享是:
下面是他对我们b站转录组视频课程的详细笔记
承接上节:RNA-seq入门实战(四):差异分析前的准备——数据检查,以及 RNA-seq入门实战(五):差异分析——DESeq2 edgeR limma的使用与比较
本节概览:
1.GSEA简单介绍
2.创建GSEA分析所需的geneList,包含log2FoldChange和ENTREZID信息
3.利用clusterProfiler进行GSEA富集GO与KEGG通路
4.GSEA富集结果可视化:GSEA结果图、 gsearank plot 、ridgeplot
以下对GSEA涉及的一些重要概念进行了简单介绍,详细介绍见: 一文掌握GSEA,超详细教程 - 云+社区 - 腾讯云 (tencent.com)史上最全GSEA可视化教程,今天让你彻底搞懂GSEA! - 知乎 (zhihu.com)
分子特征数据库。一般进行GSEA或GSVA使用的就是该数据库中的基因集,我们也可以自定义基因集。MSigDB所包含的基因集如下所示:
在了解了GSEA基本概念后就可以正式开始实操了,首先需要将基因按照在两类样本中的差异表达程度排序。 下面我们构建包含了geneList,里面含有从大到小排序的log2FoldChange和对应的ENTREZID信息:
rm(list = ls())
options(stringsAsFactors = F)
# library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(clusterProfiler)
library(enrichplot)
library(tidyverse)
library(ggstatsplot)
setwd("C:/Users/Lenovo/Desktop/test")
source('FUNCTIONS.R')
load(list.files(path = "./3.DEG",pattern = 'DEG_results.Rdata',full.names = T))
dir.create("5.GSEA_kegg_go")
setwd("5.GSEA_kegg_go")
## 物种设置
organism = 'mmu' # 人类'hsa' 小鼠'mmu'
OrgDb = 'org.Mm.eg.db'#人类"org.Hs.eg.db" 小鼠"org.Mm.eg.db"
#### 按照需要可选择不同的DEG方法数据集 ####
need_DEG <- DEG_DESeq2
need_DEG <- need_DEG[,c(2,5)] #选择log2FoldChange和pvalue(凑成数据框)
colnames(need_DEG) <- c('log2FoldChange','pvalue')
need_DEG$SYMBOL <- rownames(need_DEG)
##### 创建gsea分析的geneList(包含从大到小排列的log2FoldChange和ENTREZID信息)####
#转化id
df <- bitr(rownames(need_DEG),
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = OrgDb) #人数据库org.Hs.eg.db 小鼠org.Mm.eg.db
need_DEG <- merge(need_DEG, df, by='SYMBOL') #按照SYMBOL合并注释信息
geneList <- need_DEG$log2FoldChange
names(geneList) <- need_DEG$ENTREZID
geneList <- sort(geneList, decreasing = T) #从大到小排序
clusterProfiler包内的gseGO()和gseKEGG()函数可以很方便地对GO与KEGG通路进行GSEA, 再使用DOSE::setReadable转化id 。
##### gsea富集 ####
KEGG_kk_entrez <- gseKEGG(geneList = geneList,
organism = organism, #人hsa 鼠mmu
pvalueCutoff = 0.25) #实际为padj阈值可调整
KEGG_kk <- DOSE::setReadable(KEGG_kk_entrez,
OrgDb=OrgDb,
keyType='ENTREZID')#转化id
GO_kk_entrez <- gseGO(geneList = geneList,
ont = "ALL", # "BP"、"MF"和"CC"或"ALL"
OrgDb = OrgDb,#人类org.Hs.eg.db 鼠org.Mm.eg.db
keyType = "ENTREZID",
pvalueCutoff = 0.25) #实际为padj阈值可调整
GO_kk <- DOSE::setReadable(GO_kk_entrez,
OrgDb=OrgDb,
keyType='ENTREZID')#转化id
save(KEGG_kk_entrez, GO_kk_entrez, file = "GSEA_result.RData")
GSEA的可视化主要是GSEA结果图、 gsearank plot和ridgeplot山脊图。 同样也可以进行其他可视化如barplot、dotplot、cnetplot等等,详见RNA-seq入门的简单实战(六):GO、KEGG富集分析与超全可视化攻略或者参阅说明书Chapter 15 Visualization of functional enrichment result | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top),这里就不再进行展示啦
下面选取KEGG通路的富集结果进行gseaplot绘图示范
##选取富集结果
kk_gse <- KEGG_kk
kk_gse_entrez <- KEGG_kk_entrez
###条件筛选
#一般认为|NES|>1,NOM pvalue<0.05,FDR(padj)<0.25的通路是显著富集的
kk_gse_cut <- kk_gse[kk_gse$pvalue<0.05 & kk_gse$p.adjust<0.25 & abs(kk_gse$NES)>1]
kk_gse_cut_down <- kk_gse_cut[kk_gse_cut$NES < 0,]
kk_gse_cut_up <- kk_gse_cut[kk_gse_cut$NES > 0,]
#选择展现NES前几个通路
down_gsea <- kk_gse_cut_down[tail(order(kk_gse_cut_down$NES,decreasing = T),10),]
up_gsea <- kk_gse_cut_up[head(order(kk_gse_cut_up$NES,decreasing = T),10),]
diff_gsea <- kk_gse_cut[head(order(abs(kk_gse_cut$NES),decreasing = T),10),]
#### 经典的GSEA图
up_gsea$Description
i=2
gseap1 <- gseaplot2(kk_gse,
up_gsea$ID[i],#富集的ID编号
title = up_gsea$Description[i],#标题
color = "red", #GSEA线条颜色
base_size = 20,#基础字体大小
rel_heights = c(1.5, 0.5, 1),#副图的相对高度
subplots = 1:3, #要显示哪些副图 如subplots=c(1,3) #只要第一和第三个图
ES_geom = "line", #enrichment score用线还是用点"dot"
pvalue_table = T) #显示pvalue等信息
ggsave(gseap1, filename = 'GSEA_up_1.pdf', width =10, height =8)
#### 合并 GSEA通路
gseap2 <- gseaplot2(kk_gse,
up_gsea$ID,#富集的ID编号
title = "UP_GSEA_all",#标题
color = "red",#GSEA线条颜色
base_size = 20,#基础字体大小
rel_heights = c(1.5, 0.5, 1),#副图的相对高度
subplots = 1:3, #要显示哪些副图 如subplots=c(1,3) #只要第一和第三个图
ES_geom = "line",#enrichment score用线还是用点"dot"
pvalue_table = T) #显示pvalue等信息
ggsave(gseap2, filename = "GSEA_up_all.pdf",width =12,height =12)
下面解释一下GSEA图的含义:
gsearank()展示特定基因集的排序,横坐标为基因排序,纵坐标为ES值,利用cowplot和ggplot2包可以批量出图。
## gsearank plot 绘制出属于特定基因集的基因排序列表
##绘制up_gsea前3个富集通路
library(cowplot)
library(ggplot2)
pp <- lapply(1:3, function(i) {
anno <- up_gsea[i, c("NES", "pvalue", "p.adjust")]
lab <- paste0(names(anno), "=", round(anno, 3), collapse="\n")
gsearank(kk_gse,
up_gsea$ID[1],
title = up_gsea$Description[i]) +
xlab(NULL) +
# ylab(NULL) +
annotate("text", 10000,
up_gsea[i, "enrichmentScore"] * .75,
label = lab,
hjust=0, vjust=0)
})
rankp <- plot_grid(plotlist=pp, ncol=1)
ggsave(rankp, filename = "gsearank_up.pdf",width=8,height=10)
展示富集通路的核心富集基因的表达分布,x轴为富集通路的核心富集基因表达变化的log2倍,值为正值表示表达上调,值为负值表示表达下调。
## ridgeplot
ridgep <- ridgeplot(kk_gse_entrez,
showCategory = 15,
fill = "p.adjust",
core_enrichment = TRUE,
label_format = 30, #设置轴标签文字的每行字符数长度,过长则会自动换行。
orderBy = "NES",
decreasing = F)
ggsave(ridgep,filename = 'ridgeplot.pdf',width =10,height =10)
(之前运行报错解决方法见ridgeplot报错:Error in ans[ypos] <- rep(yes, length.out = len)[ypos] : replacement has ... )
dotplot cnetplot emapplot treeplot heatplot upsetplot 详见RNA-seq入门的简单实战(六):GO、KEGG富集分析与超全可视化攻略
GSEA分析和可视化到这就结束啦,下一节介绍GSVA的使用
参考资料
📖 Introduction | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top)
一文掌握GSEA,超详细教程 - 云+社区 - 腾讯云 (tencent.com)
史上最全GSEA可视化教程,今天让你彻底搞懂GSEA! - 知乎 (zhihu.com)
GitHub - jmzeng1314/GEO
【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili
【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili