一種分析全基因組上位性的新方法

李澤軍+陳敏+曾利軍

摘 要：傳統(tǒng)基于單位點的全基因組關聯(lián)研究存在重復性低、難以解釋性等缺陷，而采用基于機器學習的上位性分析中面臨計算復雜度高、預測準確度不足等問題.本文提出一種分析全基因組上位性的新方法，該方法采用二階段框架的上位性分析方法，它包含特征過濾階段以及上位性組合優(yōu)化階段，在特征過濾階段提出了多準則融合策略，從多個不同角度評價遺傳變異位點，以保證易感的弱效位點能被保留，然后采用多準測排序融合策略剔除與疾病狀態(tài)關聯(lián)程度低的遺傳變異，進一步在上位性組合優(yōu)化階段采用貪婪算法啟發(fā)式地搜索組合空間，以降低時間復雜度，最后采用支持向量機作為上位性評價模型.實驗中采用不同的連鎖不平衡參數(shù)與經(jīng)典算法SNPruler與ACO的性能進行對比，實驗結果表明：本文方法能有效保留弱效位點，一定程度上提高了疾病預測的正確度.

關鍵詞：全基因組關聯(lián)研究；上位性；復雜疾?。恢悄苡嬎?/p>

中圖分類號：TP39 文獻標識碼：A

一種分析全基因組上位性的新方法

Abstract：Traditional units of genome-wide association studies have serious defects such as low repeatability， difficulty to interpret， and epistasis analysis based on machine learning has troubles such as high computational complexity and insufficient prediction accuracy. This paper presented a new approach for the analysis of genome-wide epistatic. This method uses the framework of two-phase epistatic analysis method. It includes a filtering stage and an epistatic combinatorial optimization stage. The characteristics of the filtering stage presents a multicriteria fusion strategy for the evaluation of genetic loci from multiple perspectives to ensure that the weak effect of susceptibility loci can be retained， and then， this method uses the multiple criteria sorting fusion strategy to eliminate the low degree of genetic variation associated with disease states. Epistatic combinatorial optimization phase uses the greedy algorithm combination of heuristic search space in order to reduce the time complexity. Finally， a support vector machine was used as the epistatic evaluation model. Experiments with different parameters of linkage disequilibrium SNPruler with classical algorithms were compared with the performance of the ACO， and the experiment results show that the method can effectively keep weak effect locus and improve disease forecasting accuracy considerably.

Key words： GWAS（Genome-Wide Association Study）；epistasis； complex diseases；intelligent computing

復雜疾病如癌癥等嚴重威脅著人類的健康，它的形成和發(fā)展通常是由多種基因變異所導致，因此對不同患者采用相同的治療措施可能產(chǎn)生不同的治療反應.目前，腫瘤等復雜疾病通常主要采用兩種高通量的數(shù)據(jù)，一種是基因表達譜數(shù)據(jù)，另一種是全基因組單核苷酸多態(tài)性數(shù)據(jù)（SNP）.為提高疾病分析的效率，數(shù)據(jù)挖掘、機器學習等方法被廣泛用于復雜疾病分析中.令人驚訝地是機器學習方法在基因表達譜分析中應用非常廣泛，而在SNP數(shù)據(jù)中的應用則較為少見[1].研究證實，利用機器學習方法提取基于樣本狀態(tài)密切相關的特征基因，然后以此構建預測模型，其性能優(yōu)于僅使用先驗信息證實的候選位點或者通過顯著性檢驗的候選位點 [2-4].然而從高維的全基因組SNP數(shù)據(jù)或者下一代測序數(shù)據(jù)中提取具有預測效果的信號仍然是一種挑戰(zhàn)，所面臨的實驗或者計算難題與傳統(tǒng)的基因表達譜數(shù)據(jù)也有差別[5-6].同時，傳統(tǒng)的基于單位點的復雜疾病分析方法忽略了位點之間的上位性相互作用，導致關聯(lián)研究難以解釋且預測準確度低等不足.為了構建更準確可靠的基于全基因組遺傳變異的復雜表達預測模型，有必要將位點之間上位性效應融入疾病預測模型中[7-8].

