数据分析、矩阵运算基因富集分析

2022-04-17

2022-04-17  本文已影响0人  玄武Ken

GSEA缩写

GSEA的全称是Gene Set Enrichment Analysis, 中文翻译就是基因富集分析,最早提出GSEA的文章是下面的这篇文献:

基因富集分析是分析基因表达信息的一种方法,富集是指将基因按照先验知识,也就是基因组注释信息进行分类(注:我个人理解先验知识就是先于经验的知识,例如我们从书本上知道地球是圆的,也就是先验知识,我自己没有亲自去测量,也无法观察到,就知道地球是圆的。基因组注释就是对碱基序列的注释,给出一段碱基序列,根据前人研究的结果,我知道它是什么基因,发挥了什么功能,以上就是对先验知识的理解)。最早提出GSEA的文章是下面这篇文献:

1. Subramanian, A., et al. (2005). "Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles." Proceedings of the National Academy of Sciences 102(43): 15545-15550.

常规分析的局限

使用芯片或测序技术已成了生物研究的常规手段,获取某个样本的整个表达谱已经不再困难,但是,如何从这些表达谱中发掘有用的信息,以及对生物学现象进行解释则是一个挑战。在一个典型的研究中,例如我们研究某个抗肿瘤药物对于两种肿瘤细胞(一个对药物敏感,一个不敏感)的效应,然后检验这两类肿瘤细胞的表达谱,可以获得数千个基因的变化信息。我们根据这两类细胞之间基因表达差异的情况,例如根据它们之间的差异倍数,对它们进行排序,得到一个列表L,就可以从中获得一些信息。常规的做法就是关注那些表达量差异极大的基因,也就是列表L顶部或底部的基因,用于发现这两类细胞不同变化的一些线索,不过这种方法有一些局限。

  1. 经过多重假设检验校正,或许没有一个基因满足统计学显著性的阈值,因为相关的生物学差异相对于实验误差来说,并不明显,也就是不好区分出基因的差异与实验误差。
  2. 也许会有很多差异基因,但是这些差异基因无法与生物学内容整合起来,也就是说用这些差异基因不太好解释生物学现象。
  3. 单基因分析方法也许会错过许多对通路有重要影响的基因。细胞过程经常会协同影响一批基因。例如,编码一个代谢通路的一批基因的表达量都升高20%,就会造成非常大的影响,这种影响有可能比一个基因表达量增加20倍都有效。
  4. 当两个不同团队对同一个生物学体系进行研究时,他们的到的差异基因或许就不太一致。

GSEA就是为了解决上述的问题而提出来的。此方法可以在基因集(也就是多个基因)水平上来处理表达谱数据。这里所说的基因集(gene sets)是指基于当前积累的知识,例如关于生物通路知识或以前的到的共表达数据来定义的一组基因。GSEA主要研究一个基因集S(这里的S指的是我们定义的基因集)出现在目标基因列表L(也就是实验得到的按表达倍数排列的基因列表L)的顶部或底部是否与实验中的不同表型有关。

因此,我们就可以知道,GSEA分析有三个特点:1、分析的基因集合而不是单个基因;2、将基因与预定义的基因集进行比较;3、富集分析。
与GO(Gene Ontology)和KEGG pathway分析相比,GSEA分析的主要优势在于:

一般的差异分析(GO和Pathway)往往测重于比较两组间的基因表达差异,集中关注少数几个显著上调或下调的基因,这容易遗漏部分差异表达不显著却有重要生物学意义的基因,忽略一些基因的生物学特性、基因调控网络之间的关系及基因功能和意义等有价值的信息。而GSEA不需要指定明确的差异基因阈值,算法会根据实际数据的整体趋势,为研究者们提供一种合理的解决目前芯片分析瓶颈问题的方法,即使在没有先验经验存在的情况下也能在表达谱整体层次上对整条基因进行分析,从而从数理统计上把表达谱芯片数据与生物学意义很好的衔接起来,使得研究者们能够更轻松,更合理的解读芯片或测序结果。

GSEA的原理

首先是对每个样本里面的基因表达值在样本内部进行排序,这里可以理解为,实验组与对照组的差异基因进行排序。但是,这个差异如何量化,有多种方法,可以是Signal2noise,或者是T检验的值,或者是flod change的值,或者是logFC等。而GSEA默认的就是signal-to-noise metric来对基因进行排序。如果要想计算Signal2Noise,每个group必须要有3个及以上的samples。

现在我们定义一下各种参数。

