知识分享| 转录组个性化分析(2)—

上一期中给大家介绍了拟时序分析的意义及具体的分析过程，本期继续给大家带来转录组个性化分析——GSEA。废话不多说，干货直接奉上！

1 GSEA基本概念

基因集富集分析（Gene Set Enrichment Analysis，GSEA）：用一个预先定义的基因集中的基因来评估在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的贡献。

2 GSEA原理

(1) 背景基因排序：将全部基因按照某种指标（差异分析p值，表型相关性，表达量等）进行排序，比如log2FC排序。

(2) 目标基因富集：将某个特定类型的基因在排序表中标出，目标基因可以是某个通路或GO terms的基因等。

(3) 计算富集分数：使用加权法，计算ES值变化。对位于中部（与性状相关性低）的部分采用较小的权值，所以越集中在两端，与表型的相关性越高。ES曲线最大值为富集分数（Enrichment Score）。

(4) Permutation test：对基因集的ES值进行显著性检验及多重假设检验，从而计算出显著富集的基因集。

3 GSEA分析的作用

GSEA和常规的GO、KEGG的差异在于，GSEA使用的是基因集，传统的富集分析不需要考虑基因表达量的变化趋势，其算法的核心只关注这些差异基因的分布是否和随机抽样得到的分布一致，即使后期在可视化时，我们在通路图上用不同颜色标记了上下调的基因，但是由于没有采用有效的统计学手段去分析这条通路下所有差异基因的总体变化趋势，这使得传统的富集分析结果无法回答如下问题：一个富集到的通路下，既有上调差异基因，也有下调差异基因，那么这条通路总体的表现形式究竟是怎样呢，是被抑制还是激活？或者更直观点说，这条通路下的基因表达水平在实验处理后是上升了呢，还是下降了呢？

想要知道进行了差异分析的两组有什么功能和通路的差别，手上有大部分的功能分子以及对应的值，这个值可以是 logFC。可以用这个 logFC作为分子的排序，从而来评估在预先定义的基因集中是否显著富集。

4 GSEA分析代码

