基于機(jī)器學(xué)習(xí)方法的擬南芥基因組DNA復(fù)制時(shí)間預(yù)測(cè)研究

李 椰，李東維，李昭宏，楊若林

(西北農(nóng)林科技大學(xué) 生命科學(xué)學(xué)院，陜西 楊凌 712100)

真核生物基因組的復(fù)制在時(shí)空上受到嚴(yán)格的調(diào)控，細(xì)胞內(nèi)DNA復(fù)制發(fā)生在細(xì)胞周期S期，基因組不同位點(diǎn)DNA在各自開始的時(shí)間點(diǎn)精確復(fù)制，其間不僅發(fā)生DNA復(fù)制，也涉及染色質(zhì)狀態(tài)、表觀信息、染色質(zhì)空間結(jié)構(gòu)的復(fù)制[1-2]；同時(shí)，該復(fù)制過程也伴隨著其他生物學(xué)過程如核小體重塑[3]。根據(jù)DNA復(fù)制時(shí)間早晚特征，可將基因組分為CTR區(qū)(constant timing regions,CTR)和TTR區(qū)(timing transition regions, TTR)[4]。CTR區(qū)是指由多個(gè)相鄰復(fù)制子構(gòu)成的一段DNA區(qū)域，該區(qū)域內(nèi)的復(fù)制子同步激活，出現(xiàn)在S期早期(早期復(fù)制區(qū)CTR)或晚期(晚期復(fù)制區(qū)CTR)[5]。TTR區(qū)則對(duì)應(yīng)于連接上下游2個(gè)具有不同CTR的DNA區(qū)域[1]。目前，利用Repli-Seq(replication timing sequencing)測(cè)序技術(shù)，可以在全基因組水平對(duì)DNA的復(fù)制時(shí)間進(jìn)行系統(tǒng)研究。該技術(shù)通過5-溴-2-脫氧尿苷(bromouridine triphosphate,BrdU)在體內(nèi)替換胸腺嘧啶，從而對(duì)新合成的DNA鏈進(jìn)行標(biāo)記，根據(jù)不同時(shí)期DNA含量的不同，利用熒光激活細(xì)胞分選技術(shù)(fluorescence-activated cell sorting,FACS)對(duì)不同復(fù)制時(shí)間的DNA分選，用BrdU抗體免疫沉淀BrdU標(biāo)記的DNA鏈，對(duì)該DNA進(jìn)行測(cè)序[6]。通過Repli-Seq技術(shù)對(duì)一些生物的DNA復(fù)制時(shí)間進(jìn)行鑒定分析，發(fā)現(xiàn)早晚期復(fù)制的區(qū)域存在著明顯差異，例如復(fù)制早期的DNA區(qū)域的基因密度較高，具有更高的轉(zhuǎn)錄活性，表達(dá)量更高，且富集激活型的染色質(zhì)標(biāo)記[7]，復(fù)制早晚期也與結(jié)構(gòu)變異類型及突變率存在聯(lián)系[8]。同時(shí)，DNA復(fù)制早晚的異常及改變與人類的癌癥有關(guān)[8]。然而，目前研究人員對(duì)DNA復(fù)制時(shí)間調(diào)控機(jī)制的了解尚不完善。

DNA復(fù)制在不同物種內(nèi)也存在差異。早期的研究表明，異六聚體起始識(shí)別復(fù)合物(origin recognition complex, ORC)能特異性地結(jié)合到酵母基因組的特定位點(diǎn)[9]，這對(duì)于保證酵母基因組的準(zhǔn)確復(fù)制有著重要的作用。然而，在高等真核生物中，僅僅依據(jù)結(jié)合位點(diǎn)的序列特征并不能完全解釋該復(fù)合物與DNA的結(jié)合，這意味著有DNA序列特征之外的其他因素也參與了DNA復(fù)制起始的調(diào)控[10]。近年來，針對(duì)其他物種的研究發(fā)現(xiàn)，基因組局部的染色質(zhì)狀態(tài)、表觀遺傳修飾和染色質(zhì)三維結(jié)構(gòu)變化會(huì)影響DNA復(fù)制的起始時(shí)間。例如，人類甲基化的RNA基因復(fù)制時(shí)間較晚[11]。在小鼠中，激活型組蛋白修飾(H3K4me1、H3K4me2、H3K4me3、K3K20me1、H3K36me3、H3K9ac、H3K27ac)與早期復(fù)制相關(guān)[12]；果蠅中也有相似的發(fā)現(xiàn)[10]。另外，DNA的乙酰化也會(huì)影響復(fù)制時(shí)間[1,13]。在植物中，玉米的DNA復(fù)制與動(dòng)物存在差異[14]。但有關(guān)植物全基因組水平上DNA復(fù)制的時(shí)空調(diào)控機(jī)制目前尚少有報(bào)道。

