背景
之前的一些推文,大部分收录专题于生物信息学,目的是帮助大家入门生物信息学的领域。本次开设新专题,“富集分析”,了解富集分析的各种手段,学会十八般武艺。
不过,因为公众号没有留言功能,大家交流学习讨论不知从“何”下手,为了解决这个问题,大家可以移步论坛发帖。链接如下:
代码语言:javascript复制http://bioinfoer.com/
一、GSEA 基本介绍
GSEA 全称是 gene set enrichment analysis 基因富集分析,是博劳德研究所 broad institute 研究团队开发的一个针对全基因组表达谱数据进行分析的工具,免费注册后即可进行下载。
它能够在大量基因或蛋白质中识别与疾病表型有关的过表达基因或蛋白质,从而判断某项干预与某一表型的关系。
GSEA 分析所适用的主场景之一:它能帮助生物学家在两种不同的生物学状态中,判断某一组有特定意义的基因集合的表达模式更接近于其中哪一种。因此 GSEA 是一种非常常见且实用的分析方法,可以将数个基因组成的功能基因数据集与测序及芯片得到的全部数据做出简单而清晰的关联分析。
gsea特点:
1. 无需差异基因;
2. 分析的是基因集合而非单个基因(GO)或少数基因(Pathway);
3. 富集分析。
怎么理解这个富集分析?
想要理解它首先要知道单基因分析,对实验组和对照组进行高通量测序或基因芯片检测获得的数据直接进行比对分析,发现基因表达发生了变化,到此为止就是单基因分析,单基因分析未考虑基因间的相互作用,因此很难对基因的表达变化做出解释。那么,将获得的两组数据进行一定处理后与按先验知识归类的基因集合比对分析,将某个干预和某个生物学功能变化联系起来,这个过程就叫富集分析。
4. 将基因与 MSigDB 中预定义的基因集合进行比较。使用GSEA分析结果发表文章时注意引用网站上的文献。
二、分析前准备
进行分析之前需要准备 3 个文件:表达数据集、样品分组信息和基因数据集。
表达数据集是测序或芯片获得的表达谱信息、样品分组信息是自己构建的文件、基因数据集是已知功能的基因集合。
表达数据集按表达丰度排列,也就是上图中热图所展示的。功能基因数据集中出现在表达数据集当中的基因所处的位置用黑色竖线表示。
富集分数ES是从排序后的表达数据集的第一个基因开始,如果表达数据集中的基因出现在基因数据集中则加分,反之,不在基因数据集中则减分。所以,富集分数ES是动态变化的。这条绿色曲线正是富集得分的体现,绿色曲线的峰值就是最大富集分数。若 ES 为正值说明在顶部富集如 A,若为负值说明在底部富集如 B。
若研究的基因数据集的成员显著聚集在表达数据集的顶部或底部,说明基因数据集中的基因在表达数据集中高表达或低表达;若随机分配说明表达数据集与基因数据集对应的表型无关。Leading-eage subset 即最大 ES 值前对应的基因,即对富集得分贡献最大的基因成员。
三、关键概念
富集得分 ES,在计算过程中不断加分或减分,因此它是个动态的数值,一般的数据图片中往往不呈现 ES 值,有的话也应该是最大 ES 值。
校正后富集得分 NES,是个常数,用来比较表达数据集在不同功能基因数据集中的富集程度。
名义 p 值,因为它没有进行功能基因子集大小和多重假设检验校正,因此称为名义 P 值,它描述的是针对某一功能基因子集得到的富集得分的统计显著性,p 值越小富集性越好。
错误发现率 FDR,该指标进行了功能基因子集大小和多重假设检验矫正,用于判断假阳性率。
总体错误率 FWER,用于检验出尽可能多的候选变量的同时将错误发现率控制在一个可以接受的范围。
Leading-eage subset 即最大ES 值前对应的基因,即对富集得分贡献最大的基因成员(又称为核心基因)。
四、与传统富集分析的区别
GO 富集分析通过分析差异基因在生物学过程,分子功能、细胞组成中的富集定位,从而对基因进行注释和分类,它通过设定 cut-off 值选出差异表达基因,对它们进行 GO terms 富集度统计学分析,计算出差异基因 GO terms 的 p-value 及 p-value 的 FDR 值(q-value),定位差异基因最可能相关的 GO terms,也就找出了该组差异基因最相关的功能或生物学过程或细胞定位等。
KEGG 通路分析和 GO 富集分析类似,选出差异基因,通过统计学分析判断差异基因可能和哪些通路相关。这两个分析方法都需要筛选出差异基因,忽略对结果有贡献但没有落在差异显著范围内的基因。而 GSEA 是利用测序或芯片获得的全基因组表达谱进行分析,不需要指定差异基因阈值,得出的结果更加可靠。
我们看上面这个图,基因数据集显然在表达数据集高表达区富集,而进行 GO 分析时通过 P 值或矫正后的 P 值 FDR 值筛选后只剩下左右两边少量的差异表达基因,与这里显示的基因数据集进行比对时显然没有明显富集,因此就会漏掉部分富集的数据集。
五、GSEA 的操作流程
分为 4 部分
软件下载:GSEA 是基于 JAVA 环境运行,因此保持 JAVA 为最新版本是 GSEA 软件运行的基础。
数据准备与结果解读下文进行详细阐述,
参数设置与软件运行具体内容详见后续课程。
结果分析是第四部分。
准备格式无误的文件是富集分析成功的关键,需要准备 4 种文件,分别为表达数据集、样品分组信息、基因数据集,基因芯片注释。
表达数据集,在生信时代 GEO、TCGA 等数据库是丰富的数据来源,但其本质都是通过各种测序或芯片平台获得的数据,格式多种多样,均可被 GSEA 识别。以 GCT 文件为例,excel 表头以#1.2 为固定格式出现,占据第一行第一列,第二行第一列是基因个数,第二行第二列为样本数,基因 ID 根据测序或芯片平台而有不同,需要在数据分析参数选择界面选择匹配的平台。Description 为基因描述,可空着。
样品分组信息,是我们自制的一个表,告诉 GSEA 分析软件我的数据集里有多少样本以及是如何分组的,使用 EXCEL 构建,按照下图所示构建即可。
功能基因数据集是某一特定功能/表型所包含的所有基因的集合,用来判断表达数据集是否有某种功能聚集最重要的文件;GMT 文件则需要在 GSEA 网页中下载,MSigDB将基因分为各种子集,各取所需。需要注意的是 MSigDB 中只包含人的基因序列。
链接:
代码语言:javascript复制https://www.gsea-msigdb.org/gsea/msigdb/
基因芯片注释:不同测序或芯片平台所用的基因代码不一样,因此需要基因芯片注释文件来说明基因代码究竟是哪个基因,可以不用下载,GSEA 分析参数选择时直接选择相应的平台,如果你的表达数据库是从 GEO 等数据库下载来的,选择平台时注意和 GEO 中的 platform 一致。
六、结果解读
分析完成后会获得一个文件夹,里面有许许多多的图表。其中有个html文件是总体的分析结果,
红圈5 超链接里点进去可以看到 12 个高表达功能基因子集基因在该组别中的位置,以及详细的重要参数等。
接下来聚焦到具体每个功能基因数据集的结果:
当|NES|>1,P 值<0.05,FDR<0.25 这三者同时满足时,结果才有意义,一定要注意名义 P 值没有经过矫正而 FDR 值经过了功能基因子集大小和多重假设检验矫正,因此这两者都要关注。
下图是 GSEA 分析中最核心的图,由三个部分组成
最后是一个文献实例
如图,但该图最好要体现 NES,P 值,FDR 三个指标,使读者对该结果更加信服。