表觀遺傳是指在基因遺傳序列不變的前提下，外界和內在環境因素對基因產生修飾作用，導致基因表達水平或功能發生改變，進而對多種表型或疾病結局產生影響，近年來愈發受到關注。其中，DNA 的甲基化修飾已被證實與人類發育和疾病發生發展密切相關。而全基因組尺度的甲基化檢測所產生的高維度基因組學數據，能夠較為全面地體現基因組層面整體和局部的表觀遺傳修飾情況，已成為該領域的主要研究內容之一。本文基于全基因組甲基化芯片數據，總結該類組學數據的質控流程，常見表觀遺傳組學關聯性的統計分析方法和思路，以及基于 SAS JMP Genomics 10 軟件的主要結果的可視化實現，為同類研究提供參考。
引用本文： 樊娟娟, 吉昕宇, 沈思鵬, 張汝陽, 魏永越, 陳峰. 全表觀基因組關聯研究的數據質控、分析流程和可視化. 中國循證醫學雜志, 2021, 21(6): 721-728. doi: 10.7507/1672-2531.202012149 復制
在繁雜的生命系統中，生命進程中的各種事件不僅是由其基因型決定，同時還受到外界環境因素的共同調控。表觀遺傳學（epigenetics）是指在生物基因的 DNA 序列不變的前提下，在外界和內在環境因素的驅動下，基因表達水平與功能發生改變，近年來受到越來越多的關注。表觀遺傳學可以通過對細胞分裂的基因調控作用而對生物體的表現型產生影響，因此需要通過與全基因組關聯研究類似的大規模、系統性的研究—全表觀基因組關聯研究（epigenome-wide association studies，EWAS）來分析。EWAS 研究和全基因組關聯分析（genome-wide association study，GWAS）研究類似，是探索表觀遺傳變異與常見疾病/表型之間關聯的有效工具[1]，能夠為我們提供表觀基因組擾動與人類疾病相關性狀之間關聯的信息。表觀遺傳修飾中，DNA 甲基化（DNA methylation）與人類發育和腫瘤疾病密切相關，已成為表觀遺傳學的主要研究內容[2]。哺乳動物的 DNA 甲基化通常發生在 DNA 鏈中胞嘧啶-鳥嘌呤雙核苷酸（cytosine-phosphate-guanine pairs of nucleotides，CpGs）的胞嘧啶上，在 DNA 甲基化轉移酶的催化作用下生成 5-甲基胞嘧啶（5-mC）[3]。DNA 甲基化改變的是基因表達水平而不是 DNA 序列，但其在人類疾病的發生發展中起重要作用[4]。它能夠影響許多生物學進程，如發育和細胞分化，也在包括癌癥和哮喘等復雜疾病的進展中發揮關鍵作用[5]。因此，探索這種表觀遺傳變異與疾病/表型之間的關聯在全基因組研究中越來越受到青睞。本文將從甲基化數據質控和數據分析思路兩部分來具體論述。
1 甲基化數據的質控流程
DNA 甲基化檢測可通過基于全基因組測序[6]和全基因組微陣列芯片[7]兩種技術來實現。本文以 Illumina Human Methylation 450 BeadChip 芯片（簡稱甲基化 450k 芯片）為例，闡述 EWAS 芯片數據質控流程。
讀取芯片檢測的原始圖片數據，計算各位點的甲基化水平。由于甲基化水平近似 Beta 分布，因此常用 β 值來表示該位點甲基化的探針強度占總探針強度的比例[8]，其分布特點是雙峰、非對稱且方差不齊，取值范圍為 0～1，β 值越大，甲基化比例越高。計算公式為：
 |