隨著二代測(cè)序技術(shù)的普及和發(fā)展，植物基因組中尤其是擬南芥基因組中大量的表觀遺傳修飾被測(cè)定[15-25]，為探索DNA的復(fù)制時(shí)間與表觀遺傳修飾之間的關(guān)系提供了有力的數(shù)據(jù)支撐。與此同時(shí)，機(jī)器學(xué)習(xí)方法目前已在圖像分類、語音識(shí)別等領(lǐng)域取得了重大突破和進(jìn)展[26-27]。近年來，研究者將機(jī)器學(xué)習(xí)方法與生物組學(xué)數(shù)據(jù)結(jié)合，更為全面、系統(tǒng)地揭示了DNA序列與表觀遺傳修飾的聯(lián)系[28]，以及表觀遺傳修飾與基因表達(dá)之間的生物學(xué)機(jī)制[29-30]。本研究以模式植物擬南芥為對(duì)象，系統(tǒng)地收集和整合染色質(zhì)開放狀態(tài)(DNase-Seq)數(shù)據(jù)以及表觀遺傳修飾特征(ChIP-Seq)數(shù)據(jù)，利用機(jī)器學(xué)習(xí)方法，構(gòu)建擬南芥基因組DNA復(fù)制時(shí)間的預(yù)測(cè)模型；在此基礎(chǔ)上發(fā)現(xiàn)可能參與調(diào)控DNA復(fù)制時(shí)間的表觀遺傳修飾，以期為后續(xù)進(jìn)一步探究擬南芥DNA復(fù)制時(shí)間的表觀遺傳調(diào)控機(jī)制提供參考。

1 材料與方法

1.1 試驗(yàn)材料

從NCBI數(shù)據(jù)庫收集擬南芥全基因組水平染色質(zhì)開放狀態(tài)的DNase-Seq數(shù)據(jù)和多種表觀遺傳修飾的ChIP-Seq數(shù)據(jù)，相關(guān)信息如表1所示[15-16,18-25,29,31]。擬南芥基因組從Ensembl plant網(wǎng)站上下載(http://plants.ensembl.org/index.html)。

表1 擬南芥表觀遺傳修飾特征和染色質(zhì)狀態(tài)數(shù)據(jù)集信息Table 1 Information on datasets of investigated epigenetic features and chromatin status of Arabidopsis thaliana

擬南芥的全基因組復(fù)制時(shí)間數(shù)據(jù)來自于PRJNA330547。該研究將基因組劃分為長度為1 kb的區(qū)間，通過測(cè)序?qū)⑦@些區(qū)間鑒定和劃分為復(fù)制早期、中期或晚期，并提供了每個(gè)區(qū)間的信號(hào)值[7]。本研究據(jù)此進(jìn)行DNA復(fù)制時(shí)間與各種基因組表觀遺傳修飾信號(hào)之間關(guān)系的分析。

1.2 方 法

1.2.1 DNase-Seq和ChIP-Seq數(shù)據(jù)處理 用fastq-dump(v2.8.2)將下載的擬南芥DNase-Seq和ChIP-Seq原始數(shù)據(jù)(sra格式)轉(zhuǎn)換為fastq格式文件；利用fastqc(v0.11.4)對(duì)reads序列的測(cè)序質(zhì)量進(jìn)行檢查，使用Trimmomatic(v0.36)去除不合格reads兩端低質(zhì)量序列和接頭序列；使用軟件bowtie2將合格的reads序列比對(duì)到擬南芥參考基因組上獲得bam文件，并將不同組織相同類型的表觀修飾數(shù)據(jù)視為不同的特征。最后共得到31個(gè)特征，其中包括29個(gè)表觀遺傳修飾特征和2個(gè)染色質(zhì)開放狀態(tài)特征。

