clusterProfiler|GSEA富集分析及可视化
GSEA（Gene Set EnrichmentAnalysis），即基因集富集分析，无需设定阈值来区分上调下调基因，使用所有的基因进行分析。
GO 和 KEGG 可参考：R|clusterProfiler-富集分析
一 准备数据
1.1，加载R包，数据
library(org.Hs.eg.db)
library(clusterProfiler)
library(pathview)
library(enrichplot)
data <- read.csv("limmaOut.csv",header=TRUE)
head(data)GSEA分析只需要两列信息，基因列和logFC（不同软件的差异分析这一列的名字会有差别）。
1.2 ID转化
基因名是symbol，将之转换为entrezid格式
gene <- data$SYMBOL
#开始ID转换，会有丢失
gene=bitr(gene,fromType="SYMBOL",toType="ENTREZID",OrgDb="org.Hs.eg.db")
#去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)#合并data 和 entrezid
data_all <- data %>%
inner_join(gene,by="SYMBOL")
dim(data_all)
head(data_all)
1.3 基因列表排序
将基因按照logFC进行从高到低排序，只需要基因列和logFC即可
data_all_sort <- data_all %>%
arrange(desc(logFC))
head(data_all_sort)
geneList = data_all_sort$logFC #把foldchange按照从大到小提取出来
names(geneList) <- data_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(geneList)#结果为排序的logGC，names为ENTREZID 345557 1308 55220 8788 6770 8745 2.425345 1.938050 1.831625 1.825617 1.786812 1.767371
二 GSEA分析
2.1 GSEA分析
需要gmt文件，http://www.gsea-msigdb.org/gsea/downloads.jsp路径下载，选择合适的
kegg_gmt <- read.gmt("c2.cp.biocarta.v7.4.entrez.gmt") #读gmt文件
gsea <- GSEA(geneList,
TERM2GENE = kegg_gmt) #GSEA分析
head(gsea)
其中：
- ID ：信号通路
- Description ：信号通路的描述
- setSize ：富集到该信号通路的基因个数
- enrichmentScore ：富集分数，也就是ES
- NES ：标准化以后的ES，全称normalized enrichment score
- pvalue：富集的P值
- p.adjust ：校正后的P值
- qvalues ：FDR （false discovery rate）错误发现率
- rank ：排名
- core_enrichment：富集到该通路的基因列表。
可视化-点图
dotplot(gsea)
2.2 GSEA-GO分析
gse.GO <- gseGO(
geneList, #geneList
ont = "BP", # 可选"BP"、"MF"和"CC"或"ALL"
OrgDb = org.Hs.eg.db, #人 注释基因
keyType = "ENTREZID",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",#p值校正方法
)
head(gse.GO,2)
2.3 GSEA-KEGG分析
gse.KEGG <- gseKEGG(geneList,
organism = "hsa", # 人 hsa
pvalueCutoff = 1,
pAdjustMethod = "BH",) #具体参数在下面
head(gse.KEGG)
#为方便展示，最后一列暂不展示
head(gse.KEGG)[1:10]
三 GSEA可视化
使用gseaplot2函数进行可视化
3.1 简单可视化
gseaplot2(gse.KEGG,
1) #展示第一个通路
可以进行一些调整以接近文献
1）修改GSEA线条颜色
2）添加P值的table
3）展示指定的通路
4）展示多个通路
5）只展示上两部分
3.2 展示指定通路
gseaplot2(gse.KEGG,
title = "Olfactory transduction", #设置title
"hsa04740", #绘制hsa04740通路的结果
color="red", #线条颜色
base_size = 20, #基础字体的大小
subplots = 1:2, #展示上2部分
pvalue_table = T) # 显示p值
3.3 展示多个GSEA结果
A：使用数字的方式
gseaplot2(gse.KEGG,
1:3, #绘制前3个
pvalue_table = T) # 显示p值
B：使用向量指定通路
gseaplot2(gse.KEGG,
c("hsa04740","hsa05310"), #指定通路向量
pvalue_table = T) # 显示p值
C：点 形式
gseaplot2(gse.KEGG, 1:5, #按照第一个作图 ES_geom = "dot", base_size = 20,
pvalue_table = T) 
应该和文献中的图很接近了，剩下的就是用自己的数据去尝试了。
四 参考资料：
http://yulab-smu.top/clusterProfiler-book/
PS：有个交流的讨论组，想沟通交流的，后台回复”入群“。
◆ ◆ ◆ ◆ ◆
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2021-06-22，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录