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