1.2.2 樣本的特征信號(hào)值計(jì)算 對(duì)于每種表觀遺傳修飾特征，利用其對(duì)應(yīng)的bam文件計(jì)算基因組每1 kb區(qū)間的表觀遺傳修飾信號(hào)值。以表觀修飾特征H3K27ac為例，計(jì)算過程為：統(tǒng)計(jì)每個(gè)區(qū)間中的reads數(shù)目；為防止某些區(qū)域沒有reads覆蓋，對(duì)所有區(qū)間的reads數(shù)目均加1；對(duì)所有區(qū)間內(nèi)的reads數(shù)目取對(duì)數(shù)，利用Z-score方法對(duì)壓縮后的reads數(shù)目標(biāo)準(zhǔn)化，得到每個(gè)區(qū)間的信號(hào)值。具體計(jì)算公式如下：

1.2.3 t分布鄰域嵌入(t-distributed stochastic neighbor embedding，t-SNE)分析 t-SNE利用非線性降維算法將高維數(shù)據(jù)投影到低維空間，不僅可以實(shí)現(xiàn)很好的可視化效果，而且可以反映這些特征能否將不同類別的數(shù)據(jù)區(qū)分開。t-SNE的核心思想是數(shù)據(jù)在高維空間內(nèi)的歐氏距離與降維后低維空間內(nèi)的歐式距離能夠最大程度地保持一致。利用KL散度(Kullback-Leibler divergence)衡量高低維歐式距離的差值，通過最小化KL散度可得到最優(yōu)降維。該方法既能聚焦于高維數(shù)據(jù)的局部結(jié)構(gòu)，又能兼顧全局的結(jié)構(gòu)，通過R(v3.6)軟件包Rtsen完成計(jì)算。

1.2.4 34個(gè)特征的相關(guān)性 本研究計(jì)算了34個(gè)特征，包括DNA復(fù)制早期、中期、晚期的信號(hào)和原本31個(gè)表觀遺傳修飾特征兩兩之間的線性關(guān)系。皮爾遜相關(guān)系數(shù)(Pearson correlation coefficient,PCC)能夠很好地度量兩個(gè)隨機(jī)變量的線性關(guān)系，其取值范圍為[-1，1]，其絕對(duì)值的大小反映兩個(gè)變量線性相關(guān)性的強(qiáng)弱，而正負(fù)方向反映變量之間存在正相關(guān)或負(fù)相關(guān)，通過R(v3.6)軟件中的基礎(chǔ)包c(diǎn)or.test函數(shù)，可計(jì)算34個(gè)特征兩兩之間的PCC顯著性。

1.2.5 DNA復(fù)制期早晚預(yù)測(cè)模型以及表觀遺傳修飾特征的重要性評(píng)估 用隨機(jī)森林、多類別邏輯回歸模型和支持向量機(jī)3種機(jī)器學(xué)習(xí)分類器進(jìn)行建模分析，分別由R包e1071、nnet和randomForest執(zhí)行。其中支持向量機(jī)的核函數(shù)選用徑向基函數(shù)(radial basis function)，其余參數(shù)均為默認(rèn)參數(shù)。隨機(jī)森林選擇1 000棵決策樹；隨機(jī)選取80%數(shù)據(jù)作為訓(xùn)練集，剩余20%數(shù)據(jù)作為測(cè)試集用于評(píng)估模型的性能。接受者操作特性曲線(receiver operating characteristic curve，ROC)和ROC曲線下的面積(area under curve,AUC)，以及每個(gè)模型的表觀遺傳修飾特征的重要性，均用對(duì)應(yīng)的R函數(shù)完成計(jì)算。

十折交叉驗(yàn)證：將訓(xùn)練集樣本隨機(jī)分為10份，輪流選取9份構(gòu)建模型，剩余1份用于評(píng)價(jià)模型性能。計(jì)算每個(gè)類別的10次AUC平均值來評(píng)估模型的預(yù)測(cè)性能。

2 結(jié)果與分析

2.1 擬南芥基因組表觀數(shù)據(jù)的可預(yù)測(cè)性及相關(guān)性

