蔣璐凱

摘 要： 在“后基因組”時(shí)代，對(duì)于DNA功能元件的注釋，尤其是啟動(dòng)子這類關(guān)鍵的調(diào)控元件的鑒定是進(jìn)一步理解人類基因組繁雜調(diào)控網(wǎng)絡(luò)的重要研究內(nèi)容。本文基于高通量測(cè)序數(shù)據(jù)對(duì)細(xì)胞系H1-hesc中的基因啟動(dòng)子進(jìn)行識(shí)別分類，利用數(shù)據(jù)挖掘軟件Weka基于啟動(dòng)子組蛋白修飾特征建立分類模型，比較各分類算法性能優(yōu)劣，以期應(yīng)用最佳分類器在其它細(xì)胞系中識(shí)別分類啟動(dòng)子。

關(guān)鍵詞： 啟動(dòng)子；高通量測(cè)序；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開發(fā)環(huán)境的機(jī)器學(xué)習(xí)軟件，其全稱是懷卡拓知識(shí)分析環(huán)境（Waikato Environment for Knowledge Analysis）。這款開源的數(shù)據(jù)挖掘軟件發(fā)展到現(xiàn)在已由最初應(yīng)用于農(nóng)業(yè)領(lǐng)域而擴(kuò)展到更多不同領(lǐng)域，尤其是以教育和研究為主的技術(shù)科研范疇中。其優(yōu)勢(shì)在于：在GNU（General Public License）準(zhǔn)則下免費(fèi)使用，幾乎可以在任何一個(gè)現(xiàn)代計(jì)算平臺(tái)上運(yùn)行，集數(shù)據(jù)預(yù)處理和預(yù)測(cè)模型建立功能于一身，還有易用的圖形交互界面[1]。Weka可以實(shí)現(xiàn)多樣的數(shù)據(jù)挖掘任務(wù)，具體包括：數(shù)據(jù)預(yù)處理、聚類、分類、回歸分析、結(jié)果可視化以及特征提取。

本文基于Weka軟件平臺(tái)的分類算法，展開啟動(dòng)子類型識(shí)別相關(guān)研究。隨著對(duì)于人類基因組的研究進(jìn)入“后基因組時(shí)代”，基因組學(xué)的研究重心已經(jīng)由揭示生命體的遺傳信息和密碼轉(zhuǎn)移到對(duì)分子整體水平的功能研究上來[2]。作為控制基因轉(zhuǎn)錄調(diào)控起始的關(guān)鍵DNA元件—啟動(dòng)子，是基因組學(xué)的研究熱點(diǎn)，對(duì)于其類型識(shí)別預(yù)測(cè)等層面的深入研究有助于理解基因的表達(dá)調(diào)控機(jī)制等生物學(xué)特性，為疾病診治增加了新方法，為進(jìn)一步構(gòu)建生物表達(dá)調(diào)控網(wǎng)絡(luò)提供基礎(chǔ)[3]。

1 基于高通量測(cè)序技術(shù)的啟動(dòng)子分類

隨著生物信息學(xué)領(lǐng)域高通量測(cè)序技術(shù)的發(fā)展，新一代測(cè)序技術(shù)為啟動(dòng)子識(shí)別引入了新的數(shù)據(jù)支持，極大程度上促進(jìn)了啟動(dòng)子區(qū)域的定位和啟動(dòng)子功能的定性。啟動(dòng)子存在于基因的轉(zhuǎn)錄起始位點(diǎn)附近，一般是上游區(qū)域（靠近5端），是一段能夠引導(dǎo)特異性基因表達(dá)活動(dòng)的DNA序列[4]。啟動(dòng)子作為一個(gè)特殊的調(diào)控元件，在其區(qū)域會(huì)有潛在的RNA聚合酶在DNA上的初始結(jié)合位點(diǎn)以及特異性的組蛋白修飾信號(hào)，根據(jù)全基因組分析的結(jié)果表明，包括組蛋白H3第4位賴氨酸（H3K4）甲基化和組蛋白H3第9位賴氨酸乙酰化（H3K9ac）在內(nèi)的若干組蛋白修飾都會(huì)在啟動(dòng)子區(qū)域富集[5]。本文基于以上啟動(dòng)子區(qū)域特點(diǎn)信息，利用RNA-seq數(shù)據(jù)和ChIP-Seq數(shù)據(jù)進(jìn)行啟動(dòng)子類型的識(shí)別及獲取組蛋白修飾特征數(shù)據(jù)。

1.1 數(shù)據(jù)獲取及預(yù)處理

