GSEA教程

写在前面的话：GSEA对于生物信息学技能的要求不高，以至于不需要太多的生信知识就能学会，这也成为了生物研究生的必备技能，如果你还不会的话，这篇教程将让你100%学会，只要你从头到尾跟着我的步骤走。

我按照网上的教程学习GSEA的时候绕了不少弯路，参考了很多不同的教程，并且不断地Debug，花了一周时间才跑出结果。之后我再做GSEA就如行云流水般顺畅。这篇教程，保证看完包会，文末我甚至给出了可以直接拿来跑的文件，不耽误大家时间，100%能重复出结果！

下面进入正题：

GSEA（Gene Set Enrichment Analysis）,基因集富集分析，与GO分析类似，是一种基因差异表达分析。但是在原理上，与GO分析有很大区别。GO分析是针对差异基因，即看差异表达的基因在哪些GO分类或者pathway中富集，从而说明这些差异基因影响了哪些pathway。GSEA分析是针对所有基因的表达（不能用差异基因来跑GSEA,否则GSEA就失去了意义），这些基因的表达矩阵称为dataset。dataset一般是经过处理的RNA测序结果的基因表达矩阵，或者是从数据库中下载得到的表达矩阵。这些表达矩阵必须手动分组，比如分为control组与treat组。然后GSEA会自动分析出control组与treat组之间每个基因的差异并按照差异大小给每个基因算分，然后GSEA会根据这些基因的分值以及他们在某些基因集（geneset)中的富集程度来计算geneset富集分数（EnrichmentScore)并绘图（geneset由GSEA软件自带也可自行到数据库下载（后面由链接）），这样就得到了geneset的富集信息，从而推断出control组与treat组究竟在哪些geneset所代表的通路中存在差异。以上为GSEA简单的原理，算法什么的我不太懂，也不需要去了解，只要会用GSEA进行数据分析就可以了。

一. GSEA的前期准备

1.电脑一台，华为matebook，win10，8G内存。

2.GSEA软件：

http://software.broadinstitute.org/gsea/downloads.jsp

右边选择内存大小，最好选小于电脑内存的最大内存，比如我电脑8G内存，我就选择4GB(for 64-bit Java only)。内存大小决定了GSEA能处理的数据大小以及计算的精细程度。这里有两个坑我曾经栽过——请注意：第一是Java环境是必须的，请到java官网下载64位java8 （https://www.java.com/zh_CN/，官网只会推荐最适合你系统的java，没有其他版本的下载选项，但这不一定是我们所需要的版本！）或者用我文末分享的网盘链接进行下载，GSEA软件不支持java9或更高版本。第二是：如果选择2GB以上内存的java需要64位系统，如果你的电脑是32位系统，请选择1GB，win7-8系统一般会安装64与32两个系统，而win10不再有64与32位的区别，win10用户可放心下载。点launch进行下载，下载完成后可直接运行，如无法运行或没有相应请检查java是否安装正确。

下载完成后的图标为下图所示，如果图标未显示，请检查java是否安装正确。

这样，软件就下载完成了。

3. 实验数据（dataset）的准备。

NCBI的GEO数据库收录的绝大多数已发表文章的测序数据，而且下载起来很方便。本次使用的是人肝脏单细胞测序的数据，GEO编号GSE115469。

由于该测序数据集太大，EXECEL无法打开，需要用R语言进行整理，所以我将样本数缩小为100个，基因缩小为5000个，重新构建了一个表达矩阵，可以直接在EXCEL里进行排版（注意啊，我是方便教大家用excel操作才这样搞的，实际操作千万不能随意缩小数据啊，如果有几万基因几十个样本就用R处理）。这个矩阵我放在网盘里，文件名“human liver dataset”，文末有下载链接。打开之后如下图所示，第一列是基因名，第一行是样本名。这种基因名（ID)的专业术语叫做Gene symbol，我们平时接触到的Actin等基因名都是Gene symbol，此外还有很多其他的基因名，比如Ensemble ID，不同的测序平台对基因的编号也不同，这里不详细讨论Gene ID的转换问题，本篇主要带大家重复一遍GSEA的流程。下图是normalize后的单细胞测序表达矩阵，可用作GSEA的dataset。事实上，如果不是单细胞的数据，普通转录组测序的数据样本数一般只有10个以内，基因数至少2万个，完全可以用EXCEL进行处理。

对dataset的下一步处理是分组。我们这次就研究一下APOA2基因高表达的细胞与低表达的细胞有什么差异吧，首先我们将数据按照APOA2的表达水平排序。搜索到APOA2，在1454行

选择除基因名外的所有列，然后点击工具栏的“排序”

在新的对话框中，点击选项，选择按行排序，确定。

选择主要关键字“更多行”

点击1454行标，点击确定。

然后升序排列，点击确定。

排序后的APOA2