其中我们把自己测出来的差异基因的排序列表称为目标基因列表L,把根据先验知识预先定义的基因集称为功能基因集S(这个基因集可以是与免疫有关的,可以是与某个代谢通路有关的,也可以是与miRNA有关的,具体的看GSEA官网),把这个基因集中的成员称为s(前者是大写,后者是小写)。GSEA的运行原理就是判断功能基因集S里面的成员s目标基因列表L里面是随机分布的,还是聚集在目标基因列表L的顶部或底部。如果我们研究的功能基因集S的成员显著聚集在目标基因列表L的顶部或底部,那么就说明此基因集成员对实验干预造成的结果有作用哦,就是我们要关注的基因集。

以下是GSEA计算中的几个概念。

image

其中目标基因列表L表示的就是差异基因已经排序后的序列,S指的是功能基因集S,我们看左边的Ranked List,它右边是FC,即表达倍数(明显能看出来是从大到小的排列),红色表示的是Hit,也就是说这个基因在功能基因集S里面,如果在,就加分,也就是说在右边的表格里的Contribution to running sum for ES加上相应的分数。如果是蓝色,表示这个基因不在,就减去相应的分数。

GSEA分析所需要的文件

GCT文件

GCT格式室友分割符分开的一个数据文件,它记录的是基因表达的数据集,也就是我们要尽兴GSEA计算的原始数据,它的格式如下图所示:


image

从左到右,顺时针看到下,GCT文件包含以下这些信息:

  1. Always "#1.2" :这个单元格一直标注的是#1.2,不用管它。
  2. # of samples :这个是样本的数目,在上面的表格中,这个数目是130。
  3. Third column onwards are sample names. These must be UNIQUE, 它们是样本的名称,样本名称必须保证是唯一的,不能重复出现。
  4. Data stars on line: 第四行,这是每个基因具体的数值。
  5. Each column contains expression values from 1 sample, Missing values are allowed(leave empty): 基因的表达数值可以是空。
  6. Column 2 : Row descriptions.: 第2列,对基因的描述,可以为空(填上NA)。
  7. Column 1 : Row identifiers: 这一列是基因的标识,通常是探针集的ID(probe set ids)或克隆号,此列数据必须唯一,不能出现重复。
  8. A2单元格: The # of rows:探针的数目。

CLS文件

CLS的全称是Categorical (e.g tumor vs normal) class file format, 它表示的是分组的信息,下面是一个CLS文件的内容:


image

以上面的图片为例说明,第1行有3个数字,分别为38,2,1.其中38是样本的数目,2是分组的数目,1不用管。第2行,分别是ALL和AML,这表示分组的信息,即有2个分组,分别是ALL和AML。第4行,分别是0和1,0表示ALL,1表示AML,它们的数目表示各自的样本有多少个。再看一个案例:


image

在这个案例中,我们可以看到一共有9个样本,分了3组;这3组分别是KRAS_MUT, WT, MYC_MUT, 第3行是具体的分组数目。

常用的定义基因集的文件主要有GMT与GMX。

GMT文件

GMT文件的全称是Gene Matrix Transposed file format, 中文是"基因矩阵转换文件格式",缩写为GMT。GMT文件是由制表符分割的文件,它用于描述基因集,可以理解为高通量测序的注释文件。下图是GMT文件的一个案例,每行代表一个基因集:


image

在上面的图片中,我们可以看到,第1行表示的一个是基因集,第1列表示的是基因集的名称,这一列不能出现重复的基因集。第2列是对基因集的一个简单描述,可以选择不填(写上NA),第一列的长度可以不等,如果要在Excel表格中编辑GMT文件,一定要注意Excel的自动纠正功能(例如,Excel有可能吧SEP8基因自动纠正为8-Sep)。

GMX文件

GMX文件是另外一种描述基因集的文件格式,下图是GMX的一个案例,它与GMT格式正好相反,它的每列代表一个基因,如下所示:


image

RNK文件

RNK文件的全称是Ranked list file format,这个文件是一个有序的基因列表,第1列是基因名称(A1单元格忽略),第2列是数值(可以表达倍数差异,也可以使用类似于t检验进行的排序),如下所示:


image

GSEA使用案例

这个案例使用的是GSEA官网的数据,数据是p53突变型与野生型的肝细胞系的表达谱,这里使用的是芯片的表达数据,芯片现在已经不怎么用了,这里只是讲一下GSEA官网这个软件的使用。

第1案例

第一步,去官方地址下载GSEA软件,如下所示:

image

