GASE相关知识（Linux）

（Gene Set Enrichment Analysis，GSEA）

GSEA概念：

GSEA：(基因集富集分析)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的贡献。其输入数据包含两部分，一是已知功能的基因集 (可以是GO注释、MsigDB的注释或其它符合格式的基因集定义)，一是表达矩阵，软件会对基因根据其于表型的关联度(可以理解为表达值的变化)从大到小排序，然后判断基因集内每条注释下的基因是否富集于表型相关度排序后基因表的上部或下部，从而判断此基因 集内基因的协同变化对表型变化的影响。

• 个人目前简单理解为: 已知的基因集分别在两组或者多组材料中的表达富集情况。

• GSEA分析似乎与GO分析类似但又有所不同。GO分析更加依赖差异基因，实则是对一部分基因的分析 (忽略差异不显著的基因)，而GSEA是从全体基因的表达矩阵中找出具有协同差异 (concordant differences)的基因集，故能兼顾差异较小的基因。因此二者的应用场景略有区别。另外GO富集是定性的分析，GSEA考虑到了表达或其它度量水平的值的影响。所以，对于时间序列数据或样品有定量属性时，GSEA的优势会更明显，不需要每个分组分别进行富集，直接对整体进行处理。可以类比于WGCNA分析。

• GO富集分析是先筛选差异基因，再判断差异基因在哪些注释的通路存在富集；这涉及到阈值的设定，存在一定主观性并且只能用于表达变化较大的基因，即我们定义的显著差异基因。而GSEA则不局限于差异基因，从基因集的富集角度出发，理论上更容易囊括细微但协调性的变化对生物通路的影响，尤其是差异倍数不太大的基因集。

GSEA分析原理：

给定一个已经排序的基因列表L和一个预先定义的基因集合（比如某个通路的所有基因）基因集S ，计算S里面的成员s在L里面是随机分布还是主要聚集在L的顶部（排序的前面）或底部（排序的后面），如下图。这些基因排序的依据是可以是不同表型状态下的表达差异也可以是某个自己定义的排秩关系，若研究的基因集S的成员显著聚集在L的顶部或底部，则说明此基因集成员对表型的差异有贡献，也是我们关注的基因集。

原理图

GSEA中几个关键概念：

1、计算富集得分 (ES, enrichment score)： ES反应基因集成员s在排序列表L的两端富集的程度。计算方式是，从基因集L的第一个基因开始，计算一个累计统计值。当遇到一个落在s里面的基因，则增加统计值。遇到一个不在s里面的基因，则降低统计值。每一步统计值增加或减少的幅度与基因的表达变化程度（更严格的是与基因和表型的关联度）是相关的。富集得分ES最后定义为最大的峰值。正值ES表示基因集在列表的顶部富集，负值ES表示基因集在列表的底部富集。
2、评估富集得分(ES)的显著性。通过基于表型而不改变基因之间关系的排列检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。若样品量少，也可基于基因集做排列检验 (permutation test)，计算p-value。
3、多重假设检验矫正。首先对每个基因子集s计算得到的ES根据基因集的大小进行标准化得到Normalized Enrichment Score (NES)。随后针对NES计算假阳性率。（计算NES也有另外一种方法，是计算出的ES除以排列检验得到的所有ES的平均值）
4、Leading-edge subset，对富集得分贡献最大的基因成员。

GSEA分析的目的：

• GSEA分析的目的是要判断S集基因（基于先验知识的基因注释信息，某个关注的基因集合）中的基因是随机分布还是聚集在排序好的L基因集的顶部或底部（这便是富集分析）。

GSEA的安装：

GSEA和Cytoscape的安装都需要java环境，所以首先要搭建java环境。我使用的是conda建立的虚拟环境。

conda activate
conda create -n java    #建一个java环境，将java环境与外部环境隔绝起来
sudo dpkg -l | grep openjdk   可以查看可供安装的java版本741sudo apt-get install openjdk-8-jre   #安装java
java -version     #可以查看java的版本

java环境构建好之后就可以安装GSEA，去官网下载GSEA软件包，我下载的软件包版本如下：

GSEA
然后解压缩，运行gsea.sh文件既可以完成GSEA软件的安装。

GSEA分析自定义功能注释集：

我的数据是小麦转录组数据，以此为例，待分析的数据是小麦的两组实验 (mut_vs_wt)，每组各2个生物重复，故一共4个样本。各个样本有唯一的名称：CS1、CS2、Dasypyrum1、Dasypyrum2。
1、准备数据：
需要三个输入文件分别是：基因表达矩阵.gct文件；样品分组信息.cls文件；功能基因集.gmt文件（gene sets）

a 、基因表达文件：

基因表达文件不用差异基因表达文件。可以直接使用Excel文件生成(测序公司一般都会提供这个文件)，也可以自己构建基因表达数据。例如使用DESeq2里面的counts()函数对dds对象提取标准化之后的数据，经过简单的处理之后，就可以输入GSEA进行分析了。
数据的前期处理是使用Hisat2对转录组数据进行索引，然后使用HTSeq统计表达量counts，之后使用DESeq2对起进行数据标准化。

NormDats <- counts(dds, mormalized = TRUE)