當前已經(jīng)有一些基于機器學習的上位性分析方法，比如將復雜性狀（如血壓、身高等）的預測看作是回歸問題或者是疾病狀態(tài)預測問題，根據(jù)學習模型的準確率以評價所選擇的特征屬性（易感位點組合）與疾病狀態(tài)相關的程度[9-10].以較早的機器學習應用為例，Wei等人[11] 在2009年利用支持向量機及L2范數(shù)正則邏輯回歸構造了一個具有較高的預測性能的風險模型，并選擇500個左右的SNP用以預測1型糖尿?。╰ype 1 diabetes T1D）.與之相比，僅考慮少數(shù)已知的T1D易感位點的上位性，使得預測性能大大降低.該方法采用了5折交叉驗證，實驗結果顯示該方法AUC值在0.9左右.但是，簡單的交叉驗證使得兩個特征選擇階段之間存在信息泄露，從而導致過于樂觀的驗證結果.2013年Wei等人再次從15個歐洲國家獲取更大規(guī)模的樣本數(shù)據(jù)集（大于10，000個體），然后對克羅恩?。–D）和潰瘍性結腸炎（UC）進行風險預測[12].該研究利用了定制的具有更高分辨率的SNP分型芯片，可以獲取常見變異以及第一個階段GWAS研究中所忽略的稀罕變異，然后采用相對寬松的閾值（P<10-4）預測選定的10 000個左右SNP，緊接著采用L1范數(shù)正則化邏輯回歸對稀疏遺傳風險建模.該研究的結果與以往一些重要研究結果[13-15]相互印證，從而說明風險預測性能是與樣本規(guī)模、稀罕變異以及機器學習模型密切相關的.

為了解釋身高性狀的遺傳缺失性，Yang等人[16]改進了兩階段框架并且采用簡單線性回歸模型對294 831個變異位點進行分析.該研究采用的是全基因組預測方法，并沒有采用任何特征選擇方法，實驗中對4 000個歐洲后裔人群的身高表型變異進行分析，結果表明，所識別的易感位點能解釋其中45%的表型變異.Makowsky[17]等在歐洲人群中進行全基因組預測模型訓練，然而在完全獨立數(shù)據(jù)集上采用10折交叉驗證R2值大幅下降.這兩個研究進一步闡釋了缺乏特征選擇或者有效模型驗證將產(chǎn)生過度擬合現(xiàn)象.吳蓉暉等[18]利用多種統(tǒng)計準則綜合評價每個遺傳標記，然后利用蟻群算法搜索上位性組合空間，盡管該方法一定程度提高了疾病分類的準確度，但是它存在兩點不足：1）多種準則之間的互補性低，準則之間存在重疊，因此可能導致準則偏差，并且這些準則僅考慮致病等位基因在不同類樣本中的分布差異，而忽略了類內的穩(wěn)定性；2）蟻群算法具有一定隨機性，在不同的運行環(huán)境下得到一致的病組合可能存在差別，因此降低了算法魯棒性.

對以上研究中存在的不足，本文提出一種基于二階段框架的上位性識別方法，它包含特征過濾階段以及組合優(yōu)化階段.在特征過濾階段中，本文提出了利用多準則融合策略從多個互補角度更準確地評價每個位點，以避免有效的弱效位點被剔除；在組合優(yōu)化階段，為了尋找與疾病狀態(tài)關聯(lián)程度最大的多位點上位性組合，本文采用貪婪算法啟發(fā)式地搜索高階組合空間.為驗證易感上位性組合的分類準確度，本文采用留一交叉法評估學習模型的性能.多組數(shù)據(jù)集上實驗結果表明了本文方法的優(yōu)越性.

1 二階段框架的上位性分析

1.1 特征過濾階段

特征過濾階段的主要目的是刪除大量位點中的噪聲位點以及冗余位點，以顯著減少SNP的數(shù)量，從而使得后續(xù)上位性分析階段的組合爆炸現(xiàn)象得到一定程度的緩解.傳統(tǒng)的過濾準則都是采用某一種度量方法如統(tǒng)計檢驗法從單個位點在患病對照組中的分布差異進行評價，導致一些真實易感的弱效位點被剔除，從而使得后續(xù)上位性分析的準確度不足.本文從多個角度同時考察每個位點，以更準確地保留弱效位點.

1.1.1 多個單位點評價準則

1）互信息準則

信息熵理論被廣泛用于SNP數(shù)據(jù)分析中[19-20]，熵作為一種平均自信息度量方法，可以定量度量信息源中信息量，可以描述隨機變量中不確定性程度.在數(shù)據(jù)集噪聲測量中，熵可以直接測量數(shù)據(jù)集的冗余程度或者是噪聲信息含量，其中冗余程度是指數(shù)據(jù)集中各個元素之間的相互依賴程度.假設X表示一個隨機變量空間，可以表示為X={x1，x2，x3…xn}，p（X=x）表示值x出現(xiàn)的頻率，因此隨機變量X信息熵可以表示為：

本文利用互信息公式度量遺傳變異與疾病之間的依賴程度，如式（2）所示.