下载后,打开,界面如下所示:


image

GSEA的官网提供了一些数据的案例,我们可以拿来练习一下,地址点击这里,这些数据如下所示:

image
这里以第2行的p53为例说明一下,这一行中有三个数据文件,分别为P53-hug95av2.gct, P53_collapsed.gctP53.cls,都下载下来。其中GCT文件为基因表达文件,其中带有collapsed的文件为表达谱文件中的标记符(通常是一串数字)已经被基因名(symbol)替换掉了,另外一个文件是原始文件,如下所示:
image
使用Excel打开P53-hug95av2.gct这个数据文件(其实也可以使用Notepad++ PilotEdit Lite这类文本编辑器来打开GCT文件),这个文件是分级以及表达量的数据,如下所示:
image

现在用Excel打开P53.cls这个文件,这个文件是分组信息,如下图所示:

image

从注释信息中我们可以得知,有50个样本,分为两组,分别为2和1表示这两组,2表示MUT组,1表示WT组。将上述数据导入到GSEA中,过程如下所示:
打开Load data, 如下所示:

image

此时右侧出现新的页面,打开Browse for file,找到原始数据文件夹,如下图所示:

image

导入要分析的文件,分别为cls,gct文件(gmt文件后面可以直接在软件中下载),如下所示:

image

此时,单击左侧的Run GSEA,如下所示:

image

此时出现了一些参数选择页面,如下所示:

image

在第一行Expression dataset中选择表达的数据文件,即P53_collapsed_symbols.gct

在第二行的Gene sets database中选择基因集文件,即gmt文件,单击右侧省略号,此时出现基因集文件的选择框,如下所示:

image

根据官网信息,选择c2.all.v7.0.symbols.gmt这个基因集,如下所示:

image

第三行, Number of permutation, 此处设为10(我开始设置为1000,没跑出来);

第四行,Phenotype labels,选择MUT_versus_WTWT_versus_MUT,顺序不同会导致最终生成的ES图不同,我开始的时候,设为MUT_versus_WT, 发现ES得分的图是倒的,不好看,现在就设为WT_versus_MUT

第五行,Collapse dataset to gene symbols, 此处填false,这个参数用于说明是否要把探针的编号转为基因ID(gene symbols),由于P53_collapsed.gct这个文件本身就已经是使用基因ID了,因此不需要转换,设为false;

第六行,Permutation type: 填入 phenotype。

第七行,Clip platform, 芯片类型,前面Collapse dataset to gene symbols已经是false了,因此这里不用选,再看第二个页面,也就是Basic fields,填入要分析的样本名称,保存结果的文件件,如下图所示:

image

填完之后,单击最下面的Run,如下所示:

此时,GSEA开始运行,如下所示:

image

运行结束后,系统有提示,如下所示:

image

现在进入保存结果的文件夹,如下所示:

image

其中,我们选择index这个html文件,就可以看到分析整体结果,如下所示:

image

打开,如下所示:

GSEA结果解读

image

第1部分

现在分别查看这些结果的信息,先看第一部分:

image

这一部分,从上到下,可以看到:

image

现在单击第一张图,就可以查看富集分析的结果,如下图所示:

image

再回到主页面,如下所示:

image

单击enrichment results in html就可以查看所有的富集分析结果,如下所示:

image

现在把标题栏放大,如下所示:

image

其中,GS follow link to MSgnDB表示具体的功能基因集S的名称;
GS DETAILS 可以查看具体的某个功能基因集S的富集分析结果;
SIZE表示基因集里的基因数;
ES表示富集分数的ES值;
NES表示校正后的ES值;
NOM p=val是名义p值;
FDR q-val用FDR法校正后的p值,即Qvalue;
FWER p-val用FWER法(Bonferonni校正)后的P值;
RANK AT AMX表示ES值达到最大进度,对应的通路基因的排名;
LEADING EDGE,显示了用于定义Leading edge subset的一些参数,分别有Tags, List, Signal。

where N is the number of genes in the list and Nh is the number of genes in the gene set.If the gene set is entirely within the first Nh positions in the list, then the signal strength is maximal or 100%. If the gene set is spread throughout the list, then the signal strength decreases towards 0%.

现在单击一下GS DETAILS 这一栏中的Details..,就会看到某个基因集的详细结果,如下所示:
拉到中间区域,就会看到这样的结果,即:

image

现在解释一下这个列表:

如果在GSEA软件中运行Leading Edge Analysis,就会出现下面结果,如下所示:

image

其中,先看第1张图,如下所示:

image

这是一张热图,它显示了Leading edge subset中基因的情况,在热图中,基因的表达值按5种颜色进行区分,分别为red,pink, light blue, dark blue,它们代表基因表达水平为high, moderate,low,lowest。

再看第2张图,即右上,这是Set-to-Set图,如下所示:

image

这些颜色显示了subset之间的重叠情况,如果颜色块越深(仔细看,有些颜色非常浅),表明subsets之间的重叠性越高,一个小方格中A与B的密度对应了X/Y的比例,其中X是A中leading edge genes的数目,Y表示, A与B的并集。颜色越深,表明,A与B含有相同的Leading edge genes越多,如果方格是白色的,表明A与B不含相同的Leading edge genes。

再看第3张图,即左下的图,如下所示:

image

这张图的最底部是基因的名称,纵坐标表示,含有这个基因的功能基因集S的数目。

再看第4张图,也就是右下图,如下所示:

image

这个图是直方图,其中横坐标Jaccard表示,每对leading edge subset的交集除以并集(可以看到,它就是个小数,毕竟交集的数目要小于并集的数目)。Number of Occurrences表示在一个置顶的bin内,配对的Leading edge subset数目,在这个案例中,大多数subset配对的Jaccard结果都位于0到0.02之间,可以说,大部分配对的Subsets都没有重合的部分。

下面是原文:

The last plot is a histogram, where the Jacquard is the intersection divided by the union for a pair of leading edge subsets. Number of Occurrences is the number of leading edge subset pairs in a particular bin. In this example, most subset pairs have no overlap (Jacquard = 0).

GSEA的主要观察指标

GSEA主要看四个指标,分别为ES,NES,名义p值,FDR值,现在分别解释一下。

image
  1. ES,全称为Enrichment Score,即富集分数,一开始都是0,就是绿色长蛇状粗线,这个绿色的曲线救市ES的值计算过程,从左到右每遇到一个基因,就计算出一个ES值,连成线,最高峰为富集得分(ES)。
  2. 中间的长得像条形码的图形,图片中长得像条形码的东西就是功能集S在目标基因列表L中所在的位置。
  3. 红色直线是绿色曲线的顶点,它左侧的数据是Leading edge subset,也就是对富集得分贡献最大的基因成员,若富集得分为正值,则表示这基因集位于L的顶部,如果富集得分为负值,则位于L的底部。
  4. 当用功能基因集S从上到下,遍历排序好的目标基因列表L时,此时最下面的灰色区域就是不同基因排序的结果,它与分组情况有关,排序的结果如果从正值到负值进行排列,正值是与第1个分组有关(例如WT组),负值与第2个分组有关(例如MUT组),它是按从高到低进行排序的基因次序,排序用的是Signal2noise。

这是校正后的ES数值。由于ES的计算原理是,一个目标基因列表L中的基因是否在一个功能基因集S中出现来计算的,但是各个功能基因集S中所包括的基因数不同,而且不同的功能基因集S与要分析的目标基因列表L的相关性也不同,置换次数也不同,因此在比较数据集在不同功能基因集中的富集程度,需要对ES进行标准化处理。

通过校正原始ES值,NES可以解释功能基因集S大小的差异,以及功能基因集S目标基因列表L之间的相关性,因此NES可以用于比较基因集之间的结果,GSEA定义的NES公式如下所示:
[图片上传失败...(image-47ce29-1650292588849)]%7D)
NES基于所有数据集置换(permutation)的基因集富集分数进行计算的,因此通过改变置换方法,及置换数目,或者是基因表达数据量的大小,都能影响NES。例如,现在我们看2个分析:
第一,你分析一个表达数据集,GSEA会产生一个有序列表,并对这个列表进行分析;
第二,你使用GSEA对上面这个有序列表重排序,并分析。
如果你使用相同的参数设置(这里的相同是指相同的基因集S,相同的目标基因列表L),那么你得到的ES就是相同的;但是如果你在两次分析时,设定了不同的置换次数(permuation)的话,那么校正后的富集分数,即NES则不同。这两次分析能反映出用于置换的数据集的很大不同(表达数据集与排列的基因列表),如下所示(第1次的置换次数是149,第2次是10):

image

