在基因组和转录组迅速发展的今天，如何从庞大的基因表达数据库里挖掘出有价值的信息，并做出科学的诠释，是基因分析领域的重要挑战和关键问题。生物统计学家们利用t检验、比值法、固定效应模型、混合效应模型等多种统计分析方法，对筛选出的差异表达的单基因进行分析比较，推测和预测其可能的功能^[1]。然而，单个基因分析在组间基因表达量差别较小、差异基因数量过大且没有统一分类信息时，容易出现假阴性和无法解释的现象。之后，生物学家们总结了单基因分析的缺陷，提出了以基因集(Gene set)为基础的基因富集分析法(Gene set enrichment analysis, GSEA)^[1]。

从Disease Ontology(DO), Gene Ontology(GO), Kyoto encyclopedia of genes and genomes(KEGG)，到Reactome Pathway等多种基因富集性分析方法不断被推出。科学家们借助各种生物信息学的数据库和分析工具进行统计分析, 对不同层次、不同来源的数据进行归纳与整合，力图找到目标基因与库中有相同或类似功能基因集之间的相关性^[2-8]。具体是对一组基因(共表达或者差异表达)中某个功能类别的显著性利用超几何分布分析进行检验，这个功能类别可能就是导致样本性状差异的最重要的功能差别，其功能特征会阐明样本变化的内在生物学意义。常用的统计学方法主要是基于Fisher test^{[2, 3, 7]}和chi-square test^[8]检验。

目前做富集分析已经开发了很多软件工具：有GFINDer, GOALIE, GOTM, L2L, Ontology Traverser, ProfCom, SeqExpress, STEM等。但是绝大多数都是通过web和其对应的数据库中的数据进行分析，对于个人的或者特殊的物种无法分析。其中GSEA是最有名的软件，但是它基于芯片的分析软件，老版本对现今的RNA-seq数据不适用，最新的版本用起来并不是很方便。DAVID可用Gene ID或symbol做GO分析，但是无法对其数据库中没有的基因ID进行分析。本文中介绍的Practical Extraction and Report Language(perl)脚本可以利用研究者现有的基因表达差异和注释数据，进行分析。

1 数据库

基因集是先验的生物学知识的集合，即已被证实了的生物通道、基因共表达信息等。目前常用的定义基因集的基因注释数据库有GO和KEGG等。

GO是基因本体联合会建立的数据库，旨在构建一个可以随着细胞中基因和蛋白质功能知识的累积和变化，进行统一描述的语义词汇标准。它是根据不同实验目的，筛选出差异表达基因，分析差异表达基因在GO数据库中的分布情况，从而诠释实验中样本差异在基因功能层面上的差异体现。GO共有3层结构的定义方式：分子功能，生物过程和细胞组成^[9]。在GenBank中，基因与其编码的蛋白都有着特殊的ID编号一一对应，也有能通过序列注释的方法找到与之对应的GO号^[10]，找到检索的全部分子功能所处的位置和关系信息。

KEGG是整合了生物化学、基因组学以及系统功能组学的信息的基因功能系统分析知识库，这个数据库有助于科研工作者从分子通路水平来解释生物系统高层次功能的数据库。KEGG不仅可以提供出色的整合代谢途径查询，还可以利用Java的图形工具访问基因组图谱，对基因组图谱、表达图谱和序列进行图形比较和通路计算。KEGG是一个综合数据库，它们大致分为系统信息、基因组信息和化学信息三大类。进一步可细分为16个主要的数据库^[11]。科学家们将属于同一KEGG基因通路(KO分类)的基因定义为一个基因集。在生物体内，不同的生物学功能是通过不同基因间的相互协调互作而实现的，KEGG能够通过比较通路富集的情况，推测差异表达基因可能参与的主要生化代谢和信号转导路径，从而推测其可能的路径数目及功能^[12]，这对寻找出不同样品的差异表达基因可能与哪些生物学通路的改变有关具有提示作用。

2 方法 2.1 在windows操作系统下的安装Perl

Perl语言是拉里·沃尔(Larry Wall)于1987年12月18日发表的。它借鉴采用了C、awk、sed、shell等其他编程语言的特性以及优点。Perl语言可以在UNIX、Linux、MAC和Windows等操作系统下跨系统平台运行，并由第三方代码库CPAN^[13](Comprehensive Perl Archive Network)不断更新和维护，它包含了极丰富perl写成的软件包和其源文件^[14]。

在Enrich_analysis.pl脚本运行前，电脑上需事先安装一个可执行perl语言编辑与运行的解释器。一般的UNIX(含Linux与Mac OS)系统都含有perl解释器，因此在终端输入“perl-v”命令，就可以查看perl的版本。Windows系统一般可与“ActivePerl”和“Strawberry Perl”兼容。ActivePerl^[15]是ActiveState公司开发的一个Perl语言运行环境，是一套面向各个平台、集成度高、易于使用的Perl脚本解释器。ActivePerl目前最新版本为5.26，可在http://www.activestate.com/activeperl/downloads页面下下载，其中“msi”文件是安装程序。安装完成之后，打开cmd窗口输入“perl-v”如有返回perl版本信息则表示安装成功^[16]。

2.2 数据的格式

输入两个文件：第1个文件是各个基因对应的GO分类或KEGG分类; 第2个文件是差异表达基因的名称。

Enrich_analysis.pl要求输入两个数据文件。各物种GO注释数据可以通过Blast2GO获得，或者从http://geneontology.org/page/download-annotations与http://bioinfo.cau.edu.cn/agriGO/download.php网页下载(KEGG数据目前收费)，之后用Excel打开，删除掉多余信息，只留下两行，第一行为基因编号或基因名，第二列是对应的注释数据。数据结构见表 1。

