基于Weka平台的分类算法在启动子识别中的应用

蒋璐凯

摘 要： 在“后基因组”时代，对于DNA功能元件的注释，尤其是启动子这类关键的调控元件的鉴定是进一步理解人类基因组繁杂调控网络的重要研究内容。本文基于高通量测序数据对细胞系H1-hesc中的基因启动子进行识别分类，利用数据挖掘软件Weka基于启动子组蛋白修饰特征建立分类模型，比较各分类算法性能优劣，以期应用最佳分类器在其它细胞系中识别分类启动子。

关键词： 启动子；高通量测序；Weka；分类算法

Abstract：The definition of DNA functional elements （especially promoters） is an important research topic in understanding the regulatory network of human genome. This paper identifies types of gene promoters in H1-hesc cell line based on the high-throughput data and then builds classifiers between different types of promoters according to the data of histone modification feature. Finally the paper compares four classifiers' performance and it is expected to apply the best model in the predicting and identifying promoters in other cell lines.

Key words： promoter；high-throughput sequencing；Weka；classification algorithms

引言

Weka是一款基于Java开发环境的机器学习软件，其全称是怀卡拓知识分析环境（Waikato Environment for Knowledge Analysis）。这款开源的数据挖掘软件发展到现在已由最初应用于农业领域而扩展到更多不同领域，尤其是以教育和研究为主的技术科研范畴中。其优势在于：在GNU（General Public License）准则下免费使用，几乎可以在任何一个现代计算平台上运行，集数据预处理和预测模型建立功能于一身，还有易用的图形交互界面[1]。Weka可以实现多样的数据挖掘任务，具体包括：数据预处理、聚类、分类、回归分析、结果可视化以及特征提取。

本文基于Weka软件平台的分类算法，展开启动子类型识别相关研究。随着对于人类基因组的研究进入“后基因组时代”，基因组学的研究重心已经由揭示生命体的遗传信息和密码转移到对分子整体水平的功能研究上来[2]。作为控制基因转录调控起始的关键DNA元件—启动子，是基因组学的研究热点，对于其类型识别预测等层面的深入研究有助于理解基因的表达调控机制等生物学特性，为疾病诊治增加了新方法，为进一步构建生物表达调控网络提供基础[3]。

1 基于高通量测序技术的启动子分类

随着生物信息学领域高通量测序技术的发展，新一代测序技术为启动子识别引入了新的数据支持，极大程度上促进了启动子区域的定位和启动子功能的定性。启动子存在于基因的转录起始位点附近，一般是上游区域（靠近5端），是一段能够引导特异性基因表达活动的DNA序列[4]。启动子作为一个特殊的调控元件，在其区域会有潜在的RNA聚合酶在DNA上的初始结合位点以及特异性的组蛋白修饰信号，根据全基因组分析的结果表明，包括组蛋白H3第4位赖氨酸（H3K4）甲基化和组蛋白H3第9位赖氨酸乙酰化（H3K9ac）在内的若干组蛋白修饰都会在启动子区域富集[5]。本文基于以上启动子区域特点信息，利用RNA-seq数据和ChIP-Seq数据进行启动子类型的识别及获取组蛋白修饰特征数据。

1.1 数据获取及预处理

本文首先从UCSC基因组浏览器上获取人类基因组g19版本的注释基因数据，其主要包含信息见表1。研究时，对注释基因可根据以下条件进行预处理：转录起始位点唯一且转录起始位点上下游各10 kbp的区域内不包含其它基因任何位点的基因，最后得到7 732个符合条件的基因。然后依然从UCSC中下载细胞系H1-hesc的2个全细胞RNA-seq测序数据文件以及该细胞系的RNA聚合酶II的ChIP-Seq数据（版本号为wgEncodeEH000563）。最后，从基因表达综合数据库（Gene Expresion Omnibus，GEO）中下载细胞系H1-hesc的6种组蛋白修饰（与活跃启动子相关的H3K9ac、H3K27ac 和H3K4me1/2/3以及与非活跃启动子相关的H3K27me3）数据。由于从GEO中直接下载的组蛋白修饰数据的BED文件都是比对到人类基因组g18的，而本文其余的數据都是基于g19的，因此这里需要对组蛋白修饰数据利用UCSC的LiftOver工具设计进行不同版本之间的基因组坐标转换。

1.2 启动子分类

启动子是DNA调控元件，是基因转录活动“开关”。启动子是否具有生物活性，可以根据基因是否出现转录活动，即采用基因的表达水平进行衡量。为此，本文利用RNA-seq数据计算7 732个注释基因在细胞系H1-hesc的表达情况，衡量指标为RPKM（Reads Per Kilobase per Million mapped reads），其计算公式如下：

RPKM是每百万reads中来自某个基因每一千碱基区域上的reads数量，能够有效地反映基因真实的表达水平[6]。由于有2个数据文件，将2个计算结果取平均值作为基因的RPKM值。为了更好地拟合真实情况，尽可能减小避免测序误差带来的影响，在此人为规定RPKM值大于0.1以上的基因为表达基因，根据计算结果细胞系H1-hesc中61%的基因是表达基因。