式中：H（X）為變量X的熵，從該式看它描述變量取值不同時對于熵的影響，也即不確定性程度.通常某種研究背景下，變量值的分布呈現(xiàn)某種概率分布如正態(tài)分布等.信息熵作為一種不確定性度量方法，其值越小則表明隨機變量的不確定越小，即確定性越大.X為不同遺傳位點，Y為疾病狀態(tài)標簽變量；p（xi）為位點上不同等位基因的分布頻率；H（X）為位點的熵，因此式（2）表示單個位點與疾病狀態(tài)之間的關聯(lián)程度，如果關聯(lián)程度越大則表明該弱效位點是易感位點可能性越大.

2）頻率差異準則

易感遺傳位點上不同等位基因與疾病狀態(tài)之間存在明顯的關聯(lián)，當某個樣本攜帶有易感位點上的致病等位基因時，那么該樣本則患病的風險更大.因此，某個等位基因在患病組和對照組中頻率分布差異非常顯著，那么該等位基因更有可能與疾病相關，那么該位點即為易感位點.頻率差異準則如式（3）所示：

式中：Fcontrol表示等位基因在對照組中的頻率；Fcase表示在患病組中的頻率；m表示樣本的數(shù)量.可以發(fā)現(xiàn)FS∈[0，1].

3）類間類內平方和比準則

類間差異表示保證位點在不同類樣本中的分布差異最大，而類內一致表示位點上等位基因在同組樣本中的分布差異應該保持最小化.當某個遺傳位點上等位基因在不同類中的差異越大，并且在相同類的變化越小，那么表明該位點越可能與疾病分類相關，因此打分越高，如式（4）所示：

1.1.2 多準則融合策略

每個遺傳位點在經(jīng)過多個準則評價度量后，需要將不同尺度的打分值進行融合，以確定遺傳位點的不同性能.本文首先利用不同準則對每個位點進行打分排序，然后將不同準則的打分排名融合，如圖1所示.

該階段分析中，首先各準則基于各自度量的含義基礎上，按每個位點成為致病位點的可能性程度進行排序，比如互信息準則中互信息值越大，則該位點對于疾病狀態(tài)的影響越大，因此其排名越高.多準則融合過程中，可以對不同準則的排名順序作加權融合，然后選擇總排名靠前的位點構成候選易感位點.由于本文中，對各種準則沒有任何先驗信息，因此將不同準則看作是等權重的.

1.2 上位性分析階段

基于單個位點的遺傳特征評價方法的計算成本低， 因此更適合用于全基因組上的SNP數(shù)據(jù)分析過濾階段.同時，該過濾準則不與學習模型相嵌套，因此進一步降低時間復雜度.該過濾準則的最大弊端在于忽略了位點之間的相互作用即上位性，因此，本文在第二階段采用學習模型進一步考慮多個位點的相互作用.

1.2.1 基于支持向量機的疾病狀態(tài)預測

支持向量機（Support vector machine， SVM）是一種基于有監(jiān)督的學習預測模型，該模型最大優(yōu)勢在于維度無關性，即屬性的維度不影響模型的復雜度，因此特別適用于高維的全基因組SNP數(shù)據(jù).概括而言，SVM目標是構造最優(yōu)超平面保證不同類樣本間隔最大化，同時保證泛化誤差盡可能小.

將候選的上位性位點作為特征屬性組合，利用留一交叉法驗證候選上位性組合的疾病預測性能.留一交叉法將每個樣本單獨作為一份測試集，其余樣本作為訓練樣本，利用訓練樣本上特征屬性組合訓練支持向量機，然后將所預測的測試樣本類別號與其真實的類別號進行比較，反復迭代訓練、測試N次，然后計算平均錯誤率.

1.2.2 貪婪算法

從候選易感SNP集合中選擇一個最優(yōu)的上位性SNP組合，使其數(shù)量最少并且疾病預測準確度最高，該問題是一個NP問題，尤其當候選易感SNP集合中SNP數(shù)量仍然較大時，難以搜索到最優(yōu)解.因此，本文采用貪婪策略尋找近優(yōu)解.貪婪算法見表1（又稱貪心算法），首先將待求解的問題劃分為若干個子問題，然后分別得到子問題的最優(yōu)解，最后將所有子問題的解合成原問題的解.在決策過程中并不是從全局考慮每個階段的策略，而僅僅從當前情況下選取最有利的策略，可以看出，該算法并不能保證求解的全局最優(yōu)，但是當解集空間過大時，采用該算法尋找近優(yōu)解也是一種較好的替代策略.

