基于U統計量和集成學習的基因互作檢測方法

郭穎婕 劉曉燕 吳辰熙 郭茂祖,3 李 傲

1(哈爾濱工業大學計算機科學與技術學院 哈爾濱 150001)2(Rutgers大學數學系 美國新澤西洲皮斯卡特維 08854)3 (建筑大數據智能處理方法研究北京市重點實驗室(北京建筑大學) 北京 100044) (yjguo0625@gmail.com)

復雜疾病如癌癥、糖尿病等，因其發病率高、治愈率低的問題，長期危害人類健康.對此類疾病遺傳變異致病機理的研究是生命科學界的重大課題.全基因組關聯研究(Genome-Wide Association Studies，GWAS)作為一種發現遺傳變異——單核苷酸多態性位點(single nucleotide polymorphism, SNP)的新興研究模式，自2005年《Science》雜志發表首例關于GWAS的研究成果——老年性黃斑變性以來，GWAS的研究已經超過1 800項，并發表了關于237種不同疾病研究成果，共計24 000余個SNP被證實與這些疾病顯著相關[1].但對于大多數復雜遺傳相關疾病而言，這些成果僅僅闡釋了遺傳因素很小的一部分.例如，已發現的與老年性黃斑變性、高膽固醇、阿茲海默癥疾病相關的致病風險基因，分別只占相應遺傳致病因素的50%，25%及23%.這一現象被稱為“失蹤遺傳(the missing heritability)”，是目前GWAS研究的熱點.該現象的發生是由于復雜疾病由多對微效基因與環境因素共同作用所致，即基因之間存在復雜的相互作用關系.而傳統的GWAS主要基于單位點模型，獨立檢測單個SNP與給定表型之間的關系，從而導致一些存在基因互作的性狀相關基因無法被單獨檢測出來.因此，識別與分析基因相互作用，成為解決“失蹤遺傳”的可行方案之一.盡管利用基因相互作用不能完全解決“失蹤遺傳”，但此類研究依然可以通過輔助構建新的基因通路為復雜性狀的遺傳與調控機理提供生物學解釋.

基因相互作用又稱為“上位效應”(epistasis).此概念最早由Bateson在1909年提出[2],后由Fisher進一步給出定義[3]：上位效應是指2個位點之間相互作用形式，使用2個位點建模產生的性狀值與單獨使用每個位點得到的性狀值的和之間產生的差值來衡量上位效應的強度.后續研究中，對上位效應的定義變得更為寬泛，不同的方法對于相互作用也有不同的定義.初期的上位效應研究主要著眼于2個SNP位點間的相互作用關系.在精神分裂癥、阿茲海默癥、乳腺癌等復雜疾病中均有關于2個SNP之間交互作用的報道.因此產生了許多以2個SNP、或者3個SNP乃至多個SNP為研究單位的上位效應檢測方法.統計類的檢測方法通過設計一些表征相互作用強度的統計量，檢測顯著的相互作用關系，例如基于優勢比(odds ratio,OR)的統計量、基于連鎖不平衡(linkage disequilibrium, LD)的統計量以及基于熵的統計量等.另一類方法則采用計算機方法的思想，例如采用降維技術的多因子降維方法[4]、基于樹模型的TEAM(tree-based epistasis associa-tion mapping)方法[5]、通過優化存儲策略加速計算的BOOST(Boolean operation-based screening and testing)方法[6],以及BEAM(Bayesian epistasis asso-ciation mapping)，pRF等.這些基于位點的檢測方法面臨最大的挑戰是維數災難.由于必須考慮所有的SNP或SNP組，成對或者高階的相互作用關系檢測次數隨互作關系階數呈指數級增長，隨之而來的對統計顯著性的校正會導致統計效力的弱化.因此，本文研究以基因為單位，將一個基因中的所有SNP看作一個組來檢測基因之間的相互作用.

基于基因的相互作用研究有以下3個明顯的優勢：

1) 基于基因的方法可以大大減少所需的檢測次數，例如，對于20 000個基因，成對檢測基因之間相互作用是可行的；而對于300萬SNP，成對檢測顯然不現實.

