本發(fā)明屬于DNA甲基化檢測芯片的擴(kuò)展
技術(shù)領(lǐng)域：
，更為具體地講，涉及一種DNA甲基化芯片數(shù)據(jù)的擴(kuò)展方法。
背景技術(shù)：
：作為人類基因組最為典型的表觀遺傳現(xiàn)象，DNA甲基化在多種關(guān)鍵生理活動中扮演重要角色。其甲基化狀態(tài)與各種疾病，特別是癌癥的發(fā)生密切相關(guān)。DNA甲基化檢測方法最常用的方法是亞硫酸鹽測序、亞硫酸氫微陣列和基于濃縮的方法。亞硫酸鹽測序擁有最全面的基因組覆蓋，但高測序深度使它變的非常昂貴。雖然亞硫酸氫微陣列和基礎(chǔ)的濃縮的方法的成本相對較低，但是基于亞硫酸氫微陣列的平臺被先驗?zāi)繕?biāo)區(qū)域所約束，基于濃縮的方法具有相對低的分辨率和高敏感性實驗偏差。在450K芯片檢測已經(jīng)廣泛用于人體組織，尤其是在幾千患者樣本的DNA甲基化檢測中，因此將450K芯片預(yù)測的覆蓋范圍的擴(kuò)展將是非常有效的最大化使用現(xiàn)有的數(shù)據(jù)的方案。技術(shù)實現(xiàn)要素：本發(fā)明的目的在于克服現(xiàn)有技術(shù)的不足，提供一種DNA甲基化芯片數(shù)據(jù)的擴(kuò)展方法，通過DNA甲基化芯片顯著的擴(kuò)大檢測的覆蓋范圍，有效的最大化使用現(xiàn)有的數(shù)據(jù)。為實現(xiàn)上述發(fā)明目的，本發(fā)明一種DNA甲基化芯片數(shù)據(jù)的擴(kuò)展方法，其特征在于，包括以下步驟：(1)、數(shù)據(jù)獲取從公共數(shù)據(jù)庫中獲取現(xiàn)有的任意一種組織的每一個待預(yù)測CpG位點的上游和下游d1長度范圍內(nèi)基于DNA甲基化芯片測得的甲基化值1及此組織的DNA序列，以及基于全基因組亞硫酸氫鈉測序法測得的每一個待預(yù)測CpG位點甲基化值2；(2)、整合DNA甲基化芯片測得的甲基化值1及此組織的DNA序列作為預(yù)測模型的特征對甲基化值1進(jìn)行加權(quán)，并將此加權(quán)值作為預(yù)測模型的特征1；在DNA序列中，統(tǒng)計出每一個待預(yù)測CpG位點相鄰區(qū)域d2長度范圍內(nèi)的1聚體和2聚體的DNA堿基對出現(xiàn)次數(shù)，作為預(yù)測模型的特征2；統(tǒng)計NpN的產(chǎn)生頻率，作為預(yù)測模型的特征3；最后將特征1、特征2和特征3共同作為預(yù)測模型的輸入特征；(3)、特征變換(3.1)、對輸入特征進(jìn)行標(biāo)準(zhǔn)化處理；設(shè)輸入特征是一個n行p列的矩陣X＝(xi,j)n×p進(jìn)行標(biāo)準(zhǔn)化變換：其中，表示第j列的平均值，sj表示第j列的標(biāo)準(zhǔn)差，對輸入特征進(jìn)行標(biāo)準(zhǔn)化處理后得到標(biāo)準(zhǔn)化矩陣Z；(3.2)、計算標(biāo)準(zhǔn)化矩陣Z的相關(guān)系數(shù)矩陣R：(3.3)、按照表示預(yù)設(shè)閾值，計算相關(guān)系數(shù)矩陣R的特征方程|R-λIp|＝0的特征根λk，Ip單位矩陣，確定出m個主成分分量；對每個特征根λk解方程Rbk＝λkbk，求得單位特征向量bk；(3.4)、將標(biāo)準(zhǔn)化矩陣Z轉(zhuǎn)換為主成分Uk＝ZTbk，k＝1,2,…,m，Uk表示第k個主成分；(4)、構(gòu)建隨機(jī)森林回歸模型先利用m個主成分分量和甲基化值2構(gòu)建多個決策樹，再將多個決策樹按照其產(chǎn)生方式組合成隨機(jī)森林回歸模型；(5)、獨立樣本預(yù)測(5.1)、按照步驟(1)所述方法獲取任意組織的每一個待預(yù)測CpG位點的上游和下游d1長度范圍內(nèi)基于DNA甲基化芯片測得的甲基化值及此組織的DNA序列；(5.2)、按照步驟(2)-(3)的方法整合DNA甲基化芯片測得的甲基化值及此組織的DNA序列得到輸入特征，然后對輸入特征進(jìn)行特征變換得到m個主成分；(5.3)、將m個主成分輸入步驟(4)訓(xùn)練好的隨機(jī)森林模型進(jìn)行預(yù)測得到每個待預(yù)測CpG位點的甲基化值即完成對DNA甲基化芯片數(shù)據(jù)的擴(kuò)展。本發(fā)明的發(fā)明目的是這樣實現(xiàn)的：本發(fā)明一種DNA甲基化芯片數(shù)據(jù)的擴(kuò)展方法，通過預(yù)測DNA甲基化芯片未覆蓋的CpG位點來實現(xiàn)DNA甲基化芯片數(shù)據(jù)的擴(kuò)展。具體的講，對于待預(yù)測的CpG位點，先基于DNA甲基化芯片測得的數(shù)據(jù)和DNA序列提取特征，然后進(jìn)行特征變換并結(jié)合全基因組亞硫酸氫鈉測序法測得的待預(yù)測點的甲基化值訓(xùn)練隨機(jī)森林回歸模型，最后使用訓(xùn)練好的隨機(jī)森林回歸模型對新數(shù)據(jù)進(jìn)行預(yù)測。同時，本發(fā)明一種DNA甲基化芯片數(shù)據(jù)的擴(kuò)展方法還具有以下有益效果：(1)、本發(fā)明充分利用了現(xiàn)有芯片數(shù)據(jù)，對于挖掘出與疾病相關(guān)的重要信息具有重大意義。(2)、我們的模型取得了令人滿意的性能，它優(yōu)于現(xiàn)有的常用方法EpiGRAPH。此外，在新的細(xì)胞類型的甲基化預(yù)測中我們展示了預(yù)測模型的泛化能力。由于公開發(fā)表的DNA甲基化芯片檢測數(shù)據(jù)較多，因此該模型可以實現(xiàn)對多個組織的全基因組甲基化水平的預(yù)測。附圖說明圖1是一種DNA甲基化芯片數(shù)據(jù)的擴(kuò)展方法流程圖；圖2是預(yù)測模型特征選擇簡圖；圖3是預(yù)測模型的相關(guān)系數(shù)及一致性示意圖；圖4是預(yù)測模型在H9上的預(yù)測結(jié)果圖。具體實施方式下面結(jié)合附圖對本發(fā)明的具體實施方式進(jìn)行描述，以便本領(lǐng)域的技術(shù)人員更好地理解本發(fā)明。需要特別提醒注意的是，在以下的描述中，當(dāng)已知功能和設(shè)計的詳細(xì)描述也許會淡化本發(fā)明的主要內(nèi)容時，這些描述在這里將被忽略。實施例為了方便描述，先對具體實施方式中出現(xiàn)的相關(guān)專業(yè)術(shù)語進(jìn)行說明：CpG位點：CG堿基同時出現(xiàn)的位點；450K甲基化芯片：一種DNA甲基化芯片；bp：用來表示DNA長度的單位,也就是我們常說的堿基數(shù)；pearson相關(guān)系數(shù)：皮爾森相關(guān)系數(shù)；Spearman相關(guān)系數(shù)：斯皮爾曼相關(guān)系數(shù)；1聚體：物質(zhì)單一狀態(tài)時具有特定的性質(zhì)或功能,不隨多少改變，這里指A、C、G、T四種堿基；2聚體：相同或同一種類的物質(zhì),以成雙的型態(tài)出現(xiàn),可能具有單一狀態(tài)時沒有的性質(zhì)或功能，這里指A/C/G/T兩兩出現(xiàn)的組合；NpN：與CpG相同的意思，N＝A/C/G/T，這里指A/C/G/T兩兩出現(xiàn)的組合；圖1是一種DNA甲基化芯片數(shù)據(jù)的擴(kuò)展方法流程圖。在本實施例中，采用了450K甲基化芯片數(shù)據(jù)，并結(jié)合CpG位點周圍500bp的序列中選取1聚體和2聚體的堿基特征和NpN比率，對全基因組的CpG位點甲基化水平進(jìn)行了預(yù)測。如圖1所示，本發(fā)明一種450K甲基化芯片數(shù)據(jù)的擴(kuò)展方法，具體包括以下步驟：(1)、數(shù)據(jù)獲取從公共數(shù)據(jù)庫中獲取胚胎干細(xì)胞(embryonicstemcells，ESC)中的H1的每一個待預(yù)測CpG位點的上游和下游d1＝1000bp長度范圍內(nèi)基于450K甲基化芯片測得的甲基化值1及此組織的DNA序列，以及基于全基因組亞硫酸氫鈉測序法測得的每一個待預(yù)測CpG位點甲基化值2；圖2(a)所示，待預(yù)測位點表示為空心點，450K甲基化芯片測得的位點為實心點；其中，每個CpG位點的甲基化水平表示為甲基化率；每個CpG位點的的甲基化值被描述成一個beta值從0(非甲基化)到1(完全甲基化)。(2)、整合450K甲基化芯片測得的甲基化值1及此組織的DNA序列作為預(yù)測模型的特征對甲基化值1進(jìn)行加權(quán)，并將此加權(quán)值作為預(yù)測模型的特征1；在DNA序列中，統(tǒng)計出每一個待預(yù)測CpG位點相鄰區(qū)域d2＝500bp長度范圍內(nèi)的1聚體和2聚體的DNA堿基對出現(xiàn)次數(shù)，作為預(yù)測模型的特征2；統(tǒng)計NpN的產(chǎn)生頻率，作為預(yù)測模型的特征3；最后將特征1、特征2和特征3共同作為預(yù)測模型的輸入特征；在本實施例中，對于每一個的CpG位點，收集相鄰區(qū)域中1000個堿基對的450K芯片檢測數(shù)據(jù)，將450K芯片檢測的此1000個堿基對中的CpG位點的加權(quán)甲基化值作為一個特征；將一個CpG位點相鄰區(qū)域中500bp范圍內(nèi)的1聚體和2聚體的DNA堿基對共20個特征，以及NpN(N＝A/C/G/T)共16個特征，如圖2(b)所示，最后共同構(gòu)成37個特征作為模型構(gòu)建的輸入特征。(3)、特征變換(3.1)、對輸入特征進(jìn)行標(biāo)準(zhǔn)化處理；設(shè)輸入特征是一個n行p列的矩陣X＝(xi,j)n×p進(jìn)行標(biāo)準(zhǔn)化變換：其中，表示第j列的平均值，sj表示第j列的標(biāo)準(zhǔn)差，對輸入特征進(jìn)行標(biāo)準(zhǔn)化處理后得到標(biāo)準(zhǔn)化矩陣Z；(3.2)、計算標(biāo)準(zhǔn)化矩陣Z的相關(guān)系數(shù)矩陣R：(3.3)、按照計算相關(guān)系數(shù)矩陣R的特征方程|R-λIp|＝0的特征根λk，Ip單位矩陣，確定出15個主成分分量；對每個特征根λk解方程Rbk＝λkbk，求得單位特征向量bk；(3.4)、將標(biāo)準(zhǔn)化矩陣Z轉(zhuǎn)換為主成分Uk＝ZTbk，k＝1,2,…,15，Uk表示第k個主成分；(4)、構(gòu)建隨機(jī)森林回歸模型先利用m個主成分分量和甲基化值2構(gòu)建多個決策樹，再將多個決策樹按照其產(chǎn)生方式組合成隨機(jī)森林回歸模型；下面進(jìn)行詳細(xì)說明：(4.1)、構(gòu)建決策樹設(shè)訓(xùn)練樣本向量的維度是F維，即訓(xùn)練集有F個屬性。在構(gòu)建開始之前選定一個參數(shù)f，滿足f＜＜F，在構(gòu)建每個內(nèi)部節(jié)點的過程中，都需要從訓(xùn)練集中采用隨機(jī)抽樣的方法從他的所有F個屬性選取f個屬性，然后從f個屬性中根據(jù)信息增益比，選出一個最優(yōu)的屬性充當(dāng)分裂屬性，進(jìn)而是決策在此節(jié)點產(chǎn)生分裂。信息增益比的計算采用如下公式：其中S為全部樣本集合，value(T)是屬性T所有取值的集合，v是T的其中一個屬性值，Sv是S中屬性T的值為V的樣例集合，|Sv|為Sv中所含樣例數(shù)。Entropy(Sv)即表示信息增益，他的計算采用如下公式：其中，就是類別的總數(shù)，類別C是變量，它的取值是C1,C2,...,Cn，而每一個類別出現(xiàn)的概率分別是P(C1),P(C2),...,P(Cn)。(4.2)構(gòu)建隨機(jī)森林回歸模型隨機(jī)森林其實是由很多決策樹組合而成的，但是其回歸模型的性能往往依賴于組合是按照一種什么樣的準(zhǔn)則進(jìn)行的，每一棵決策樹的樣本都取自一個訓(xùn)練集，都依賴于一個其產(chǎn)生方式的規(guī)則，這個規(guī)則往往用一個隨機(jī)向量來表示。它的產(chǎn)生方式的遞歸描述如下：首先為第q棵決策樹生成一個決定了其生成過程的隨機(jī)向量θq，θq需要滿足與之前生成的q-1個隨機(jī)向量獨立同分布。然后在原始訓(xùn)練樣本X中利用隨機(jī)抽樣方法進(jìn)行抽樣，這棵決策樹的生成過程中的數(shù)據(jù)都取自這次抽樣的結(jié)果Xq。采用上述策略，隨機(jī)向量θq和抽樣結(jié)果Xq就能夠生成第q棵決策樹h(Xq,θq)。在對樣本集經(jīng)過q輪抽樣之后，一共可以得到q棵決策樹。當(dāng)輸入一個新樣本時，隨機(jī)森林中的每一顆決策樹分別進(jìn)行判斷，最后對每顆決策樹的結(jié)果取平均值。(4.3)、模型性能的評估；交叉檢驗法是十分常用的模型驗證方法。其原理是將訓(xùn)練樣本分成容量相同的w個子集，并對模型訓(xùn)練w次。在第u次(u＝1,2,L,w)訓(xùn)練時,要用除了第u個子集的所有子集訓(xùn)練模型，再用得到的模型對第u個子集計算誤差。以w次誤差的平均數(shù)值作為模型推廣能力的近似數(shù)值。對于預(yù)測模型性能指標(biāo)我們采用相關(guān)系數(shù)(Spearman相關(guān)系數(shù)和Spearman相關(guān)系數(shù))、一致性、特異性(SP)、靈敏度(SE)、準(zhǔn)確性(ACC)和馬修相關(guān)系數(shù)(MCC)來進(jìn)行評估。兩個變量之間的Pearson相關(guān)系數(shù)定義為這兩個變量的協(xié)方差與二者標(biāo)準(zhǔn)差積的商，即：其中，和μY分別是和Y的平均值，和σY分別是和Y的標(biāo)準(zhǔn)差，和Y分別代表擬合的甲基化值和實際WGBS記錄的甲基化值。Spearman相關(guān)系數(shù)的計算公式為：其中，n為樣本集大小，n行甲基化值和經(jīng)過等級排序后的值為一致性是預(yù)測值與實際值之間的差值小于0.25的百分比。SE，SP，ACC和MCC的計算公式如下：其中TP表示預(yù)測正確的正樣本(true-positive)；TN表示預(yù)測正確的負(fù)樣本(true-negative)；FP表示預(yù)測錯誤的負(fù)樣本(false-positive)；FN表示預(yù)測錯誤的正樣本(false-negative)。通過使用3倍交叉驗證測試20次，獲得預(yù)測模型在每條染色體上的平均性能。預(yù)測值和真實值的Pearson相關(guān)系數(shù)和Spearman相關(guān)系數(shù)示于圖3(a)，一致性示于圖3(b)。在圖3(a)和3(b)中結(jié)果最好的是在18號染色體上，相關(guān)系數(shù)為0.91(0.82)，一致性為88％；結(jié)果最差的是在Y染色體上，相關(guān)系數(shù)為0.84(0.73)，一致性為82％。(4.4)、方法對比；目前預(yù)測甲基化水平最常用的方法是EpiGRAPH，用我們的模型與其進(jìn)行對比。把預(yù)測值經(jīng)過二值化處理，當(dāng)CpG位點的預(yù)測的甲基化值大于0.5的時候，我們認(rèn)為的其甲基化狀態(tài)為+1，反之當(dāng)預(yù)測值小于0.5時，我們認(rèn)為其甲基化之為-1。在結(jié)果最好與最壞的染色體上，用我們的模型與現(xiàn)有常用方法對比。如表1所示，我們構(gòu)建的模型在18號染色體上，結(jié)果優(yōu)于現(xiàn)有常用方法EpiGRAPH。如表2所示，我們構(gòu)建的模型在Y染色體上，結(jié)果也優(yōu)于現(xiàn)有常用方法EpiGRAPH。表1是預(yù)測模型在18號染色體上與現(xiàn)有方法EpiGRAPH的對比結(jié)果；Chr18SESPACCMCCRFinEpiGRAPH0.900.820.860.73RFinourmodel0.940.920.930.86表1表2是預(yù)測模型在Y染色體上與現(xiàn)有方法EpiGRAPH的對比結(jié)果；ChrYSESPACCMCCRFinEpiGRAPH0.960.260.820.33RFinourmodel0.920.730.850.74表2(5)、獨立樣本預(yù)測(5.1)、按照步驟(1)所述方法獲取H9組織的每一個待預(yù)測CpG位點的上游和下游d1＝1000bp長度范圍內(nèi)基于DNA甲基化芯片測得的甲基化值及此組織的DNA序列；(5.2)、按照步驟(2)-(3)的方法整合450K甲基化芯片測得的甲基化值及此組織的DNA序列得到輸入特征，然后對輸入特征進(jìn)行特征變換得到15個主成分；(5.3)將15個主成分輸入步驟(4)訓(xùn)練好的隨機(jī)森林模型進(jìn)行預(yù)測得到每個待預(yù)測CpG位點的甲基化值即完成對DNA甲基化芯片數(shù)據(jù)的擴(kuò)展。為了驗證結(jié)果，下載對應(yīng)的基于全基因組亞硫酸氫鈉測序法測得的每一個待預(yù)測CpG位點甲基化值，在H9組織上預(yù)測的性能指標(biāo)如圖4所示，除了X染色體以外，預(yù)測值和真實值的Pearson相關(guān)系數(shù)為0.80±0.01，如圖4(a)所示一致性為88％±1％，如圖4(b)所示。預(yù)測的甲基化值在經(jīng)過二進(jìn)制處理之后，我們計算該預(yù)測結(jié)果的SE，SP，ACC和MCC,結(jié)果如圖4(c)所示，得到的結(jié)果是令人滿意的。盡管上面對本發(fā)明說明性的具體實施方式進(jìn)行了描述，以便于本
技術(shù)領(lǐng)域：
的技術(shù)人員理解本發(fā)明，但應(yīng)該清楚，本發(fā)明不限于具體實施方式的范圍，對本
技術(shù)領(lǐng)域：
的普通技術(shù)人員來講，只要各種變化在所附的權(quán)利要求限定和確定的本發(fā)明的精神和范圍內(nèi)，這些變化是顯而易見的，一切利用本發(fā)明構(gòu)思的發(fā)明創(chuàng)造均在保護(hù)之列。當(dāng)前第1頁1 2 3