基因具有表達水平，说明存在着具有生物活性的启动子引导了基因的转录活动。而启动子能够调控起始基因转录，需要结合特异性的RNA聚合酶II，因此利用其ChIP-Seq数据去识别具有RNA聚合酶II富集的候选启动子区域。在此，将具有RNA聚合酶II信号的表达基因的启动子分类为活跃启动子，将具有RNA聚合酶II信号的、但基因RPKM值介于0～0.1之间的启动子分类为弱启动子，将具有RNA聚合酶II信号的、但基因RPKM值为0的启动子分类为预备启动子。最终，分类结果如图1所示。其中，活跃启动子1 260个，弱启动子705个以及预备启动子81个。

1.3 启动子组蛋白修饰特征

组蛋白修饰会在启动子区域富集，具有显著的局部的峰和广泛的分布，而在不同类型的启动子中各个组蛋白修饰特征分布又会存在一定的差异性，因此本文研究细胞系H1-hesc的6个组蛋白修饰数据在启动子区域的分布情况。一般认为，基因启动子主要是在转录起始位点上游1 kbp范围内。为此，可将基因转录起始位点上下游各1 kbp的区域作为候选启动子区域，进行组蛋白修饰信号特征的提取。将2 kbp区域划分为10个连续且不重叠的、长度为200 bp的小bins，然后统计每个基因bins上的各个组蛋白修饰read的分布情况。每个组蛋白修饰数据均有2个实验数据文件，为此取二者统计结果的平均值作为组蛋白修饰read落于某个bins内的数目。研究可得组蛋白修饰分布情况如图2所示。

2 分类算法及性能比较

基于先前的工作，已经得到了各个类型启动子的组蛋白修饰特征数据，每个启动子是60维的特征向量，需要处理的是一个三分类问题。对于有些分类算法，如支持向量机在设计时针对的是二分类问题。为此本文采取的方法是一对一策略，即在每两类之间建立分类器，那么三分类问题中会建立3个分类器，对于新的未知样例将根据3个分类器的投票结果来判定其类别[7]。Weka提供了多分类的分析环境，在Classify目录下选择meta中的MultiClassClassifier，然后根据实验对象选择合适的分类算法就可以实现多分类。本文选择了4种分类算法进行比较，分别是：基于C4.5决策树学习算法的J48、随机森林（Random Forest）、基于径向基核函数的LibSVM以及朴素贝叶斯网络。采取10折交叉验证的方法进行分类器评估及选择。分类器相关参数都是默认值。研究中，各分类算法的性能比较结果可见表2。

从表2的结果对比中，综合各个指标可以看出：随机森林分类算法在启动子识别分类中的性能较为优异。在一般分类预测问题中，随机森林可以胜任预测类问题，尤其是多分类问题的第一选择。图3即是随机森林算法在该分类预测中结果的混淆矩阵及分类器的ROC曲线（曲线1为基于活跃启动子、曲线2基于弱启动子、曲线3为基于预备启动子）。可以进一步看出，分类器对于3个类型的启动子预测准确率都在70%以上，这在多分类不平衡问题中是一个较好的结果。因此，可以应用这一经过训练的分类器在其他细胞系中去识别预测启动子类型。

3 结束语

本文主要研究了基于Weka数据挖掘平台的分类算法在启动子识别分类中的应用。基于第二代测序技术的实验数据对细胞系H1-hesc中的基因启动子进行了识别分类并提取组蛋白修饰特征，然后对比Weka中的4种分类算法在启动子分类预测上的性能优劣，得到随机森林分类算法能较好对启动子进行分类预测，今后的相关研究工作将进一步优化分类器，从而提高模型性能。

参考文献

[1] WITTEN I H FRANK E HALL M A. Data mining： Practical machine learning tools and techniques[M]. 3rd ed. 李川，张永辉，译. 北京：机械工业出版社，2014.

[2] GIFFORD C A ZILLER M J GU Hongcang et al. Transcriptional and epigenetic dynamics during specification of human embryonic stem cells[J].Cell 2013 153（5）：1149-1163.

[3] RAZIN S V GAVRILOV A A ULYANOV S V. Transcription-controlling regulatory elements of the eukaryotic genome[J]. Molecular Biology 2015 49（2）：185-194.

[4] Davari K Lichti J Gallus C et al. Rapid genome-wide recruitment of RNA polymerase II drives transcription splicing and translation events during T cell responses[J]. Cell Reports 2017 19（3）：643-654.

[5] BARSKI A CUDDAPAH S CUI K et al. High-resolution profiling of histone methylations in the human genome[J]. Cell 2007 129（4）：823-837.

[6] MORTAZAVI A WILLIAMS B A MCCUE K et al. Mapping and quantifying mammalian transcriptomes by RNA-seq[J]. Nature Methods 2008 5（7）：621-628.

[7] XU Wenxuan ZHANG Li. Human promoter recognition based on single nucleotide statistics and support vector machine ensemble[J]. Journal of Computer Applications 2015 35（10）：2808-2812.