2) 2組基因之間可能存在多對SNP的相互作用，組內的SNP之間也可能存在連鎖不平衡關系，這些同時存在的作用關系會隱性地呈現在以基因為單位的模型中，更利于相互作用的檢測.

3) 基于基因的方法可以更好地利用已有的生物學背景知識，縮小研究范圍.例如可以檢測那些蛋白質互作網絡(protein-protein interaction， PPI)中已經呈現互作關系的蛋白質編碼基因之間的關系，或者某個通路內基因之間的相互作用關系.

目前，在以基因為單位的相互作用研究中，Peng等人[7]在疾病組與對照組中分別對2個基因進行典型相關性分析(canonical correlation analysis, CCA)，并設計統計量CCU(canonical correlation-based U statistic)來度量2個基因在疾病與對照組中相關性指標的差異程度，用于表征相互作用的強度.該方法的局限性在于CCA只能度量2個基因之間的線性關系.Yuan等人[8]和Larson等人[9]針對Peng等人的方法存在的問題，將CCU擴展到KCCU(kernel canonical correlation-based U statistic)，在做典型相關性分析之前，將核函數作用在疾病和對照組中2個基因的數據上，從而增強模型對非線性關系的解釋能力.Li等人[10]提出了一種基于熵的非參數假設檢驗方法GBIGM(gene based informa-tion gain method).通過分析2個基因共同作用時與考慮只有單個基因時熵的變化(即信息增益)，并利用隨機置換類標簽的方式獲得相互作用的顯著性p值.Emily[11]開發了AGGrEGATOr(a gene-based gene-gene interaction)方法，該方法首先計算兩基因間所有SNP對的Wald統計值，并將一組Wald統計值結合成為一個顯著性p值用于度量2個基因之間是否存在相互作用.此前，Ma等人[12]成功地將這一策略用于數量型性狀的基因互作檢測中.

本文的主要貢獻有4個方面:

1) 提出了一種基于U統計量與集成學習結合的假設檢驗方法GBUtrees(gene-based gene-gene interaction detection based on U statistics and tree-based ensemble learners)，將基因相互作用問題轉化為檢驗2組SNP與疾病性狀的對數優勢比(log odds ratio, LOR)之間是否滿足加性模型關系；

2) 利用以樹為基分類器的集成學習方法作為底層學習方法，充分發揮該類方法對非線性關系的建模能力，有效學習2組基因之間可能存在的相互作用關系；

3) 通過多次采樣子樣例集使得預測模型的結果滿足U統計值框架，利用其滿足漸進正態分布的性質設計表征基因相互作用強度的統計量；

4) 利用8種不同相互作用模式疾病模型生成的模擬數據集與真實人類類風濕關節炎數據的實驗結果表明，我們提出的GBUtrees方法可以有效挖掘基因-基因之間不同形式的相互作用關系.

1 基于U統計量的度量集成學習不確定性方法

1.1 基因互作檢測問題描述

本文以基因作為基本研究單位.假設基因G1和G2分別包含p和q個SNP，每個SNP位點只有高頻等位基因(major allele)與低頻等位基因(minor allele)兩種堿基形式，分別記為A和a，則對于人類二倍體而言，每個SNP位點存在3種可能的基因型AA，Aa，aa.數據中使用0,1,2分別對應上述3種基因型.則基因G1可表示為

G1=(s1,s2,…,sp),

(1)

其中si∈{0,1,2}，i=1,2,…,p.

同理，包含q個SNP的基因G2可表示為

G2=(s1,s2,…,sq),

(2)

其中si∈{0,1,2},i=1,2,…,q.

本文方法主要針對疾病-對照數據，因此我們的目標值y∈{0,1}，其中1表示疾病組，0表示對照組.如果對于G1，G2與y的LOR之間存在以下的關系：

(3)

其中F1，F2可以是任意形式的函數.我們則稱G1，G2與y之間滿足加性關系，即G1，G2之間沒有相互作用關系.換言之，當基因G1中變量發生改變時，并不會改變基因G2對性狀的影響.因此，想要判斷2個基因之間是否有相互作用關系，可以通過建立如下假設檢驗來衡量觀測數據是否可以被建模成為加性模型：