假設候選易感SNP集合I0中含有s個SNP，那么在第1次迭代情況下，要從s個SNP中選擇一個位點剔除，找出一個包含有s-1個位點的子集I1，該子集的疾病預測準確度增加的最多，該次迭代的計算復雜度為O（s），然后在I1的子集中選擇一個位點滿足式（5），依次迭代直到滿足退出條件.該優(yōu)化表達式如下：

Max：Acct-1-Acct. （5）

Acct表示第t次迭代過程中所對應的子集具有的最大預測準確度.貪婪算法的結束條件是，從當前子集中刪除任何一個位點，預測準確度都會降低，即子集中已經(jīng)沒有冗余位點.

從算法流程來看，當最終上位性組合中包含的SNP數(shù)目為n時，該貪婪算法的時間復雜度為O（n*s*T），其中T表示預測過程所消耗的時間.貪婪算法的流程圖如圖2所示.

2 數(shù)據(jù)集及評價指標

為驗證本文多準則策略的有效性，將利用模擬數(shù)據(jù)集評價本文方法.模擬數(shù)據(jù)集中包含2 000個樣本（1 000個病例樣本和1 000個對照樣本），每個數(shù)據(jù)集包含100個標記（其中2個致病標記和98個非致病標記）.首先，生成2個致病標記的基因型分布，要設置4個參數(shù)：疾病外顯率p（D），邊際效應λ，連鎖不平衡r2，次等位基因頻率MAF.通過這4個參數(shù)的值計算得到基因效應θ和基線效應α的值.再根據(jù)基因效應θ和基線效應α的值生成相應的仿真數(shù)據(jù)集.本文將這4個參數(shù)分別設置為p（D）=0.1；λ=0.3；r2=0.5， 0.7， 1；MAF=0.05， 0.1， 0.2和0.5構造模擬數(shù)據(jù)集，數(shù)據(jù)集分別標記為sim1，sim2和sim3.

本文將采用預測準確度指標.如果某個上位性組合是真實的致病上位性組合，那么利用它理論上可以很好地預測個體的患病狀態(tài).預測準確度如式（6）所示：

3 實驗結果及分析

圖3～圖5比較了ACO算法[18]和SNPRuler[21]這幾種方法的預測準確度.

從圖3～圖5可以發(fā)現(xiàn)，隨著上位性位點數(shù)目的增加，本文方法的預測準確度也在增加，這是因為上位性位點增加，能引入更多與疾病分類狀態(tài)相關的易感位點.同時，由于采用了貪婪策略以保證每輪迭代中預測準確度至少保持不變，因此，以上圖中本文方法所對應的曲線基本上是遞增的.圖中也可以看出，其余兩種方法的準確度不斷增加，這是因為兩種方法都能從候選子集中選擇與疾病狀態(tài)相關的易感位點.

上位性分析的時間復雜度為Cmn，其中n表示數(shù)據(jù)集中包含的位點數(shù)量，m則表示上位性的階數(shù)，因此搜索全局最優(yōu)解的成本過于昂貴.盡管本方法中貪婪算法難以保證搜索到全局最優(yōu)解，但該方法適用于在較大數(shù)據(jù)集上搜索高階上位性，一定程度提高了上位性分析策略的適用性.

從圖3～圖5中發(fā)現(xiàn)，本文方法的預測準確度總體上看來是高于其它兩個方法的.這表明本文方法能選出與疾病更為相關的易感上位性組合.其它兩種方法SNPruler與ACO的性能存在一定波動性，在圖3和圖4中，后兩種方法的準確度比較接近，但是在圖5中，ACO方法優(yōu)于SNPruler方法.

將圖3，圖4和圖5進行比較發(fā)現(xiàn)，隨著連鎖平衡值的增加.3種方法的準確度都有所增加，以本文方法為例，圖3中本文方法的評價準確度為0.72，圖4的平均準確度為0.83，而圖5中的平均準確度為0.86.進一步分析該現(xiàn)象可以推測，當致病位點與其它標記之間的連鎖不平衡性增加，那么間接表明該致病位點更可能被其它標記所描述，因此，致病位點所包含的信息也能被學習模型所利用，從而提高了預測準確度.

4 結 論

傳統(tǒng)基于單位點的全基因組關聯(lián)研究具有計算簡便性，但是過度簡化了復雜疾病的致病模型，從而導致研究結果的重復性低.本文首次提出了一種多準則的上位性分析方法，它從多個角度互補地評價了遺傳位點，從而避免了真實的弱效位點被剔除.在上位性組合優(yōu)化階段，本文采用貪婪算法作為優(yōu)化策略，盡管該策略仍然不能全局最優(yōu)性，但是該算法的魯棒性高，避免了研究結果的隨機性，從而更適用于臨床應用.