library(DO.db)
require(DOSE)
library(clusterProfiler)
library(AnnotationHub)
library(readr)
library("genefilter")
library(pheatmap)
library(tidyverse)
library(DESeq2)
library(ggplot2)
library(export)
library(enrichplot)
library(Rgraphviz)
#构造图片输出函数，need input filename width height
#函数依赖export包
out_img <- function(filename,pic_width=5,pic_height=7){graph2png(file=filename,width=pic_width,height=pic_height)graph2ppt(file=filename,width=pic_width,height=pic_height)graph2tif(file=filename,width=pic_width,height=pic_height)
}
#all_entrez.csv的第一列是entrezid,第二列是FoldChange的值。
#gse需要单独做数据格式
d <- read.csv("all_entrez.csv")
geneList <- d[,2]
names(geneList)=as.factor(d[,1])
geneList <- sort(geneList,decreasing = TRUE)#gseGO进行GSEA分析
#参考连接https://yulab-smu.github.io/clusterProfiler-book/chapter12.html
###gseBP <- gseGO(geneList=geneList,ont="BP",OrgDb=maize,keyType = 'ENTREZID',nPerm = 50000,minGSSize = 100,maxGSSize = 6000,pvalueCutoff = 0.05,verbose = FALSE)############# GSEA CC 模式 start
ego3 <- gseGO(geneList = geneList,OrgDb = maize,ont = "CC",nPerm = 1000,minGSSize = 100,maxGSSize = 1000,pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego3,file = "GESA-GO_CC.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego3)
out_img(filename = "ridgeplot_CC",pic_width = 12,pic_height = 12)
#只显示值最高的一组的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])#显示前4组信息
gseaplot2(ego3,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_CC",pic_width=12,pic_height = 10)#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA CC 模式end############# GSEA BP 模式 start
ego2 <- gseGO(geneList = geneList,OrgDb = maize,ont = "BP",pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego2,file = "GESA-GO_BP.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego2)
out_img(filename = "ridgeplot_BP",pic_width = 12,pic_height = 12)
#只显示值最高的一组的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])#显示前4组信息
gseaplot2(ego2,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_BP",pic_width=12,pic_height = 10)#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA BP 模式end############# GSEA MF 模式 start
ego4 <- gseGO(geneList = geneList,OrgDb = maize,ont = "MF",pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego4,file = "GESA-GO_MF.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego4)
out_img(filename = "ridgeplot_MF",pic_width = 12,pic_height = 12)
#只显示值最高的一组的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])#显示前4组信息
gseaplot2(ego4,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_MF",pic_width=12,pic_height = 10)#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA MF 模式end#gsaKEGG基因富集分析
kk2 <- gseKEGG(geneList = geneList,organism = 'zma',pvalueCutoff = 0.05,verbose = FALSE)
write.csv(kk2,file="GSEA_KEGG.csv")gseaplot2(kk2,geneSetID = 1:4,ES_geom = "dot",pvalue_table = TRUE)
out_img(filename="GSEA_KEGG",pic_width = 12,pic_height = 10)
ridgeplot(kk2)
out_img(filename="ridgeplot_GSEA_KEGG",pic_width = 12,pic_height = 12)

5 结果解读

典型结果图由上、中、下三个部分组成：

上：为富集评分的情况，如果 NES 为正，则峰出现在左侧（头部富集）（高表达组富集）基因集中核心分子主要集中在左侧高表达组中；如果 NES 为负，则尾部会出现谷（尾部富集）（低表达组富集），基因集中核心分子主要集中在右侧低表达组中。

中：每一根竖线代表基因集中一个分子，上传数据的分子根据给定的值进行排序，排序后单独提取当前基因集中的定义的分子，分子的位置情况即为中间部分的所示。

下：把上传数据分子给定的值进行归一化后的值进行可视化。下部分的结果可以不用怎么关注。

从该图中可以看出，这个基因集是在MUT这一组高表达的，在折线图中有个峰值，该峰值就是这个基因集的Enrichemnt score，峰值之前的基因就是该基因集下的核心基因。基因集下的所有基因在这个排序列表的顶部富集，我们可以说，从总体上看，该基因集是上调趋势，反之，如果在底部富集，则是下调趋势。

知识分享| 转录组个性化分析(2)——GSEA相关推荐

知识分享| 转录组个性化分析(5)——转录因子及转录因子结合位点预测
转录因子也称反式作用因子,在动植物的生长发育及其对外界环境的反应中起着重要的调控作用,已成为现在生物学研究领域的热点,其功能分析是重要的研究内容之一.因此我们可以通过转录因子注释,或者可通过表达量聚类 ...
知识分享 | 转录组个性化分析(4)——蛋白互作分析
蛋白质-蛋白质相互作用(protein-protein interaction, PPI)是指两个或两个以上的蛋白质分子通过非共价键形成蛋白质复合体(protein complex)的过程.在进行数据 ...
知识分享| 转录组个性化分析(6）——WGCNA
近年来,很多SCI高分文章中都使用了WGCNA,那么WGCNA是什么,它可以应用于哪些研究方向,又如何进行WGCNA分析,其分析结果如何看?现在就带着这些问题,跟着小易一起学习探讨吧! 1 WGCN ...
css前端知识分享—页面布局分析
今天分享下"css前端知识分享-页面布局分析"这篇文章,文中根据实例编码详细介绍,或许对大家的编程之路有着一定的参考空间与使用价值,需要的朋友接下来跟着云南仟龙Mark一起学习一下 ...
凌恩生物文献分享|转录组高级分析--植物抗性基因分析
植物在生态系统中无处不在并且占据重要地位,也是人类食物及各种产品的一个基本来源.而影响植物生存生长发育的很重要一个因素就是病虫害,全世界近40%的农作物产量因病虫害而损失,但是研究人员发现有些植物天然 ...
MPB：上海交大肖湘组分享基于基因芯片的海洋微生物转录组学分析技术
为进一步提高<微生物组实验手册>稿件质量,本项目新增大众评审环节.文章在通过同行评审后,采用公众号推送方式分享全文,任何人均可在线提交修改意见.公众号格式显示略有问题,建议点击文末阅读原文 ...
高级转录组调控分析和R语言数据可视化第十三期（线上线下，7月底开课）
福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现决定安排扩增子16S分析.宏基因组.转录组线上直播课.报名参加线上直播课的老师可在365天内选择参加同课程的一次线下课 .期待和大家的 ...
有了易生信，导师再也不用担心我的单细胞转录组整合分析啦
2019年10月9日,单细胞转录组再等Nature.题为Decoding human fetal liver haematopoiesis的研究,对受孕后4周至17周的人胚胎肝脏.卵黄囊.肾脏和皮肤组 ...
二代三代转录组测序分析实战班
本文原创"生信宝典"公众号,作者陈同. 转录组大家都很熟悉了,我们之前也有几篇介绍: 转录组分析的正确姿势 39个转录组分析工具,120种组合评估(转录组分析工具哪家强-导读版) ...

知识分享| 转录组个性化分析(2)——GSEA

知识分享| 转录组个性化分析(2)——GSEA相关推荐

最新文章

热门文章