H0:

?F1,F2,F(G1,G2)=
F1(G1)+F2(G2),?(G1,G2)∈XTEST,

(4)

H1:

?F1,F2，F(G1,G2)≠
F1(G1)+F2(G2),?(G1,G2)∈XTEST,

(5)

其中，XTEST表示測試集.

使用監督學習中的集成學習方法對目標值進行建模，同時獲得集成學習方法中的統計指標，來完成上述假設檢驗.Mentch等人[13]在2016年提出了基于U統計理論度量隨機森林方法中的不確定性，并于2017年進一步完善該方法[14].本文利用該假設檢驗框架實現對基因互作的檢測.下面將分別介紹U統計理論、集成學習方法以及基于該假設檢驗框架提出的基因互作檢測方法.

1.2 U統計量

θ=Ε(h(Z1,Z2,…,Zk))

(6)

其中,h是一個函數，有k≤n個變量，則參數θ的最小無偏估計如下：

(7)

(8)

1.3 基于樹模型的集成學習

集成學習使用多個分類器來解決同一問題，通過組合多個弱監督的基分類器以期得到一個更為全面的強監督模型.通過合理的選擇并組合基分類器，集成學習可以在顯著提高一個學習系統的泛化能力的同時提高學習器的分類精度.該方法已在多個領域成功應用[16-17].在1.2節U統計量的框架下，任何類型基分類器，只要能夠寫成式(7)的形式，均可以應用到本文的方法中.考慮到基因互作問題中存在嚴重的非線性相互作用，本文選取決策樹模型作為集成學習中的基分類器.假設訓練集有n個樣本，

(9)

其中X=(X1,X2,…,Xd)是一個特征向量，y∈是目標函數值.選取一個小于n的k值，(Xi 1,yi 1),(Xi 2,yi 2),…,(Xi k,yi k)是訓練樣本的一個子樣例集.給定一個待測樣本x*，我們把一個樹模型在該子樣例集上的預測模型記為Tx*.為個子樣例集都訓練一棵樹，那么最終待測樣本x*的預測可表示為

(10)

1.4 基于U統計理論的基因互作檢測方法GBUtrees

針對1.1節中式(4),(5)對應的假設檢驗問題，我們設計如下方法：

首先，我們要構建一個用于作為測試集的測試網格(test grid).測試網格中的數據并非真實數據，其重點在于模擬真實數據中特征的分布.例如本文研究的是2個基因相互作用，所以設計一個二維網格.網格大小為N1×N2，網格中每個位置代表G1,G2的一組取值情況.為了避免測試數據與真實數據重復，我們預先從原始數據中抽取N1+N2個樣本用于構建網格.將前N1個樣本的G1與后N2個樣本的G2分別組合，從而生成N1×N2的網格.

(11)

則有：

(12)

(13)

式(12)表示G1取第i個SNP組合時，遍歷測試網格上所有G2時的預測平均值，同理可得式(13)是G2取第j個SNP組合時，遍歷測試網格上所有G1時的預測平均值.如果G1，G2之間是加性關系，則網格中所有的點(g1i,g2j)都可以表示成Fi,j=fi ·+f·j-μ，其中μ表示測試網格中所有點預測值的真實期望.因此，1.1節中的假設檢驗也可以表示為

H0:Fi,j-fi ·-f·j+μ=0,
?(g1i,g2j)∈XTEST.

(14)

H1:Fi,j-fi ·-f·j+μ≠0,
?(g1i,g2j)∈XTEST.

(15)

令N=N1×N2，我們引入一個維度為N×N的差異矩陣D:

(16)

(17)