名义p值描述提针对某一个功能基因集S得到的富集得分的统计显著性,p值越小,越能说明基因的富集性越好。但在GSEA分析结果的四个参数中,只有FDR值进行了功能基因集S大小和多重假设的校正,而p值没有,因此,当分析结果中出现一个高度富集的功能基因集S,具有很小的名义p值,却有很大的FDR值时,往往意味着其实和其他功能基因子相比较,它的富集并不是很显著,原因可能是用于分析的芯片数据样本量较少,或者杂交信号微弱,又或是选择的功能基因子集并没有很好的反映样本的物理学意义。在GSEA的报告中,当p=0时,这就意味着世纪的p值小于1/随机组合的次数,例如,当选择随机组合数为100时,报告中的p=0实际上是说p值小于0.01,所以提高随机数,就能得到更精确的p值,通常情况下,运行GSEA的随机组合次数都会选择1000,但不要超过1000(个人电脑一般选10),否则系统内存会被耗尽。

通常情况下,我们一般认为|NES|>1, NOM p-val<0.25的通路下的基因集是有意义的。

现在再回到主页面,如下所示:

image

Detail enrichment results in excel format (tab delimited text)这一行是使用Excel表格呈现的富集分析结果;

最后一行,即Guide to interpret results,则是对结果的说明。

第2部分

再看结果的第二部分,即如下所示:

image

这是对MUT组的富集分析结果,内容与WT组的差不多。

第3部分

现在看第3部分,如下所示:

image

这里的信息提示,一共分析了10100个基因;

第4部分

看第4部分,如下所示:

image

这是功能基因集S的一些信息,第一行提示,根据min=15,max=500的标准来筛选,1329个基因集中,不能用的有430个,能用的有899个;

第二行信息说明,能用的基因集是899个(WT组有421个,MUT组有478个);

第三行说明,基因集的内容,也就是不同的功能基因集中有多少个基因。

第5部分

再看第5部分,如下所示:

image

从这一部分,我们可以看到这些信息:

第1行:一共有10100个基因;

第2行:WT组的基因是5707个,占56.5%,它所点相关关系面积的59.8%(这个59.8%我的理解就是在ES富集分析图中,最下面的那个灰色区域);

第3行:MUT组的基因是4393个,占43.5%,它所点相关关系面积的40.2%;

第4行:目标基因列表的排序,是个EXCEL表格,打开,如下所示:

image

一共是10100个,最后一部分如下所示:

image

第5行: 热图和相关关系的示意图(只显示了前50个基因),其中,热图用5种颜色来表示基因表达水平的高低,即red, pink, light blue, dark blue分别表示high, moderate, low, lowest, 如下所示:

image

(热图内容:Heat Map of the top 50 features for each phenotype in P53_collapsed_symbols.P53.cls#WT_versus_MUT)

image

(图片内容:Fig 2: Ranked Gene List Correlation Profile Ranked list correlations for P53_collapsed_symbols.P53.cls#WT_versus_MUT)

此外,在GSEA的使用的手册上发现,还有一种蝴蝶图(Butterfly plot)不过,在我的这个分析中,貌似没有找到,蝴蝶图显示的是基因顺序和排序度量得分(ranking metric score) 之间的正相关和负相关的关系。默认情况下,蝴蝶图只显示排名前100的基因,这也就是说,它只显示目标基因列表前100个和最后100个基因。不过,可以在Run GSEA Page这个页面中设置Number of marker参数来设置显示的基因数。富集图下面的部分显示了观察到的基因排列和Ranking metric score 得分之间的相关关系,蝴蝶图显示了观察到的相关关系,以及置换(1%,5%,50%)正相关,负相关的关系。蝴蝶图提供了一些哪些功能基因集的置换能够改变基因排序和Ranking metric score的信息,如下所示(这张图是从使用手册中截取的,不清楚):

image

第6部分

再看第6部分,如下所示:

image

第2行显示了全局的ES值,如下所示:

image

第7-8部分

剩下的内容如下所示:

image

Other则是这次分析的一些参数,Comments作为随机种子的时间戳。

参考资料

1. GSEA分析一文就够.Jimmy.生信技能树

2.GSEA富集分析 - 界面操作

3.基因集富集分析

4.GSEA分析是个什么鬼(上)

5.GSEA User Guide

6.Quick Tour of the GSEA Java Desktop Application

7.用GSEA来做基因集富集分析

8.使用PGSEA包进行基因集分析

9.制作自己的gene set文件给gsea软件

10.Quick Tour of the GSEA Java Desktop Application

11.冯春琼, 邹亚光, 周其赵,等. GSEA在全基因组表达谱芯片数据分析中的应用[J]. 现代生物医学进展, 2009, 9(13):2553-2557.

上一篇下一篇

猜你喜欢

热点阅读