其中，M 表示匹配到該位點的所對應的所有探針中，檢測到甲基化信號的強度，U 表示未檢測到甲基化的探針信號強度。數據質控主要針對 β 值進行，由六個部分組成：背景校正、位點質控、樣本質控、分位數標化、探針類型校正和批次效應校正，具體流程如圖 1 所示。
圖1
甲基化 450k 芯片數據質控流程
1.1 背景校正（background correction）
芯片的探針信號中普遍存在背景噪聲，需要進行背景降噪（subtract background），以提高甲基化信號的可信度。在芯片設計之初，設計者添加了針對無甲基化區域的陰性探針組，以此來估計背景噪聲。實際檢測位點的甲基化信號強度即為原始的信號強度減去該背景噪聲強度。
1.2 位點質控（quality control probes）
檢測 P 值（detection P-value）代表著一個探針在某個樣本檢測質量的好壞，P 值越小代表著探針檢測質量越好。negative 探針是一系列質量差的探針信號，用于計算 P 值。若檢測 P 值大于一定閾值（一般取 0.05 或 0.01，這里以 0.05 為例），則將此探針視為無法區分靶信號和背景信號的低質量探針。位點質控通常包括以下幾步：① 剔除在超過 5% 的樣本中檢測 P 值>0.05 的低質量 CpG 位點；② 剔除變異系數<5% 的 CpG 點，變異系數小意味著無信號波動，即 CpG 位點的甲基化程度穩定；③ 剔除 CpG 附近含有單核苷酸多態性（single nucleotide polymorphism，SNP）的探針，因其檢測出的甲基化水平可能會受到 SNP 的較大影響；④ 剔除性染色體（sex-chromosome）探針。一般僅保留 1～22 號染色體上的 CpG 位點。當分組比較的樣本混合了不同性別，則無法正確衡量性染色體上的甲基化水平的差異；⑤ 剔除交叉反應性（cross-reactive）的 CpG 位點。探針的交叉反應性：當實際靶標與預期靶標高度同源時，探針與實際靶標雜交，因此可能產生假信號，從而導致無效結論[9]。
1.3 樣本質控（quality control samples）
剔除含有超過 5% 未檢出信號或低質量 CpG 位點的樣本。
1.4 分位數標化（quantile normalization）
為消除技術變異對數據所帶來的影響，使所有樣本中探針強度的分布具有一致的分位數，使數據具有較好的可比性。將單一探針的強度用矩陣中所有有相同秩次的探針平均值來替代，使得每個樣本的分布相同。
1.5 探針類型校正（bias correction for type I/II probes）
450k 芯片包括兩種類型探針[10]，兩種探針設計方式的不同，導致兩種探針的信號分布不一致，如果將兩類探針 CpG 位點檢測結果作為整體去分析，所得結果將不準確。Teschendorff 等提出了一種 β 混合分位數校正法[11]，將 2 型探針 β 值的分布變換為近似 1 型探針的分布，以盡可能消除兩種探針之間的差異。
1.6 批次效應校正（batch effects adjustment）
批次效應是不同批次樣本組之間的系統誤差，是由與生物學無關因素的變異而導致的數據誤差。主要由樣品處理和制備方面的差異（如提取方法），DNA 處理（如擴增、標記、雜交），陣列/芯片掃描（如背景噪聲），樣品在芯片上的位置，技術偏差等原因導致，批次效應校正一般使用基于經驗貝葉斯的方法—ComBat[12]。
此外，有研究發現 450k 芯片中許多 CpGs 都存在較大的測量誤差，降低了對關聯關系的統計能力，從而導致統計效能降低[13]。Jun 等提出了一種 450k 芯片所特有的校正測量誤差的方法—CpGFilter[14]，以降低檢測技術所引入的偏倚。基于線性混合效應模型來估計組內相關系數（intra-class correlation coefficient，ICC），進而過濾測量誤差較大的 CpG。
經過以上 6 步質控，得到可供統計分析的甲基化數據集。
2 甲基化數據的分析思路
對近幾年國內外學者提出的常見甲基化數據統計分析方法進行了介紹和探討，將從研究目的角度，系統地介紹每一類研究的背景和基本統計思想。
2.1 識別差異甲基化位點和區域的統計分析方法
差異甲基化區域（differentially methylated region，DMR）指的是不同組的樣本（例如正常細胞和腫瘤細胞、不同疾病亞型等）甲基化水平具有明顯的變化，為要識別這些位點或區域[8]。常見的統計學分析方法包含：
2.1.1 基本方法
參數檢驗：假設甲基化數據服從特定參數分布，檢驗正常人和腫瘤患者（或癌組織 vs.癌旁組織）甲基化水平差別；如：t 檢驗，方差分析。非參數檢驗：不依賴于數據的分布形式，檢驗甲基化分布的差異；如：Wilcoxon 秩和檢驗。
2.1.2 基于均勻—正態混合分布模型的似然比檢驗
甲基化 β 值呈雙峰分布，在分子水平，位點的甲基化狀態有三種：未甲基化、完全甲基化和半甲基化[15]。因此，Wang[16]提出了基于均勻-正態混合分布模型的似然比檢驗方法，來識別病例和對照間的差異甲基化位點。根據甲基化位點的三種不同狀態，通過三組混合分布（兩均勻分布和一截斷的正態分布）來模擬甲基化數據。通過混合分布的概率和正態分布的均數來檢驗差異甲基化位點。
2.1.3 回歸模型
大多數研究假設甲基化數據服從正態分布，以甲基化信號作為因變量，以組別或其他潛在研究因素作為自變量，進行線性回歸，以評估研究因素和甲基化信號之間的關聯性；亦可以將組別作為因變量，將甲基化信號作為自變量，進行 logistic 回歸。亦有研究假設甲基化數據不依賴數據分布形式，可采用非參數回歸模型（如：分位數回歸）。
2.1.4 聚類
Yip 等[17]提出空間聚類（spatial clustering method，SCM）的方法來尋找基因組中與疾病有關的候選 DMR。利用位點間的距離信息，在關聯研究中提高檢驗效能，篩選 DMR。
2.1.5 貝葉斯方法
Teng 等[18]構建了 5 種經驗貝葉斯模型比較不同方法分析差異甲基化效果的好壞。絕大部分經驗貝葉斯模型基于 gamma 分布或對數正態分布的先驗假設，借助 E-M 算法進行參數估計，進而尋找出差異甲基化位點。面對超高維甲基化數據，貝葉斯方法具有符合生物學假設、精度高、檢驗效能高等特點。針對不同分組 CpG 位點甲基化水平方差不齊和樣本量較小的問題，Feng 等[19]提出了貝葉斯分層模型的方法來解決。也可采用 Beta-二項分布分層模型來解釋個體內和個體間兩方面的變異。
差異甲基化分析的結果可用火山圖（volcano plot）和瀑布圖（bar plot）等展示，見圖 2 和圖 3。
圖2
火山圖展示組間差異甲基化分析結果
圖中每個點代表一個甲基化位點；橫軸表示差異甲基化程度，即：fold change；縱軸為差異甲基化統計學假設檢驗的
圖3
瀑布圖
該圖表示甲基化位點在兩組間的差異性情況。橫軸為甲基化位點的差異大小，縱軸為各個甲基化位點。SAS JMP Genomics 10 軟件實現：圖形→圖形生成器。
2.2 結局關聯甲基化位點的識別策略
類似于 GWAS 的思想，甲基化關聯性分析旨在全基因組或部分區域尋找與結局相關的甲基化位點，全面揭示疾病發生、發展與治療相關的表觀遺傳標志物。常見統計方法有：① 單位點的關聯性分析：最基本的 EWAS 關聯性分析是評估每個 CpG 位點與研究結局的關聯性強弱。其中，Cox 比例風險模型是生存分析的常用方法，logistic 回歸模型則是關聯疾病發生風險或鑒別診斷的常用模型。其檢驗結果通常使用 Bonferroni 等方法進行多重比較校正。② CpG 集合的關聯性分析：核函數關聯檢驗是有監督的、靈活的、計算效率較高的回歸模型，采用方差成分檢驗，來提高檢驗效能。自由度的調整取決于 CpG 位點間的相關性結構，若 CpG 位點間呈現高度相關，檢驗的自由度就會大大減小[20]。主要用于識別某一區域內的遺傳變異（包括常見變異和罕見變異）與連續性或分類結局變量間的關聯，同時又便于協變量的調整[21]。GSVA（Gene Set Variation Analysis）是一種非參數的無監督分析方法，主要用于評估基于基因集的富集結果，多用于評估樣本在不同生物學通路上的富集程度[22]。ssGSEA（single sample Gene Set Enrichment Analysis）是 GSEA 方法的擴展，可以計算出每個樣本的富集水平[23]。③ 位點間交互作用分析：CpG 位點間的相關結構較為復雜且不固定，變量間還可能存在一階或高階交互作用[20]，交互作用的檢測方法有叉生分析、分層分析、回歸模型。回歸模型是分析交互作用的常用模型，回歸中交互項回歸系數的估計是以存在相乘交互作用為基礎，在分析交互作用時引入一個相乘項進入回歸模型，通過相乘項回歸系數的估計來判斷交互作用是否有統計學意義及其作用大小。關聯性分析結果可用曼哈頓圖（Manhattan plot）展示，見圖 4。
圖4
曼哈頓圖
橫軸表示甲基化位點在染色體上的位置，縱軸表示每個位點關聯性分析假設檢驗
2.3 基于甲基化數據的結局預測
對于有預測價值的 CpG 位點，為預測疾病的發生風險或預后，對臨床輔助治療、靶向藥物研發具有一定的指導意義。此時，需建立預測模型進行分析，在構建預測模型時，通常遵循著數據降維—模型構建的流程。
數據降維：① 在變量降維篩選中，懲罰方法是使用較為廣泛的一類方法，通過對未知參數的值進行壓縮，同時實現變量選擇和參數估計，具有降低估計偏差、提高預測精度和模型可解釋性的優點。在建模過程中，為防止過度擬合，也可添加懲罰項，引入正則化參數。一個基因或一個區域內含多個 CpG 位點，這些位點通常存在較強的、復雜的相關性結構，其中只有一部分是致病位點。根據這一特點，Sun 等[24]提出懲罰 logistic 回歸模型，在考慮 CpG 復雜相關性結構的前提下，篩選出和疾病相關的獨立 CpG 位點。現有的主流懲罰模型還有 Lasso[25]、Elastic net[26]等。② 主成分分析法（principle component analysis，PCA）目前已被廣泛用于變量間存在相關性的多元分析中。PCA 通常作為一種降維的手段。鑒于 PCA 是一種無監督的降維方法，沒有考慮結局變量與自變量的關系，所以無法保證主成分與結局變量是否存在真正的相關。基于此，有學者提出了有監督的主成分分析法（supervised principle component analysis，SPCA）[27]，先根據自變量與結局變量相關性進行排序，篩選相關性最高的幾個自變量提取主成分，然后評估潛在變量與結局變量間的統計學意義。核函數主成分分析法（kernel principle component analysis，KPCA）是一種非線性 PCA，能處理自變量與因變量之間的非線性關系，近年來在機器學習領域被廣泛運用[28]。③ Leo Breiman 于 2001 年提出隨機森林（random forest，RF）的算法[29]，用于處理位點數多而樣本量少的基因組數據，對高維數據進行降維，利用 RF 算法獲得數據集中每個變量的重要性評分排序，達到初篩變量的目的。其優點在于通過重要性評分排序進行了變量篩選，可大大減少噪音變量，消除冗余，因而可用來分析變量數遠大于樣本數的資料。張秋伊等[30]將 RF 算法運用于類風濕性關節炎病例對照研究的高維甲基化數據的分析，結果表明 RF 能綜合位點信息，有效剔除噪音變量，提高檢驗效能。該方法適用于當變量數遠大于樣本量的情況，可以提高檢驗效能。
此外，基于樹的方法，RF 的強化版算法—Ranger，能夠固定模型中的協變量，在實際分析中具有分析優勢。Wei 等[31]通過 Ranger 算法獲得候選 CpG 位點，再繼續進行傳統 Cox 回歸，得到與早期非小細胞肺癌存活相關的 CpG 位點。
模型構建：① 回歸類參數模型：可以對于篩選后的位點建立 logistic 或 Cox 回歸模型，用 OR 或 HR 估計其對疾病風險大小，多位點的綜合效應可以用風險評分（score）進入模型，采用 
估計其對疾病風險的總效應[32]。如果納入變量較多時直接構建參數預測模型，則十分容易導致過度擬合現象（overfitting problems）。故懲罰回歸類方法應運而生，目前已被較多地應用于變量篩選和模型構建，如 LASSO 和彈性網絡（elastic net）等。其通過壓縮回歸系數，將關聯性弱的變量系數壓縮為 0，從而保留了關聯性強的變量[33]。② 基于機器學習的預測分類模型：近些年，機器學習類方法由于模型預測表現較好，受到了極為廣泛的關注。如支持向量機（support vector machine，SVM），已被廣泛用于分類和回歸問題，是一種高效的學習模型。它基于簡單而強大的思想，是一種非常流行的分類方法[34]。RF 是可以用來處理分類、回歸等問題[35]。它有很多優點：具有極高的準確率；隨機性的引入，不容易過度擬合，具有較好的抗噪能力；能處理很高維度的數據，并且無需做特征選擇。K 近鄰作為一種有監督分類算法，是一種基于距離的分類器，在分類方面已經得到了廣泛的研究和討論。其算法主要思想是根據距離相近的鄰居類別，來判斷自己所屬類別。K 近鄰思想簡單，易于理解，易于實現，無需估計參數，無需訓練。樸素貝葉斯（Na?ve Bayes）也是經典的機器學習算法之一，是為數不多的基于概率論的分類算法[36]。其是一個基于貝葉斯定理的比較簡單的概率分類器，原理比較簡單，學習和預測的效率都很高。
建立預后預測模型后，常用生存曲線（Kaplan-Meier curve）、森林圖（forest plot）來對比不同組別生存差異，見圖 5 和圖 6。
圖5
生存曲線
該圖表示了一個甲基化位點隨著隨訪時間的增加發生終點事件概率變化情況。為了比較不同甲基化水平的生存時間差異，通常把該位點進行分組。例如三分類之后，每組的生存曲線如圖所示。SAS JMP Genomics 10 軟件實現：分析→可靠性和生存→生存。
圖6
森林圖
該圖展示甲基化位點在不同亞組中分層分析的結果，藍點和線表示效應值（如風險比）及 95% 可信區間。SAS JMP Genomics 10 軟件實現：圖形→圖形生成器，或插件“Forest Plot”。
3 討論
本文從 EWAS 的數據質控、分析及可視化展開，簡要介紹了甲基化數據規范質控流程和數據分析常用方法，并展示了可視化結果的常見圖形，各部分互相聯系，不可或缺。鑒于目前還未有系統介紹 DNA 甲基化數據處理流程分析策略的文章，本文可為類似研究提供一個參考。對于 EWAS 的關聯性分析和預測模型構建，通常需要進行獨立的外部驗證，確保結果的可重復性、穩健性。在當今的研究中，隨著研究隊列的逐漸增多，對多階段、大樣本數據的外部驗證已成為趨勢[37, 38]。
隨著生物技術的進步，DNA 甲基化的測量手段也在不斷升級。如最近推出的 Infinium Methylation EPIC BeadChip（簡稱 850k 芯片）、Whole Genome Bisulfite Sequencing（甲基化測序）等手段能夠大大提高甲基化位點的識別數量，費用也在逐漸降低，在將來的研究中或許將成為主流。此外，甲基化與基因多態性間具有緊密的聯系，稱為 meQTL；甲基化對基因表達也有著調控，例如可讓基因沉默。如何將基因多態性、DNA 甲基化、基因表達等多種組學的信息整合起來，這將是 GWAS 和 EWAS 統計分析需要進一步探討的問題[39]。
本文所介紹的這些統計分析方法都有各自的適用條件，在實際應用研究中需要結合研究目的、疾病特征、樣本量等具體因素選擇最佳分析方法。
致謝：感謝張瑛、管美建、包文駿博士對 SAS JMP Genomics 軟件可視化內容所提供的建議和幫助。
聲明：所有作者均聲明本研究不存在任何利益沖突。
在繁雜的生命系統中，生命進程中的各種事件不僅是由其基因型決定，同時還受到外界環境因素的共同調控。表觀遺傳學（epigenetics）是指在生物基因的 DNA 序列不變的前提下，在外界和內在環境因素的驅動下，基因表達水平與功能發生改變，近年來受到越來越多的關注。表觀遺傳學可以通過對細胞分裂的基因調控作用而對生物體的表現型產生影響，因此需要通過與全基因組關聯研究類似的大規模、系統性的研究—全表觀基因組關聯研究（epigenome-wide association studies，EWAS）來分析。EWAS 研究和全基因組關聯分析（genome-wide association study，GWAS）研究類似，是探索表觀遺傳變異與常見疾病/表型之間關聯的有效工具[1]，能夠為我們提供表觀基因組擾動與人類疾病相關性狀之間關聯的信息。表觀遺傳修飾中，DNA 甲基化（DNA methylation）與人類發育和腫瘤疾病密切相關，已成為表觀遺傳學的主要研究內容[2]。哺乳動物的 DNA 甲基化通常發生在 DNA 鏈中胞嘧啶-鳥嘌呤雙核苷酸（cytosine-phosphate-guanine pairs of nucleotides，CpGs）的胞嘧啶上，在 DNA 甲基化轉移酶的催化作用下生成 5-甲基胞嘧啶（5-mC）[3]。DNA 甲基化改變的是基因表達水平而不是 DNA 序列，但其在人類疾病的發生發展中起重要作用[4]。它能夠影響許多生物學進程，如發育和細胞分化，也在包括癌癥和哮喘等復雜疾病的進展中發揮關鍵作用[5]。因此，探索這種表觀遺傳變異與疾病/表型之間的關聯在全基因組研究中越來越受到青睞。本文將從甲基化數據質控和數據分析思路兩部分來具體論述。
1 甲基化數據的質控流程
DNA 甲基化檢測可通過基于全基因組測序[6]和全基因組微陣列芯片[7]兩種技術來實現。本文以 Illumina Human Methylation 450 BeadChip 芯片（簡稱甲基化 450k 芯片）為例，闡述 EWAS 芯片數據質控流程。
讀取芯片檢測的原始圖片數據，計算各位點的甲基化水平。由于甲基化水平近似 Beta 分布，因此常用 β 值來表示該位點甲基化的探針強度占總探針強度的比例[8]，其分布特點是雙峰、非對稱且方差不齊，取值范圍為 0～1，β 值越大，甲基化比例越高。計算公式為：
|
其中，M 表示匹配到該位點的所對應的所有探針中，檢測到甲基化信號的強度，U 表示未檢測到甲基化的探針信號強度。數據質控主要針對 β 值進行，由六個部分組成：背景校正、位點質控、樣本質控、分位數標化、探針類型校正和批次效應校正，具體流程如圖 1 所示。
圖1
甲基化 450k 芯片數據質控流程
1.1 背景校正（background correction）
芯片的探針信號中普遍存在背景噪聲，需要進行背景降噪（subtract background），以提高甲基化信號的可信度。在芯片設計之初，設計者添加了針對無甲基化區域的陰性探針組，以此來估計背景噪聲。實際檢測位點的甲基化信號強度即為原始的信號強度減去該背景噪聲強度。
1.2 位點質控（quality control probes）
檢測 P 值（detection P-value）代表著一個探針在某個樣本檢測質量的好壞，P 值越小代表著探針檢測質量越好。negative 探針是一系列質量差的探針信號，用于計算 P 值。若檢測 P 值大于一定閾值（一般取 0.05 或 0.01，這里以 0.05 為例），則將此探針視為無法區分靶信號和背景信號的低質量探針。位點質控通常包括以下幾步：① 剔除在超過 5% 的樣本中檢測 P 值>0.05 的低質量 CpG 位點；② 剔除變異系數<5% 的 CpG 點，變異系數小意味著無信號波動，即 CpG 位點的甲基化程度穩定；③ 剔除 CpG 附近含有單核苷酸多態性（single nucleotide polymorphism，SNP）的探針，因其檢測出的甲基化水平可能會受到 SNP 的較大影響；④ 剔除性染色體（sex-chromosome）探針。一般僅保留 1～22 號染色體上的 CpG 位點。當分組比較的樣本混合了不同性別，則無法正確衡量性染色體上的甲基化水平的差異；⑤ 剔除交叉反應性（cross-reactive）的 CpG 位點。探針的交叉反應性：當實際靶標與預期靶標高度同源時，探針與實際靶標雜交，因此可能產生假信號，從而導致無效結論[9]。
1.3 樣本質控（quality control samples）
剔除含有超過 5% 未檢出信號或低質量 CpG 位點的樣本。
1.4 分位數標化（quantile normalization）
為消除技術變異對數據所帶來的影響，使所有樣本中探針強度的分布具有一致的分位數，使數據具有較好的可比性。將單一探針的強度用矩陣中所有有相同秩次的探針平均值來替代，使得每個樣本的分布相同。
1.5 探針類型校正（bias correction for type I/II probes）
450k 芯片包括兩種類型探針[10]，兩種探針設計方式的不同，導致兩種探針的信號分布不一致，如果將兩類探針 CpG 位點檢測結果作為整體去分析，所得結果將不準確。Teschendorff 等提出了一種 β 混合分位數校正法[11]，將 2 型探針 β 值的分布變換為近似 1 型探針的分布，以盡可能消除兩種探針之間的差異。
1.6 批次效應校正（batch effects adjustment）
批次效應是不同批次樣本組之間的系統誤差，是由與生物學無關因素的變異而導致的數據誤差。主要由樣品處理和制備方面的差異（如提取方法），DNA 處理（如擴增、標記、雜交），陣列/芯片掃描（如背景噪聲），樣品在芯片上的位置，技術偏差等原因導致，批次效應校正一般使用基于經驗貝葉斯的方法—ComBat[12]。
此外，有研究發現 450k 芯片中許多 CpGs 都存在較大的測量誤差，降低了對關聯關系的統計能力，從而導致統計效能降低[13]。Jun 等提出了一種 450k 芯片所特有的校正測量誤差的方法—CpGFilter[14]，以降低檢測技術所引入的偏倚。基于線性混合效應模型來估計組內相關系數（intra-class correlation coefficient，ICC），進而過濾測量誤差較大的 CpG。
經過以上 6 步質控，得到可供統計分析的甲基化數據集。
2 甲基化數據的分析思路
對近幾年國內外學者提出的常見甲基化數據統計分析方法進行了介紹和探討，將從研究目的角度，系統地介紹每一類研究的背景和基本統計思想。
2.1 識別差異甲基化位點和區域的統計分析方法
差異甲基化區域（differentially methylated region，DMR）指的是不同組的樣本（例如正常細胞和腫瘤細胞、不同疾病亞型等）甲基化水平具有明顯的變化，為要識別這些位點或區域[8]。常見的統計學分析方法包含：
2.1.1 基本方法
參數檢驗：假設甲基化數據服從特定參數分布，檢驗正常人和腫瘤患者（或癌組織 vs.癌旁組織）甲基化水平差別；如：t 檢驗，方差分析。非參數檢驗：不依賴于數據的分布形式，檢驗甲基化分布的差異；如：Wilcoxon 秩和檢驗。
2.1.2 基于均勻—正態混合分布模型的似然比檢驗
甲基化 β 值呈雙峰分布，在分子水平，位點的甲基化狀態有三種：未甲基化、完全甲基化和半甲基化[15]。因此，Wang[16]提出了基于均勻-正態混合分布模型的似然比檢驗方法，來識別病例和對照間的差異甲基化位點。根據甲基化位點的三種不同狀態，通過三組混合分布（兩均勻分布和一截斷的正態分布）來模擬甲基化數據。通過混合分布的概率和正態分布的均數來檢驗差異甲基化位點。
2.1.3 回歸模型
大多數研究假設甲基化數據服從正態分布，以甲基化信號作為因變量，以組別或其他潛在研究因素作為自變量，進行線性回歸，以評估研究因素和甲基化信號之間的關聯性；亦可以將組別作為因變量，將甲基化信號作為自變量，進行 logistic 回歸。亦有研究假設甲基化數據不依賴數據分布形式，可采用非參數回歸模型（如：分位數回歸）。
2.1.4 聚類
Yip 等[17]提出空間聚類（spatial clustering method，SCM）的方法來尋找基因組中與疾病有關的候選 DMR。利用位點間的距離信息，在關聯研究中提高檢驗效能，篩選 DMR。
2.1.5 貝葉斯方法
Teng 等[18]構建了 5 種經驗貝葉斯模型比較不同方法分析差異甲基化效果的好壞。絕大部分經驗貝葉斯模型基于 gamma 分布或對數正態分布的先驗假設，借助 E-M 算法進行參數估計，進而尋找出差異甲基化位點。面對超高維甲基化數據，貝葉斯方法具有符合生物學假設、精度高、檢驗效能高等特點。針對不同分組 CpG 位點甲基化水平方差不齊和樣本量較小的問題，Feng 等[19]提出了貝葉斯分層模型的方法來解決。也可采用 Beta-二項分布分層模型來解釋個體內和個體間兩方面的變異。
差異甲基化分析的結果可用火山圖（volcano plot）和瀑布圖（bar plot）等展示，見圖 2 和圖 3。
圖2
火山圖展示組間差異甲基化分析結果
圖中每個點代表一個甲基化位點；橫軸表示差異甲基化程度，即：fold change；縱軸為差異甲基化統計學假設檢驗的
圖3
瀑布圖
該圖表示甲基化位點在兩組間的差異性情況。橫軸為甲基化位點的差異大小，縱軸為各個甲基化位點。SAS JMP Genomics 10 軟件實現：圖形→圖形生成器。
2.2 結局關聯甲基化位點的識別策略
類似于 GWAS 的思想，甲基化關聯性分析旨在全基因組或部分區域尋找與結局相關的甲基化位點，全面揭示疾病發生、發展與治療相關的表觀遺傳標志物。常見統計方法有：① 單位點的關聯性分析：最基本的 EWAS 關聯性分析是評估每個 CpG 位點與研究結局的關聯性強弱。其中，Cox 比例風險模型是生存分析的常用方法，logistic 回歸模型則是關聯疾病發生風險或鑒別診斷的常用模型。其檢驗結果通常使用 Bonferroni 等方法進行多重比較校正。② CpG 集合的關聯性分析：核函數關聯檢驗是有監督的、靈活的、計算效率較高的回歸模型，采用方差成分檢驗，來提高檢驗效能。自由度的調整取決于 CpG 位點間的相關性結構，若 CpG 位點間呈現高度相關，檢驗的自由度就會大大減小[20]。主要用于識別某一區域內的遺傳變異（包括常見變異和罕見變異）與連續性或分類結局變量間的關聯，同時又便于協變量的調整[21]。GSVA（Gene Set Variation Analysis）是一種非參數的無監督分析方法，主要用于評估基于基因集的富集結果，多用于評估樣本在不同生物學通路上的富集程度[22]。ssGSEA（single sample Gene Set Enrichment Analysis）是 GSEA 方法的擴展，可以計算出每個樣本的富集水平[23]。③ 位點間交互作用分析：CpG 位點間的相關結構較為復雜且不固定，變量間還可能存在一階或高階交互作用[20]，交互作用的檢測方法有叉生分析、分層分析、回歸模型。回歸模型是分析交互作用的常用模型，回歸中交互項回歸系數的估計是以存在相乘交互作用為基礎，在分析交互作用時引入一個相乘項進入回歸模型，通過相乘項回歸系數的估計來判斷交互作用是否有統計學意義及其作用大小。關聯性分析結果可用曼哈頓圖（Manhattan plot）展示，見圖 4。
圖4
曼哈頓圖
橫軸表示甲基化位點在染色體上的位置，縱軸表示每個位點關聯性分析假設檢驗
2.3 基于甲基化數據的結局預測
對于有預測價值的 CpG 位點，為預測疾病的發生風險或預后，對臨床輔助治療、靶向藥物研發具有一定的指導意義。此時，需建立預測模型進行分析，在構建預測模型時，通常遵循著數據降維—模型構建的流程。
數據降維：① 在變量降維篩選中，懲罰方法是使用較為廣泛的一類方法，通過對未知參數的值進行壓縮，同時實現變量選擇和參數估計，具有降低估計偏差、提高預測精度和模型可解釋性的優點。在建模過程中，為防止過度擬合，也可添加懲罰項，引入正則化參數。一個基因或一個區域內含多個 CpG 位點，這些位點通常存在較強的、復雜的相關性結構，其中只有一部分是致病位點。根據這一特點，Sun 等[24]提出懲罰 logistic 回歸模型，在考慮 CpG 復雜相關性結構的前提下，篩選出和疾病相關的獨立 CpG 位點。現有的主流懲罰模型還有 Lasso[25]、Elastic net[26]等。② 主成分分析法（principle component analysis，PCA）目前已被廣泛用于變量間存在相關性的多元分析中。PCA 通常作為一種降維的手段。鑒于 PCA 是一種無監督的降維方法，沒有考慮結局變量與自變量的關系，所以無法保證主成分與結局變量是否存在真正的相關。基于此，有學者提出了有監督的主成分分析法（supervised principle component analysis，SPCA）[27]，先根據自變量與結局變量相關性進行排序，篩選相關性最高的幾個自變量提取主成分，然后評估潛在變量與結局變量間的統計學意義。核函數主成分分析法（kernel principle component analysis，KPCA）是一種非線性 PCA，能處理自變量與因變量之間的非線性關系，近年來在機器學習領域被廣泛運用[28]。③ Leo Breiman 于 2001 年提出隨機森林（random forest，RF）的算法[29]，用于處理位點數多而樣本量少的基因組數據，對高維數據進行降維，利用 RF 算法獲得數據集中每個變量的重要性評分排序，達到初篩變量的目的。其優點在于通過重要性評分排序進行了變量篩選，可大大減少噪音變量，消除冗余，因而可用來分析變量數遠大于樣本數的資料。張秋伊等[30]將 RF 算法運用于類風濕性關節炎病例對照研究的高維甲基化數據的分析，結果表明 RF 能綜合位點信息，有效剔除噪音變量，提高檢驗效能。該方法適用于當變量數遠大于樣本量的情況，可以提高檢驗效能。
此外，基于樹的方法，RF 的強化版算法—Ranger，能夠固定模型中的協變量，在實際分析中具有分析優勢。Wei 等[31]通過 Ranger 算法獲得候選 CpG 位點，再繼續進行傳統 Cox 回歸，得到與早期非小細胞肺癌存活相關的 CpG 位點。
模型構建：① 回歸類參數模型：可以對于篩選后的位點建立 logistic 或 Cox 回歸模型，用 OR 或 HR 估計其對疾病風險大小，多位點的綜合效應可以用風險評分（score）進入模型，采用 估計其對疾病風險的總效應[32]。如果納入變量較多時直接構建參數預測模型，則十分容易導致過度擬合現象（overfitting problems）。故懲罰回歸類方法應運而生，目前已被較多地應用于變量篩選和模型構建，如 LASSO 和彈性網絡（elastic net）等。其通過壓縮回歸系數，將關聯性弱的變量系數壓縮為 0，從而保留了關聯性強的變量[33]。② 基于機器學習的預測分類模型：近些年，機器學習類方法由于模型預測表現較好，受到了極為廣泛的關注。如支持向量機（support vector machine，SVM），已被廣泛用于分類和回歸問題，是一種高效的學習模型。它基于簡單而強大的思想，是一種非常流行的分類方法[34]。RF 是可以用來處理分類、回歸等問題[35]。它有很多優點：具有極高的準確率；隨機性的引入，不容易過度擬合，具有較好的抗噪能力；能處理很高維度的數據，并且無需做特征選擇。K 近鄰作為一種有監督分類算法，是一種基于距離的分類器，在分類方面已經得到了廣泛的研究和討論。其算法主要思想是根據距離相近的鄰居類別，來判斷自己所屬類別。K 近鄰思想簡單，易于理解，易于實現，無需估計參數，無需訓練。樸素貝葉斯（Na?ve Bayes）也是經典的機器學習算法之一，是為數不多的基于概率論的分類算法[36]。其是一個基于貝葉斯定理的比較簡單的概率分類器，原理比較簡單，學習和預測的效率都很高。
建立預后預測模型后，常用生存曲線（Kaplan-Meier curve）、森林圖（forest plot）來對比不同組別生存差異，見圖 5 和圖 6。
圖5
生存曲線
該圖表示了一個甲基化位點隨著隨訪時間的增加發生終點事件概率變化情況。為了比較不同甲基化水平的生存時間差異，通常把該位點進行分組。例如三分類之后，每組的生存曲線如圖所示。SAS JMP Genomics 10 軟件實現：分析→可靠性和生存→生存。
圖6
森林圖
該圖展示甲基化位點在不同亞組中分層分析的結果，藍點和線表示效應值（如風險比）及 95% 可信區間。SAS JMP Genomics 10 軟件實現：圖形→圖形生成器，或插件“Forest Plot”。
3 討論
本文從 EWAS 的數據質控、分析及可視化展開，簡要介紹了甲基化數據規范質控流程和數據分析常用方法，并展示了可視化結果的常見圖形，各部分互相聯系，不可或缺。鑒于目前還未有系統介紹 DNA 甲基化數據處理流程分析策略的文章，本文可為類似研究提供一個參考。對于 EWAS 的關聯性分析和預測模型構建，通常需要進行獨立的外部驗證，確保結果的可重復性、穩健性。在當今的研究中，隨著研究隊列的逐漸增多，對多階段、大樣本數據的外部驗證已成為趨勢[37, 38]。
隨著生物技術的進步，DNA 甲基化的測量手段也在不斷升級。如最近推出的 Infinium Methylation EPIC BeadChip（簡稱 850k 芯片）、Whole Genome Bisulfite Sequencing（甲基化測序）等手段能夠大大提高甲基化位點的識別數量，費用也在逐漸降低，在將來的研究中或許將成為主流。此外，甲基化與基因多態性間具有緊密的聯系，稱為 meQTL；甲基化對基因表達也有著調控，例如可讓基因沉默。如何將基因多態性、DNA 甲基化、基因表達等多種組學的信息整合起來，這將是 GWAS 和 EWAS 統計分析需要進一步探討的問題[39]。
本文所介紹的這些統計分析方法都有各自的適用條件，在實際應用研究中需要結合研究目的、疾病特征、樣本量等具體因素選擇最佳分析方法。
致謝：感謝張瑛、管美建、包文駿博士對 SAS JMP Genomics 軟件可視化內容所提供的建議和幫助。
聲明：所有作者均聲明本研究不存在任何利益沖突。