實驗中存在一個問題，當測試網格數據樣本量超過50時，方法性能會大幅下降,因此考慮降低測試網格的維度.Johnson-Lindenstrauss定理保證了任何降維方法的精度上下限.影響降維精度的因素主要是維數、數據大小等，與降維方法本身無關.在維數差不是特別大時，用任何方法都可以保證一個較高的精度.本文采用隨機投射(random projec-tion)方法.該方法是機器學習中常見的降維方法.通過隨機生成一個p×q的投射矩陣，并對矩陣中的每一個軸(即每一行)歸一化，之后左乘數據，便可將X∈q降至p.實驗中我們構造1 000個隨機矩陣，降維的維數設為15.詳細方法描述見算法1.

算法1. 基于U統計量的基因互作檢測方法GBUtrees.

輸入：2個基因的基因型數據X=(X1,X2)、樣本的類標簽Y∈[0,1]n;

① 將原始樣本集劃分為訓練集和測試集，測試集樣本個數為N1+N2;

② 計算差異矩陣D;

③ 選取測試網格約簡的維度r=15;

④ 生成M=1 000個隨機映射矩陣R1,R2,…,RM;

⑦ Forj=1,2,…,Ndo

算法1中涉及到的主要參數分為以下2個部分：

2) 構造每一棵樹時，需要考慮每個分裂節點至少包含樣本個數(minsplit)以及樹的最大深度(maxdepth).利用網格搜索(grid search)方式，對每個minsplit與maxdepth組合做5重交叉驗證(5-fold cross-validation)，最終選取minsplit=3，maxdepth=4.

2 實驗與結果

本文實驗分為仿真模擬數據與真實數據2部分.

2.1 仿真模擬實驗

2.1.1 仿真數據來源

仿真實驗選用經典的基于重采樣仿真的SNP遺傳數據生成軟件gs2.0.該軟件以單體型數據為模板，通過設定疾病優勢比(OR)、人群患病率(population prevalence)以及樣本大小，可以高效產生大量疾病-對照SNP仿真數據.

本文使用的模板數據是Hapmap 高加索人群(CEU)數據，使用NCBI build37作為參照基因組進行注釋.該數據共有90個單體型，以1 000 Genome數據為參考基因組，通過使用SHAPEIT[18]與IMPUTE2，對缺失的位點進行基因型填充.隨后分別在一號染色體與二號染色體上各選取了一個基因[11-12]，一號染色體上基因選取了14個tag SNP，二號染色體上基因選取了10個tag SNP.2個基因內部連鎖不平衡(LD)模式見圖1.選取SNP的詳細信息見表1.

Fig. 1 LD patterns of two empirical loci used in simulation studies圖1 仿真數據2個基因模板的連鎖不平衡模式圖

IndexSNP name:positionLocus1(chr1)Locus2(chr2)1rs11589332:120273204rs12712643:396263882rs512854:120275063rs11680220:396306573rs539426:120275647rs1558854:396395724rs593911:120288463rs7598583:396472215rs536662:120292733rs10779925:396895886rs668156:120293568rs11891871:397284317rs17186233:120297163rs17467001:397354628rs1441010:120301432rs13420425:397358139rs1441009:120307515rs7585512:3973589910rs12563433:120308975rs17532603:3974010211rs3828089:12030927412rs1441008:12031029213rs753425:12031674614rs10737757:120320726

2.1.2 仿真模型

gs2.0共有8種可選的兩位點疾病遺傳風險模型.根據方法不同性能的考量，實驗分為以下2個方面：

1) 為了考察方法第1類錯誤率的表現，設定優勢比OR=1.0，此時生成不存在相互作用.分別生成樣本數為1 000，2 000，3 000，4 000，5 000的樣本集各100個，觀察第1類錯誤率在不同樣本規模時的變化規律.

2) 為了考察方法在不同交互程度時的表現，將優勢比設為1.5，2，2.5，3，3.5，4共6個取值，人群患病率(population prevalence)設為0.1，樣本由2 000疾病樣本、2 000對照樣本組成.8種疾病模型共有48種組合，每種組合生成100個數據集.

2.1.3 實驗結果分析