本文首先從UCSC基因組瀏覽器上獲取人類基因組g19版本的注釋基因數(shù)據(jù)，其主要包含信息見表1。研究時(shí)，對(duì)注釋基因可根據(jù)以下條件進(jìn)行預(yù)處理：轉(zhuǎn)錄起始位點(diǎn)唯一且轉(zhuǎn)錄起始位點(diǎn)上下游各10 kbp的區(qū)域內(nèi)不包含其它基因任何位點(diǎn)的基因，最后得到7 732個(gè)符合條件的基因。然后依然從UCSC中下載細(xì)胞系H1-hesc的2個(gè)全細(xì)胞RNA-seq測(cè)序數(shù)據(jù)文件以及該細(xì)胞系的RNA聚合酶II的ChIP-Seq數(shù)據(jù)（版本號(hào)為wgEncodeEH000563）。最后，從基因表達(dá)綜合數(shù)據(jù)庫（Gene Expresion Omnibus，GEO）中下載細(xì)胞系H1-hesc的6種組蛋白修飾（與活躍啟動(dòng)子相關(guān)的H3K9ac、H3K27ac 和H3K4me1/2/3以及與非活躍啟動(dòng)子相關(guān)的H3K27me3）數(shù)據(jù)。由于從GEO中直接下載的組蛋白修飾數(shù)據(jù)的BED文件都是比對(duì)到人類基因組g18的，而本文其余的數(shù)據(jù)都是基于g19的，因此這里需要對(duì)組蛋白修飾數(shù)據(jù)利用UCSC的LiftOver工具設(shè)計(jì)進(jìn)行不同版本之間的基因組坐標(biāo)轉(zhuǎn)換。

1.2 啟動(dòng)子分類

啟動(dòng)子是DNA調(diào)控元件，是基因轉(zhuǎn)錄活動(dòng)“開關(guān)”。啟動(dòng)子是否具有生物活性，可以根據(jù)基因是否出現(xiàn)轉(zhuǎn)錄活動(dòng)，即采用基因的表達(dá)水平進(jìn)行衡量。為此，本文利用RNA-seq數(shù)據(jù)計(jì)算7 732個(gè)注釋基因在細(xì)胞系H1-hesc的表達(dá)情況，衡量指標(biāo)為RPKM（Reads Per Kilobase per Million mapped reads），其計(jì)算公式如下：

RPKM是每百萬reads中來自某個(gè)基因每一千堿基區(qū)域上的reads數(shù)量，能夠有效地反映基因真實(shí)的表達(dá)水平[6]。由于有2個(gè)數(shù)據(jù)文件，將2個(gè)計(jì)算結(jié)果取平均值作為基因的RPKM值。為了更好地?cái)M合真實(shí)情況，盡可能減小避免測(cè)序誤差帶來的影響，在此人為規(guī)定RPKM值大于0.1以上的基因?yàn)楸磉_(dá)基因，根據(jù)計(jì)算結(jié)果細(xì)胞系H1-hesc中61%的基因是表達(dá)基因。

基因具有表達(dá)水平，說明存在著具有生物活性的啟動(dòng)子引導(dǎo)了基因的轉(zhuǎn)錄活動(dòng)。而啟動(dòng)子能夠調(diào)控起始基因轉(zhuǎn)錄，需要結(jié)合特異性的RNA聚合酶II，因此利用其ChIP-Seq數(shù)據(jù)去識(shí)別具有RNA聚合酶II富集的候選啟動(dòng)子區(qū)域。在此，將具有RNA聚合酶II信號(hào)的表達(dá)基因的啟動(dòng)子分類為活躍啟動(dòng)子，將具有RNA聚合酶II信號(hào)的、但基因RPKM值介于0～0.1之間的啟動(dòng)子分類為弱啟動(dòng)子，將具有RNA聚合酶II信號(hào)的、但基因RPKM值為0的啟動(dòng)子分類為預(yù)備啟動(dòng)子。最終，分類結(jié)果如圖1所示。其中，活躍啟動(dòng)子1 260個(gè)，弱啟動(dòng)子705個(gè)以及預(yù)備啟動(dòng)子81個(gè)。

1.3 啟動(dòng)子組蛋白修飾特征