為了初步檢查31個(gè)表觀遺傳修飾特征的聯(lián)合信號(hào)能否預(yù)測(cè)對(duì)應(yīng)DNA區(qū)域的復(fù)制時(shí)間早晚，首先使用t-SNE變換將所有擬南芥DNA片段的高維表觀遺傳修飾特征展示到二維空間，然后以3種不同顏色將對(duì)應(yīng)DNA片段的復(fù)制時(shí)間按早、中、晚類別信號(hào)顯示，結(jié)果見圖1。如圖1所示，復(fù)制晚期的DNA片段明顯聚類到一起，與早期和中期片段分開；而早期和中期的DNA片段雖然沒有明顯分開，但也呈現(xiàn)一定程度的各自聚類模式，說明這些表觀遺傳修飾特征整體上能夠有效地用于預(yù)測(cè)DNA區(qū)域的復(fù)制早晚。

31個(gè)表觀遺傳修飾特征與DNA片段復(fù)制早、中、晚3個(gè)時(shí)期信號(hào)值之間的相關(guān)性結(jié)果(圖2)顯示，一些表觀遺傳修飾特征之間存在較高的相關(guān)性，復(fù)制早期的信號(hào)與一些表觀遺傳修飾特征(如H4K16ac)呈較強(qiáng)的正相關(guān)(PCC=0.394,P<2.2×10-16)，復(fù)制晚期的信號(hào)則與之相反(PCC=-0.525,P<2.2×10-16)，而復(fù)制中期的信號(hào)與大部分表觀遺傳修飾特征的相關(guān)性較小。而對(duì)于表觀遺傳修飾特征H2AW而言，它們則與復(fù)制早期信號(hào)呈顯著負(fù)相關(guān)(PCC=-0.508,P<2.2×10-16)，而與復(fù)制晚期的信號(hào)呈正相關(guān)(PCC=0.612,P<2.2×10-16)。

圖1 擬南芥基因組DNA復(fù)制早、中、晚3個(gè)時(shí)期的樣本分群結(jié)果Fig.1 Clustering results of Arabidopsis thaliana genome-wide DNA replication timing classification

每個(gè)小格的扇形面積表示對(duì)應(yīng)兩個(gè)特征的相關(guān)性大小，×表示對(duì)應(yīng)兩個(gè)特征之間相關(guān)性不顯著(P>0.05)The fan-shaped area of each cell represents the degree of correlation between two features,× represents insignificant correlation between two features (P>0.05)圖2 擬南芥基因組表觀遺傳修飾特征與復(fù)制時(shí)間信號(hào)之間的相關(guān)性Fig.2 Correlation among all epigenetic features and signal of DNA replication timing in Arabidopsis thaliana

2.2 DNA復(fù)制時(shí)間的模型預(yù)測(cè)

擬南芥基因組DNA片段按照復(fù)制時(shí)間早晚被分為36 671個(gè)早期樣本、23 113個(gè)中期樣本和9 052個(gè)晚期樣本。將其中一類樣本作為正樣本，其他兩類作為負(fù)樣本進(jìn)行建模。從中各隨機(jī)挑選80%樣本用于模型構(gòu)建和性能的評(píng)價(jià)，剩余20%樣本作為獨(dú)立的測(cè)試集。

為了評(píng)估隨機(jī)森林、多類別邏輯回歸模型和支持向量機(jī)3種分類器的預(yù)測(cè)性能，對(duì)80%數(shù)據(jù)在3種分類器上進(jìn)行十折交叉驗(yàn)證。圖3展示了3種分類器對(duì)復(fù)制早期、中期、晚期三個(gè)類別的平均AUC以及宏平均(Macro)和微平均(Micro)指數(shù)，宏平均和微平均分別表示每個(gè)類別的真陽性平均值和所有類別總真陽性。從圖3可見，3種模型對(duì)DNA復(fù)制早期、中期、晚期預(yù)測(cè)的十折交叉驗(yàn)證平均AUC都在0.8以上，且變異很?。磺?種分類器預(yù)測(cè)性能差異較小，說明復(fù)制時(shí)間早晚與多種表觀遺傳修飾的協(xié)同作用之間可能存在密切的生物學(xué)聯(lián)系。從圖3還可見，DNA復(fù)制晚期樣本平均AUC最高，分類效果最好，說明表觀遺傳修飾特征很大程度上與該區(qū)域的復(fù)制時(shí)間有關(guān)。這與t-SNE分析結(jié)果吻合，較高的Macro和Micro指數(shù)表征了3種分類器關(guān)于DNA復(fù)制時(shí)間穩(wěn)定的綜合預(yù)測(cè)效應(yīng)。