為了驗證基于U統計框架方法的有效性，本文采用功效(power)與假陽率2個指標.功效即在每種參數設置下，可以成功檢測到模擬數據中插入的相互作用的數據集占100個數據集的比例.GBUtrees與GBIGM[10]，KCCU[8-9]，AGGrEGATOr[11]共4種方法進行了比較.其中，GBIGM與AGGrEGATOr兩種方法沒有參數限制；KCCU方法采用文獻中推薦的參數，設定修剪摺刀法(trimmed jackknife)的修剪比例為ω=0.15.

在假陽性方面，設定顯著性閾值α=0.05，隨著樣本規模增加，GBUtrees犯第1類錯誤的概率均可控制在顯著性閾值以內，可以用于后續的相互作用檢驗.

(18)

則基因型gi的外顯率為

(19)

Fig. 2 Power comparison between GBUtrees and competitive methods under four disease models with different odds ratio圖2 4種交互作用檢測方法在8種疾病模型的不同OR值下的功效比較

由此我們可以得到2個位點基因型組合的外顯率表，見表3.

Table 2 Table of Odds for Recessive-Recessive Model表2 隱性-隱性模型的疾病比率表

Fig. 3 Power comparison between GBUtrees and competitive methods under four disease models with different sample size圖3 4種交互作用檢測方法在4種疾病模型的不同樣本量下的功效比較

Table 3 Penetrance Table for Recessive-Recessive Model表3 隱性-隱性模型的疾病外顯率表

3種對比方法中，AGGrEGATOr在8種疾病模型下均表現出與GBUtrees相似的變化趨勢與性質，而KCCU與GBIGM在該模擬實驗中幾乎無法檢測出相互作用.本文中仿真實驗與GBIGM中模擬實驗[8]的區別在于選取的模板數據不同，不同的模板數據由于內部不同的連鎖不平衡結構會對方法的性能產生影響.

此外，樣本量對方法性能影響的結果如圖3所示.隨著樣本量增加，GBUtrees與AGGrEGATOr性能均呈現單調遞增的趨勢，另外2種方法性能沒有隨樣本數量的增加而發生改變.圖3只展示了4個模型下方法的表現，GBUtrees在剩余疾病模型上均表現出類似的變化趨勢，且優于其他3種方法.

2.2 真實數據實驗

2.2.1 真實數據描述及預處理

為了評估方法在處理真實疾病對照數據時的表現，我們選用了一組類風濕關節炎(rheumatoid arth-ritis， RA)的真實基因型數據.RA是一種遺傳相關的慢性的自身免疫性疾病，持續性的炎癥會影響骨重塑，進一步導致骨壞死.我們使用WTCCC(2007) 數據集，該數據使用Affymetrix GeneChip 500K進行基因型鑒定.為了驗證潛在的基因之間的相互作用，我們參考了KEGG①上的RA 通路hsa05323.該通路共包含90個基因，其中MHCII與V-ATP是2個蛋白質結合體，內部存在許多的相互作用.因此，對這2個復合體各選取了一個代表基因.通過使用NCBI發布的人類基因組序列NCBI build36的注釋文件確定每個基因在染色體上的起始終止位置，并為每個基因的上下游增加了10 kb的緩沖區，我們從RA數據中為最終剩余的48個基因確定其包含的SNP.隨后借助軟件PLINK對數據進行質量控制，共包括3個步驟：1)刪除性別與性染色體不匹配的數據；2)設置SNP缺失率(missing rate)小于10%、低頻等位基因頻率MAF>5%，以及顯著影響表型的缺失；3)刪除不滿足Hardy-Weinberg平衡的SNP.經過質量控制之后，共剩余385個SNP和4 966個樣本，其中4 966個樣本包括1 973個病例樣本以及2 993個對照樣本.

2.2.2 真實數據結果分析

為了避免人群結構、性別成為干擾因素，我們將性別作為變量加入數據.隨后使用軟件GCTA對數據做主成分分析，并將前10個分量加入數據用于校正潛在人群結構帶來的影響.

Table4RankingofSignificantGGIResultsintheWTCCCRheumatoidArthritisDatasetAnalysis

表4WTCCC人類類風濕數據集上顯著相互作用關系在不同方法中的秩次