參考文獻

[1] KRUPPA J， ZIEGLER A， KONIG I R. Risk estimation and risk prediction using machine-learning methods[J]. Human Genetics， 2012， 131（10）： 1639-1654.

[2] PAHIKALA T， OKSER S， Airola A， et al. Wrapper-based selection of genetic features in genome-wide association studies through fast matrix operations[J]. Algorithm Mol Biol，2012， 7（1）：11.

[3] OKSER S， LEHTIMAKI T， ELO L L， et al.Genetic variants and their interactions in the prediction of increased Pre-clinical carotid atherosclerosis[J]. The Cardiovascular Risk in Young Finns Study， PLoS Genet，2010， 6（9）：e1001146.

[4] KOOPERBERG C， LEBLANC M， OBENCHAIN V. Risk prediction using genome-wide association studies[J]. Genet Epidemiol， 2010， 34（7）：643-652.

[5] KRAFT P， WACHOLDER S， CORNELIS M C， et al. Beyond odds ratios： communicating disease risk based on genetic profiles[J]. Nat Rev Genet ，2009，10（4）：264-269.

[6] ASHLEY E A， BUTTE A J， WHEELER M T， et al. Clinical assessment incorporating a personal genome[J]. Lancet， 2010，375（9725）：1525-1535.

[7] MANOLIO T A. Bringing genome-wide association findings into clinical use[J]. Nat Rev Genet， 2013，14（8）：549-558.

[8] GIBSON G. Hints of hidden heritability in GWAS[J]. Nat Genet，2010， 42（8）：558-560.

[9] YANG J， BENYAMIN B， MCE VOY B P， et al. Common SNPs explain a large proportion of the heritability for human height[J]. Nat Genet，2010，42（11）：565-569.

[10]MAKOWSKY R， PAJEWSKI N M， KLIMENTIDIS YC，et al. Beyond missing heritability： prediction of complex traits[J]. PLoS Genet ，2011，7：e1002051.

[11]WEI Z， WANG K， QU H Q， et al. From disease association to risk assessment： an optimistic view from genome-wide association studies on type 1 diabetes[J]. PLoS Genet， 2009，5： e1000678.

[12]WEI Z， WANG W， BRADFIELD J， et al. Large sample size，wide variant spectrum， and advanced machine-learning technique boost risk prediction for inflammatory bowel disease[J]. Am J Hum Genetics，2013， 92（6）：1008-1012.

[13]CHATTERJEE N， WHEELER B， SAMPSONJ， et al. Projecting the performance of risk prediction based on polygenic analyses of genome-wide association studies s[J]. Nat Genet，2013， 45（4）：400-405.

[14]DUDBRIDGE F. Power and predictive accuracy of polygenic risk scores[J]. PLoS Genet ，2013，9： e1003348.

[15]DO C B， HINDS D A， FRANCKE U， et al. Comparison of family history and SNPs for predicting risk of complex disease[J]. PLoS Genet，8： e，2012，1002973.

[16]YANG J， BENYAMIN B， MCEVOY B P， et al. Common SNPs explain a large proportion of the heritability for human height[J]. Nat Genet， 2010，42（7）：565-569.

[17]MAKOWSKY R， PAJEWSKI N M， KLIMENTIDIS Y C， et al. Beyond missing heritability： prediction of complex traits[J]. PLoS Genet ，2011，7：e1002051.

[18]吳蓉暉， 盧友敏. 基于蟻群算法的復雜疾病上位性分析方法[J]. 湖南大學學報：自然科學版， 2014，42（8）： 125-131.

WU Rong-hui ， LU You-min. An epistasis analysis method of complex diseasesBased on ant colony algorithm[J]. Journal of Hunan University： Natural Sciences， 2014，42（8）： 125-131.（In Chinese）

[19]LIU Z， LIN S. Multiocus LD measure and tagging SNP selection with generalized mutual information[J]. Genetic Epidemiology， 2005， 29（4）： 353-364.

[20]LI X， LIAO B， ZHU W，et al.Informative SNPs selection based on two-locus and multilocus linkage disequilibrium： criteria of max-correlation and min-redundancy[J]. IEEE/ACM Trans Comput Biol Bioinform， 2013， 10（3）： 688-695.

[21]XIANG Wan， CAN Yang，QIANG Yang，et al. Predictive rule inference for epistatic interaction detection in genome-wide association studies[J]. Bioinformatics， 2010，26 （1）：30-37.