組蛋白修飾會(huì)在啟動(dòng)子區(qū)域富集，具有顯著的局部的峰和廣泛的分布，而在不同類型的啟動(dòng)子中各個(gè)組蛋白修飾特征分布又會(huì)存在一定的差異性，因此本文研究細(xì)胞系H1-hesc的6個(gè)組蛋白修飾數(shù)據(jù)在啟動(dòng)子區(qū)域的分布情況。一般認(rèn)為，基因啟動(dòng)子主要是在轉(zhuǎn)錄起始位點(diǎn)上游1 kbp范圍內(nèi)。為此，可將基因轉(zhuǎn)錄起始位點(diǎn)上下游各1 kbp的區(qū)域作為候選啟動(dòng)子區(qū)域，進(jìn)行組蛋白修飾信號(hào)特征的提取。將2 kbp區(qū)域劃分為10個(gè)連續(xù)且不重疊的、長度為200 bp的小bins，然后統(tǒng)計(jì)每個(gè)基因bins上的各個(gè)組蛋白修飾read的分布情況。每個(gè)組蛋白修飾數(shù)據(jù)均有2個(gè)實(shí)驗(yàn)數(shù)據(jù)文件，為此取二者統(tǒng)計(jì)結(jié)果的平均值作為組蛋白修飾read落于某個(gè)bins內(nèi)的數(shù)目。研究可得組蛋白修飾分布情況如圖2所示。

2 分類算法及性能比較

基于先前的工作，已經(jīng)得到了各個(gè)類型啟動(dòng)子的組蛋白修飾特征數(shù)據(jù)，每個(gè)啟動(dòng)子是60維的特征向量，需要處理的是一個(gè)三分類問題。對(duì)于有些分類算法，如支持向量機(jī)在設(shè)計(jì)時(shí)針對(duì)的是二分類問題。為此本文采取的方法是一對(duì)一策略，即在每兩類之間建立分類器，那么三分類問題中會(huì)建立3個(gè)分類器，對(duì)于新的未知樣例將根據(jù)3個(gè)分類器的投票結(jié)果來判定其類別[7]。Weka提供了多分類的分析環(huán)境，在Classify目錄下選擇meta中的MultiClassClassifier，然后根據(jù)實(shí)驗(yàn)對(duì)象選擇合適的分類算法就可以實(shí)現(xiàn)多分類。本文選擇了4種分類算法進(jìn)行比較，分別是：基于C4.5決策樹學(xué)習(xí)算法的J48、隨機(jī)森林（Random Forest）、基于徑向基核函數(shù)的LibSVM以及樸素貝葉斯網(wǎng)絡(luò)。采取10折交叉驗(yàn)證的方法進(jìn)行分類器評(píng)估及選擇。分類器相關(guān)參數(shù)都是默認(rèn)值。研究中，各分類算法的性能比較結(jié)果可見表2。

從表2的結(jié)果對(duì)比中，綜合各個(gè)指標(biāo)可以看出：隨機(jī)森林分類算法在啟動(dòng)子識(shí)別分類中的性能較為優(yōu)異。在一般分類預(yù)測(cè)問題中，隨機(jī)森林可以勝任預(yù)測(cè)類問題，尤其是多分類問題的第一選擇。圖3即是隨機(jī)森林算法在該分類預(yù)測(cè)中結(jié)果的混淆矩陣及分類器的ROC曲線（曲線1為基于活躍啟動(dòng)子、曲線2基于弱啟動(dòng)子、曲線3為基于預(yù)備啟動(dòng)子）。可以進(jìn)一步看出，分類器對(duì)于3個(gè)類型的啟動(dòng)子預(yù)測(cè)準(zhǔn)確率都在70%以上，這在多分類不平衡問題中是一個(gè)較好的結(jié)果。因此，可以應(yīng)用這一經(jīng)過訓(xùn)練的分類器在其他細(xì)胞系中去識(shí)別預(yù)測(cè)啟動(dòng)子類型。

3 結(jié)束語

本文主要研究了基于Weka數(shù)據(jù)挖掘平臺(tái)的分類算法在啟動(dòng)子識(shí)別分類中的應(yīng)用。基于第二代測(cè)序技術(shù)的實(shí)驗(yàn)數(shù)據(jù)對(duì)細(xì)胞系H1-hesc中的基因啟動(dòng)子進(jìn)行了識(shí)別分類并提取組蛋白修飾特征，然后對(duì)比Weka中的4種分類算法在啟動(dòng)子分類預(yù)測(cè)上的性能優(yōu)劣，得到隨機(jī)森林分類算法能較好對(duì)啟動(dòng)子進(jìn)行分類預(yù)測(cè)，今后的相關(guān)研究工作將進(jìn)一步優(yōu)化分類器，從而提高模型性能。

參考文獻(xiàn)

[1] WITTEN I H FRANK E HALL M A. Data mining： Practical machine learning tools and techniques[M]. 3rd ed. 李川，張永輝，譯. 北京：機(jī)械工業(yè)出版社，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.