圖3 基于3種不同模型的擬南芥DNA復(fù)制時(shí)間十折交叉驗(yàn)證平均AUCFig.3 Average AUC of 10-fold cross-validation with different classifier of DNA replication timing in Arabidopsis thaliana

為了進(jìn)一步驗(yàn)證模型性能，本研究利用所選取的80%訓(xùn)練集數(shù)據(jù)構(gòu)建了3種模型，再用獨(dú)立的20%測(cè)試集數(shù)據(jù)驗(yàn)證。3種模型的ROC曲線(圖4)和AUC值(表2)均表明，模型對(duì)DNA復(fù)制時(shí)間的預(yù)測(cè)性能良好，而且與十折交叉驗(yàn)證的結(jié)果一致，均表現(xiàn)為復(fù)制晚期的分類效果最佳。3種分類器中以支持向量機(jī)的預(yù)測(cè)效果最佳，其次為隨機(jī)森林和多類別邏輯回歸。由于基于獨(dú)立測(cè)試集構(gòu)建的模型使用了更多樣本，因此，它們的AUC值略高于十折交叉驗(yàn)證的AUC值，使得模型訓(xùn)練更加準(zhǔn)確。

圖4 基于獨(dú)立測(cè)試集的3種擬南芥DNA復(fù)制時(shí)間預(yù)測(cè)模型的ROC曲線Fig.4 ROC of three models in predicting DNA replication timing of Arabidopsis thaliana on an independent test set

表2 基于獨(dú)立測(cè)試集的3種擬南芥DNA復(fù)制時(shí)間預(yù)測(cè)模型的AUC值Table 2 AUC of three models in predicting DNA replication timing on an independent test set

2.3 擬南芥基因組表觀遺傳修飾特征的重要性評(píng)估

由于在3種分類器中，隨機(jī)森林(RF)和支持向量機(jī)(SVM)2種模型的預(yù)測(cè)能力更為突出，因此本研究通過隨機(jī)森林和支持向量機(jī)的特征重要性，衡量擬南芥中可能對(duì)DNA復(fù)制早晚具有影響的表觀遺傳修飾。圖5展示了隨機(jī)森林和支持向量機(jī)2種分類器特征的重要性，其中，隨機(jī)森林模型的特征重要性通過特征的平均下降精確度大小來衡量。

圖5 基于隨機(jī)森林和支持向量機(jī)分類器的特征重要性Fig.5 Feature importance of RF and SVM classifiers

圖5顯示，H3.1_leaf、H3.3_leaf、H2AW_leaf、H3K36me3_leaf、H3K4me3_leaf、H4K16ac_seedling在2個(gè)模型特征的重要性排序中都位于前10。H3K4me1_leaf、H3K9ac_leaf、H3K27me3_leaf、H3K4me2_leaf也位于RF模型特征重要性的前10，H4K5ac_seedling、polⅡ_seedling、polⅡser2p_seedling、H3K36ac_leaf也位于SVM模型特征重要性的前10，說明這些表觀遺傳修飾特征可能影響DNA的復(fù)制時(shí)間。

3 討 論

本研究通過整合擬南芥多種表觀遺傳修飾數(shù)據(jù)，構(gòu)建3種分類器(隨機(jī)森林、多類別邏輯回歸、支持向量機(jī))，實(shí)現(xiàn)對(duì)擬南芥基因組DNA復(fù)制時(shí)間的預(yù)測(cè)。結(jié)果顯示，3種預(yù)測(cè)模型均能夠準(zhǔn)確地對(duì)DNA復(fù)制早期、中期、晚期進(jìn)行預(yù)測(cè)和分類。