获得基因表达数据文件后将其整理成GSEA所需要的文件格式，如下图：

基因表达文件

• 第一行是固定写法：#1.2，表示版本号，自己准备文件时照抄就行。
• 两个数分别表示gene NAME的数量和样本数量，72748表示此表格有72748个基因，注意72748不是这个表的总行数，而是从第四行开始数的基因数，并且无重复的基因名称。4表是有4个样本(与样本设置文件中的4必须对应)。
• 第三行NAME和DESCRIPTION为固定写法。DESCRIPTION列下面可以写na或任意字符串，后面的就是基因在不同样本中标准化后的表达数据了。
• 第四行及以后的行中，基因的表达值最好是没有取过对数的。比如TPM、RPKM/FPKM或DESeq2标准化后的值均可。
• 将此文件另存为文本文件(制表符分隔)(*.txt)，然后再将后缀名称改为.gct即可。

b 、样品分组信息文件：

本质上是一个文本文件，需按照以下格式书写：

样品分组信息文件格式

• 第一行4表示共有4个样本；2 表示是两个分组；1 是固定写法，不要改。
• 第二行：以#开始，#后面没有空格，Dasypyrum和CS为2个分组的名字，可以自行更改为自己试验设定的分组名字，名字内部不只允许有字母、数字、下划线。分组信息（有几个分组便写几个，多个分组在比较分析时，后面需要选择待比较的任意2组）。
• 第三行的书写方式需要格外注意，内容是样本对应的组名，这里写的是4个字符串(代表4个样本)，样本重复必须赋予相同的名字，以便GSEA判断样品所属的组，并且其排列的相对顺序必须与稍后要准备的基因表达文件(或基因表达矩阵)中对应的样本书写顺序相同。例如相同处理的不同重复在自己试验记录里如果是Treat6h_1、Treat6h_2、Treat6h_3，但是在这里一定都要写成一样的值Treat6h，名字相同的代表样品属于同一组。写完之后注意需要将文件后缀名保存为.cls。

c 、功能基因集文件：

GSEA官网只提供了人类的数据，但是掌握了官网中基因表达矩阵和注释文件的数据格式，就可以根据自己研究的物种，在公共数据库下载对应物种的注释数据，自己制作格式一致的功能基因集文件，这样便就可以做各种物种的GSEA富集分析了。
基因集数据库文件的格式：

功能基因集文件

• 第一列是某一个生物过程、细胞组分、分子功能或信号通路等等，是一个基因集(gene set)的名称。
• 第二列是该基因集的编号或其它描述信息均可。
• 第三列及以后的每列是从属于该基因集的基因，有的是一个，有的是两个或多个。你也可以输入自己感兴趣的基因，只需满足以上格式即可。
• 列与列之间以Tab键隔开，基因与基因之间也以Tab键隔开，文件后缀名写为.gmt。
  本文分析的对象是小麦，MSigDB并无提供，因此需要自己去Ensemble plant数据库或国际小麦基因组测序联盟(International Wheat Genome Sequencing Consortium,IWGSC)上下载。再将其加工修改为一个.gmt文件。
  我在IWGSC上下载小麦的GO注释数据后整理成下图文件格式，需要对其进行整理去除重复。
  IWGSC上下载小麦的GO注释数据

awk 'BEGIN{OFS=FS="\t"}{anno=$1"\t"$2; a[anno]=a[anno]==""?$3:a[anno]"\t"$3}END{for(i in a) print i,a[i]}' ID_GOanalisiy.txt > GOanalisiy_ID_only.gmt

GSEA功能基因集格式文件

由此，GSEA分析所需要的3个文件已经整理完成。

d 、GSEA分析：

分析图1

• 点击1处的Load data后，点击2处的Browse for files；将已经准备好的样本分组文件(.cls)、基因表达文件(.gct)和基因集数据库文件(.gmt)一并上传。注意：如果分析的是人的基因数据，且调用GSEA数据库中的gmt文件，软件会自动下载，此处可以无需上传本地的gmt文件。若文件格式没问题会弹出一个提示There were no error的框（图中的3），证明文件上传成功；若出错，请仔细核对文件格式。

• 文件上传成功之后，点击4，运行GSEA，会出现如下图的界面。

  
  i分析图2
• Expression dataset： 导入表达数据集文件，点击后自动显示上一步中从本地导入软件内的文件，所以一定要确认上一步导入数据是否成功；
• Gene sets database： 基因功能集数据库，点进去可以选择从本地导入；在联网的情况下软件也可以为自动下载GSEA官网中的gene sets文件；
• 在Collapse dataset to gene symbols处选择No_Collapse，这是因为在基因表达文件(.gct)和基因集数据库文件(.gmt)中已经有对应好的基因名称，不需要再做转换。
• 在Permutation type处选择gene_set。
• 在6处即可点击Run来运行GSEA。
• 7处会显示运行的状态running，当出现success后点点击 (等同于点击Show result folder并打开文件夹中的index.html)，即可弹出分析报告。

如此GSEA的分析就完成了，后面再写对其结果的分析。

【参考】：

一文掌握GSEA富集分析
常见GSEA分析的疑问
GSEA使用实战