我们在心中把前10列定义为APOA2low，把后10列定义为APOA2high，保存EXCEL。

下一步，将.csv格式的ECXEL转换为.gct格式——可被GSEA软件识别的格式。这里使用既可以绘制热图又可以转换格式的在线软件：Morpheus。

https://software.broadinstitute.org/morpheus/

界面如下

然后上传上文处理好的表达矩阵。

点OK，会出现一个热图，无视这个热图，对我们没有帮助，我们是利用这个软件转换格式而不是画热图。然后点击工具栏File-Save dataset

命名文件，切记选择GCT1.2, 点OK，然后保存文件到本地

到这里我们对数据集(dataset)的处理就全部完成了。接下来进行最后一步准备工作。

4. 数据的Annotation文件的准备。

打开excel，输入

第一行20表示20个样本，2表示分为2组，1为识别字符（必须为1）

第二行表示分组信息，分为APOA2low，APOAhigh两组，注意#号不要丢。

第三行对应于dataset中每一列的分组，前10列输入APOA2low，后10列输入APOA2high。

然后点击保存，文件名结尾加上.cls，可被GSEA识别，储存为制表符分隔的文本文件（.txt)

打开文件检查一下，应该是这样的。

至此，GSEA前期准备就完成了。

下一步运行GSEA软件。

二. GSEA软件的运行

上传文件，点击load data, 使用Method1依次上传.gct与.cls结尾的txt文件。

成功会有以下提示，如果失败则可能是文件格式不对或者数据的排列不对。

然后我们成功上传了两个文件，如下图。

接下来我们设置GSEA output的储存位置，file-preferences。根据自己习惯设置。

2.设置参数。

我们进入run GSEA界面

第一行选择上传的dataset，第二行genesets选择c5.bp（biological process)。genesets详细介绍见MsigDB:http://software.broadinstitute.org/gsea/msigdb.

下一个选项，permutations内存够大CPU够强，就选1000，这里选100。越大越好，越大结果越精确。

phenotype二选一均可。

这里选false，因为我们的基因名已经是gene symbol了，不需要转换。

然后选phenotype，最后一行不选，如果基因名是测序平台的ID的话最后一行要选对应的测序平台并将上面的false改成ture

最后点击最下面的run.

成功的提示：

三. 结果解读

打开储存的文件夹，可以看到很多GSEA结果图

值得注意的是一个名为index的网页文件

打开它我们看以看到有446个基因sets在APOA2low中上调，1008个在APOA2high中上调，

点击enrichment results in html，得到下图的列表，GS为基因集的名字，SIZE代表该基因集下的基因总数，ES代表Enrichment score, NES代表归一化后的Enrichment score, NOM p-val代表pvalue，表征富集结果的可信度，FDR q-val代表qvalue, 是多重假设检验矫正后的p值，注意GSEA采用pvalue < 5%, qvalue < 25% 对结果进行过滤。Rank at max表示ES分数达到最大值时最后一个基因所处的位置。说的通俗一点，比如说在我们此次GSEA的5000个基因里，有3个基因对于对于cell cycle通路的ES分数有正贡献，三个基因分别位于第1，2，3位（本文开头说过，GSEA会根据基因表达差异大小对基因进行排序打分），他们对于ES分数的贡献分别为0.1，0.01，0.001，ES分数在第3位达到最大值0.111, 此时的Rank at max就为3,第三位往后的基因对cell cycle通路的ES分数是负贡献。

GSEA的结果图：

峰值的顶点为最大富集分数，下面的黑色竖条表示基因的位置，在最大富集分数对应横坐标位置左侧的黑线表示leading edge subset，意思是这些gene的排序靠前，对ES分数有主要贡献。总而言之，这张图主要是为了展示GSEA结果，说明我们确实通过GSEA分析发现了low和high两个表型在这条通路存在差异，这张图上的其他参数不必细究。

补充一下：

1. GSEA软件自带的geneset只有人的，如果是鼠的测序数据则要用鼠的geneset，把它下载下来再导入GSEA，很遗憾，我没有找到鼠的geneset。或者将鼠的基因转换为人的进行GSEA分析，这种转换方法我后续会进行介绍。

2. GSEA软件还有很多可供调整的参数，文中没有提到，还有一些其他功能，如leading edge analysis，这些功能的介绍我放在了文末的云盘中。

最后祝大家都能跑出结果！

资料链接：

数据集（dataset)以及可直接用来跑的的.cls与.gct文件，不过我还是建议用dataset从头到尾自己分析一遍：

链接：https://pan.baidu.com/s/1mTlXDXSUe9CRMZKu8P74uw

提取码：oyl2

java64位：

链接：https://pan.baidu.com/s/1qXoT_ocDA1F8zWlxbIJpRg

提取码：i1xz

参考资料与深度学习资料：

链接：https://pan.baidu.com/s/1jNc3Ho8lRoONuKyCDbgJHw

提取码：lvds