本研究構(gòu)建的3種分類器均具有較高的預(yù)測(cè)和分類性能，表明擬南芥基因組中DNA的復(fù)制時(shí)間可以利用表觀遺傳修飾進(jìn)行準(zhǔn)確地預(yù)測(cè)和分類，這與模式生物果蠅中的研究結(jié)果[10]有相似之處。通過對(duì)各種表觀遺傳修飾與DNA復(fù)制時(shí)間的相關(guān)性分析發(fā)現(xiàn)，一些表觀遺傳修飾特征與復(fù)制時(shí)間相關(guān)，如大部分的乙?；揎?、染色質(zhì)開放程度與復(fù)制早期信號(hào)呈正相關(guān)，與復(fù)制晚期信號(hào)呈負(fù)相關(guān)，這些結(jié)果與已有研究報(bào)道相吻合(如H3K4me1)[10]；相反，異染色質(zhì)的表觀遺傳修飾特征與復(fù)制早期信號(hào)呈負(fù)相關(guān)，與復(fù)制晚期信號(hào)呈正相關(guān)，如H3K9me2(與果蠅上的報(bào)道一致)和H2A.W[10,21]。模型的特征重要性分析揭示，在DNA復(fù)制時(shí)間的表觀遺傳調(diào)控過程中，存在可能扮演重要作用的組蛋白修飾以及組蛋白變體，比如H3.1、H3.3、H3K4me3等。其中H3.1是典型的H3變體，主要在S期表達(dá)，并以依賴DNA復(fù)制的方式整合到染色質(zhì)中[17]；H3.3位于常染色質(zhì)的轉(zhuǎn)錄活性區(qū)域內(nèi)，包括基因內(nèi)和啟動(dòng)子區(qū)域[32]，其是預(yù)測(cè)DNA的復(fù)制時(shí)間的重要特征。此外，復(fù)制早期信號(hào)與RNA聚合酶結(jié)合位點(diǎn)、啟動(dòng)子區(qū)域相關(guān)的特征(H3K4me3)[19]，以及激活型表觀遺傳修飾特征(H3K36me3、H3K4me3、H3K9ac、H4K16ac)呈正相關(guān)[19-20]。這可能是因?yàn)閿M南芥復(fù)制早期的DNA位于基因組內(nèi)高基因密度區(qū)域，其他物種中也存在類似現(xiàn)象[33]。

相對(duì)于動(dòng)物和人類的研究，植物的DNA復(fù)制時(shí)間以及表觀遺傳修飾的相關(guān)研究尚處于起步階段[7]，這不僅表現(xiàn)在相關(guān)的調(diào)控機(jī)制尚不明確，也表現(xiàn)在相關(guān)組學(xué)數(shù)據(jù)的匱乏。本研究雖然盡可能全面地收集并整理了擬南芥中已發(fā)表和公布的多種表觀遺傳修飾數(shù)據(jù)，但是依然無法盡數(shù)收集到針對(duì)每個(gè)特定組織的表觀遺傳修飾的數(shù)據(jù)，因此無法構(gòu)建組織特異性的預(yù)測(cè)模型。本研究所用的擬南芥復(fù)制時(shí)間數(shù)據(jù)是多個(gè)組織的平均值[7]，但是使用的表觀遺傳修飾數(shù)據(jù)卻具有組織特異性，這會(huì)導(dǎo)致模型的性能下降。但本研究通過分析所構(gòu)建的機(jī)器學(xué)習(xí)模型對(duì)應(yīng)的表觀遺傳修飾特征的重要性，獲得了一些在DNA復(fù)制時(shí)間表觀遺傳調(diào)控過程中可能扮演重要作用的組蛋白修飾以及組蛋白變體，但更為具體的DNA復(fù)制時(shí)間表觀遺傳修飾機(jī)制還需要進(jìn)一步研究。

4 結(jié) 論

本研究構(gòu)建了一種高精度的基于DNA表觀遺傳修飾預(yù)測(cè)基因組DNA復(fù)制時(shí)間的機(jī)器學(xué)習(xí)模型，結(jié)果表明，基因組DNA復(fù)制時(shí)間可以通過表觀遺傳修飾進(jìn)行準(zhǔn)確預(yù)測(cè)，且對(duì)DNA復(fù)制晚期的預(yù)測(cè)最為準(zhǔn)確。擬南芥基因組DNA復(fù)制時(shí)間可能受到精細(xì)的表觀遺傳調(diào)控，其中H3.1、H3.3、H2AW、H4K16ac、H3K36me3、H3K4me3等可以作為后續(xù)驗(yàn)證的重要組蛋白變體及表觀遺傳修飾。