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