差异基因数据文件只有一列：差异基因的基因编号或基因名。

将Excel表另存为“文本文件(用制表符“Tab”分割)(*.txt)”类型的文本文件，例如annotation.txt和DEGlist.txt。将这两个数据文件与Enrich_analysis.pl存放在同一个目录下，如C盘根目录下。

2.3 Enrich_analysis.pl的运行

通过点击任务栏上的“开始-运行”，可以运行存放在数据目录下的Enrich_analysis.pl文件。在运行框中输入“cmd”，打开“cmd命令提示符”，继而在窗口中键入“cd C:\\”后回车，转到含perl命令执行程序的目录下。在“C: >”提示符下输入“perl Enrich_analysis.pl-h”回车后输出帮助信息。在“C: >”提示符下输入“perl Enrich_analysis.pl annotation.txt DEGlist.txt > Result.txt”，回车开始运行。输出的结果Enrich_analysis.pl会在“C:\\”目录下自动生成文件“Result.txt”，并存放入其中。“Result.txt”分8列，依次表示：Enrich term、Ontology、Description、Number in input、Number in BG/Ref、p-value、Genes、FDR。

3 结果与分析 3.1 Enrich_analysis.pl的结果与已发表sci文章比较

为了证明程序的可靠性，我们利用Priest et al., 2015(Plos One)^[17]，Hu, et al., 2016(Bmc Plant Biology)^[18]，Wu, et al., 2016(Plant Physiology and Biochemistry)^[19]三篇sci文章的数据进行了比较分析(见表 2、表 3、表 4)。可以发现结果基本是相同的，但是p-value略有差异。这是由于我们使用的GO注释数据与作者使用的不完全相同。更多的比较可通过在线工具agriGO(http://bioinfo.cau.edu.cn/agriGO)^[20]进行特定物种的分析。

3.2 Enrich_analysis.pl的原理和特点

做富集分析需要两个条件：1)针对差异表达的基因(比如说基因芯片或者转录组的数据); 2)需要该物种的基因组的所有基因组注释结果作为背景。

各种软件做GO富集性分析时一般都是使用超几何分布^[21](Hypergeometric Distribution)进行计算。而Fisher's Exact Test^[22]是依据超几何概率分布得到。

在富集分析中把基因类比成非黑即白的抽球问题，用来判断和某个注释/分类是否有相关性。也可以列一个2×2的表(见表 5)，进行独立性分析。之后有很多方法可以做检验，经典的有卡方检验和fisher's exact test。但是卡方检验通常也只能做为近似估计值，特别是当样本量比较小时，计算并不准确。而fisher's exact test和超几何模式计算的p-值是一样的。通过FDR(false discovery rate，错误发现率)校正之计算得到FDR值后，以0.05为阈值，≤0.05的注释/分类为显著富集。

那么np1个差异表达基因中正好有n11个基因属于这个注释的概率为：

$ P\left( {X = n11|n{\rm{pp}}, n1{\rm{p}}, n{\rm{p}}1} \right) = \frac{{\left( \begin{array}{l} n1{\rm{p}}\\ n11 \end{array} \right)\left( \begin{array}{l} n2{\rm{p}}\\ n21 \end{array} \right)}}{{{\rm{ }}\left( \begin{array}{l} n{\rm{pp}}\\ n{\rm{p}}1 \end{array} \right)}}{\rm{ }} $

(1)

Fisher's exact test得分表示np1个差异表达基因中至少有n11个属于这个注释的概率：

$ p = 1 - \sum\limits_{{\rm{ }}i = 0}^{n11 - 1{\rm{ }}} {} \frac{{\left( \begin{array}{l} n1{\rm{p}}\\ i \end{array} \right)\left( \begin{array}{l} n2{\rm{p}}\\ n{\rm{p}}1 - i \end{array} \right)}}{{\left( \begin{array}{l} n{\rm{pp}}\\ n{\rm{p}}1 \end{array} \right)}} $

(2)

“Result.txt”为结果输出文档，Enrich_analysis.pl脚本运行操作简单、快捷，所需设置的参数少。

3.2 Enrich_analysis.pl的源代码

Enrich_analysis.pl源代码见附件。

也可在gitHub平台https://gist.github.com/xukaili/15bf411472092ebcf03b676352ec6a00网址下载获得。

4 结论

1) 基因富集分析的主要任务是筛选出功能富集基因，有效地将基因表达信息与基因注释等信息相结合，从统计上把转录组数据与生物学意义很好地衔接起来，使得研究者们能够更加轻松、合理地解读转录组结果，还能达到降维的目的，使转录组的海量基因信息得到充分利用，但目前对生物基因注释信息的准确度和精确度还有待进一步完善。

2) 本文用Perl语言写了Enrich_analysis.pl，根据基因注释信息进行基因富集分析，并利用Fisher's Exact Test做检验，提高数据的可解释性。

3) 利用3篇SCI文章的差异基因数据，对本perl脚本进行了对比测试，证明本perl脚本算法正确可靠。

4) 目前各物种的注释信息并不完善，仍有很多基因尚没有注释信息，随着研究的不断深入，这些未知的基因也在不断得到注释，使得基因数据库在不断更新。但是目前现有的富集分析数据库中的数据并没有及时更新，或者有些国外富集分析网站不能访问。利用本perl脚本对最新的基因注释数据进行富集分析，相对于使用多年不更新的网络版工具更具有优势。