Gene 1Gene 2RankingGBUtreesAGGrEGATOrKCCUGBIGMAng1Tie211099560260Tie2Flt121108685460Tie2LFA1311217391035MHCIIAng14109534270MHCIITie25858950680MCSFAng163320590MCSFTie271033485495Ang1Flt181029265500Tie2RANK91117720524TGF-betaTie2101037745580CD86CCL290011066650TGF-betaAP12702605260MCSFAng163320590CD86TLR4650450565CXCL1RANKL7705130825CCL20TGF-beta6606300350RANKLAPRIL5407270180RANKTGF-beta340887568IL23LFA17509395530IL1V-ATPase102010360905

GBUtrees檢測出最顯著的相互作用來自血管生成素(Ang1)與Tie2.該相互作用是一對RA致病機制中已被證實的互作關系[19].血管翳是RA病變過程中的一個特征性的病例產物.滑膜血管增生是RA早期關鍵事件，是促進和維持血管翳的重要因素.而RA血管生成是由滑膜組織髓細胞和纖維母細胞(包括血管生成素Ang1)釋放的促血管生成因子驅動和維持的.內皮細胞(EC)的特異性因子Ang1及其酪氨酸激酶受體Tie2在生理性和病理性的血管生成發育中至關重要.Ang1/Tie2在RA滑膜組織中上調表達.有研究表明Ang1/Tie2信號通路介導了腫瘤壞死因子(TNF-)、白介素(IL)以及toll類受體TLR2在RA中促進血管生成的作用.AGGrEGATOr方法發現的第7對是RANKL與APRIL之間的互作[20].研究表明纖維細胞樣滑膜細胞FLS是RA發病機制中的主要效應細胞，而APRIL促進了RANKL的FLS表達.KCCU與GBIGM方法給出了太多顯著性結果，因此沒有對這2種方法的結果進行進一步分析.

3 總 結

檢測基因-基因相互作用的研究在理解人類復雜疾病致病機理方面具有重要意義.本文提出一種基于U統計量與集成學習結合的假設檢驗方法GBUtrees用于基因相互作用檢測.我們將檢測基因之間的相互作用問題轉化為假設檢驗問題，定義零假設為2個基因與疾病標簽的LOR之間滿足加性關系.這一假設對基因之間互作的形式沒有限制，增強了方法對不同類型作用關系的檢測能力.仿真數據實驗結果表明，在可控第1類錯誤率的前提下，GBUtrees在模擬的8種疾病模型上均優于相比較的其他方法，且方法效力隨著作用強度的增大以及樣本量的增大呈現單調遞增.此外，在真實類風濕關節炎數據驗證實驗中，GBUtrees成功復現了已被證實的Ang1/Tie2相互作用.以上結果表明該方法在基因相互作用挖掘中的有效性.

GBUtrees采用集成學習方法作為底層學習方法，增強了方法的泛化能力.同時，以樹模型作為基分類器，發揮其強大的非線性建模能力，大大增加了方法檢出不同類型相互作用的功效.此外，通過設計重采樣的策略，使預測結果具有U統計量的漸進正態性質，從而可以設計用于表征相互作用關系強度的統計量.獲得相互作用關系強度p值對于檢驗結果的復現以及應用到后續的全基因組關聯研究中都具有至關重要的意義.仿真模擬數據以及真實人類疾病數據的實驗結果表明，本文提出的GBUtrees方法在檢測效力方面優于現有的方法，可以有效用于疾病相關基因相互作用的挖掘研究.

GuoYingjie, born in 1987. PhD candidate. Her main research interests include machine learning and bioinformatics.

LiuXiaoyan, born in 1963. PhD. Associaate professor and master supervisor. Her main research interests include bioinformatics and data mining.(liuxiaoyan@hit.edu.cn).

WuChenxi, born in 1988. PhD. Assistant professor. His main research interests include flat surface, three dimension topology.

GuoMaozu, born in 1966. Professor. PhD supervisor. Member of CCF. His main research interests include machine learning, smart city and bioinformatics.

LiAo, born in 1995. Master. His main research interests include machine learning and deep learning.