基于通路分析的遺傳交互網絡方法：BGTA方法及應用*

任志強 皮路程 劉桂炎 劉 晴 程金群 郜艷暉△

【提 要】 目的 采用BGTA(backward genotype-trait association)算法結合通路分析策略探索全基因組關聯研究中廣泛的遺傳交互作用。方法 采用GAW19(genetic analysis workshop 19)中無相關人群外顯子測序數據及真實高血壓表型數據，在KEGG(kyoto encyclopedia of genes and genomes)中選擇高血壓相關腎素-血管緊張素-醛固酮系統(renin-angiotensin-aldosterone system，RAAS)通路中的遺傳變異作為候選變異，采用兩階段的 BGTA 算法進行基因交互作用分析并構建交互網絡圖，并與隨機森林聯合決策樹方法結果比較。結果 BGTA兩階段分別篩選出76個(含61個低頻)和56個(含44個低頻)高血壓相關遺傳變異(P<0.10)。logistic回歸驗證有16對無主效應的變異間交互作用(P<0.05)。第一階段隨機森林基于MDG(mean decrease Gini)和MDA(mean decrease accuracy)分別篩選出35個(含0個低頻)和69個(含30個低頻)遺傳變異，第二階段決策樹基于MDG和MDA分別篩選出5個和7個遺傳變異，未發現低頻變異。logistic回歸驗證了7對無主效應交互(P<0.05)。結論 兩階段BGTA在探索RAAS通路遺傳變異交互作用與高血壓關聯時，比隨機森林聯合決策樹方法發現更多無主效應交互作用。將BGTA算法和生物學通路分析方法結合應用于復雜疾病的全基因組關聯研究中，可提高關聯變異的識別能力，為了解復雜疾病的遺傳結構提供線索。

近年來，GWAS(genome-wide association study)和二代測序技術(next-generation sequencing )為了解復雜疾病機制提供了許多新線索[1-2]。識別低頻或稀有遺傳變異的主效應及無主效應的交互作用(也稱上位效應)被認為是解釋復雜疾病“遺傳缺失”的重要原因[3-5]。隨機森林(random forest，RF)算法[6]通過構建多個決策樹對每個變異進行綜合評分，可以識別無強主效應的交互作用，但是其檢測交互作用的能力本質上仍取決于主效應[7]。BGTA(backward genotype-trait association)，也稱反向遺傳關聯算法[8]，通過考慮遺傳變異聯合效應來發現無主效應但存在交互作用的變異。近年來，隨著生物信息學技術迅猛發展，基于生物通路進行遺傳關聯分析在解決維度困擾、提高效能以及提高生物學解釋等方面具有顯著優勢[9]。本文利用GAW19(genetic analysis workshop 19)數據庫[10]，先基于通路分析策略篩選高血壓相關通路，再運用兩階段BGTA算法構建包括低頻遺傳變異的遺傳交互網絡，探索常見疾病低頻變異間交互作用，并與隨機森林聯合決策樹方法比較，為多遺傳變異交互作用研究提供方法學支持。

資料和方法

1.數據來源和預處理

(1)GAW19數據庫簡介

本研究利用GAW19中1851個無關聯糖尿病個體的外顯子測序數據和真實表型數據，其中外顯子測序數據包括奇數染色體上經質量控制后符合要求的約169萬個單堿基變異(single nucleotide variants，SNVs)，表型數據包括收縮壓、舒張壓、年齡和性別等。本研究將真實高血壓表型作為結局，定義為收縮壓>140mmHg或舒張壓>90mmHg，包括高血壓427例，非高血壓1424例。

(2)候選SNVs篩選

KEGG是整合基因、酶和化合物及代謝網絡信息的綜合性數據庫[11-12]。本文利用其通路子數據庫檢索出高血壓相關腎素-血管緊張素通路；基于疾病子數據庫檢索出原發性高血壓相關基因，并定位到醛固酮合成與分泌通路及醛固酮調節鈉吸收通路；最后將三條通路整合成RAAS(renin-angiotensin-aldosterone system)，也稱腎素-血管緊張素-醛固酮系統通路，共包含96個高血壓相關基因，其中53個位于奇數染色體，編碼26個蛋白或基因產物。

其后在ensembl(http：//grch37.ensembl.org/index.html)中提取奇數染色體上53個基因的位置信息，并與GAW19數據庫中的SNVs位置匹配，共獲得12251個SNVs。進一步根據條件(1)哈迪溫伯格平衡檢驗P≥0.05；(2)最小等位基因(minimum allele frequency，MAF)>1%；(3)在連鎖不平衡r2>0.8的連鎖域采用標簽策略選擇標簽SNVs；共篩選出248個SNVs納入分析，其中包含110個低頻變異(MAF在1～5%)。為便于表達，本研究對所有SNVs按1，2，…，248順序編號，重要SNVs的位置及基因信息見表1。

2.BGTA 算法原理

(1)GTD(genotype-trait distortion)和GTA(genotype-trait association)得分

假設有k個SNVs，每個SNV有三種基因型，則共有3k個多變異組合形式。設nD和nU分別表示病例組和對照組人數，nD，s和nU，s分別表示病例組和對照組中k個變異集的聯合基因型頻數，可計算病例組和對照組的聯合基因型頻率差異即GTD得分(式(1))，表示多個變異集與疾病相關的程度。

(1)

計算k個變異集{M1，M2，…，Mk}的GTDk得分與無Mr的k-1個變異集{M1，M2，…，Mr-1，Mr+1，…，Mk}的GTDk-1得分差異，可反映單變異Mr與疾病關聯的程度，即GTA得分(式(2))。

(2)

(2)基于隨機子集的BGTA篩選

為將交互作用SNVs同時選出，Lo等人[13-14]采用重復性地隨機選取變異子集的方式篩選，最大化地利用數據信息。具體流程如下：

①從K個SNVs的數據集中隨機選擇k個變異，一般k為2～10。

②計算每個變異Mr的GTA得分，若得分<0，則Mr保留；若得分≥0，則剔除GTA得分最大的SNV，保留k-1個SNVs。

③重復②，直至子集中GTA得分均為負，定義為關聯子集，并計算關聯子集的GTD得分。

④重復①到③B次。

⑤計算B個隨機子集中保留SNVs的出現頻次及關聯子集的GTD得分并排序。

可以看到，BGTA算法考慮隨機子集的聯合效應，無主效應時聯合效應考慮為交互效應。對較為重要的SNVs(有主效應或交互效應)，BGTA算法有較高的保留頻次。

(3)隨機子集變異數k和重復次數B的確定

考慮樣本量和計算量，k一般取2～10。k越大，同時分析的變異越多[8]。重復次數B主要取決于變異總數K及交互作用。假設k個變異集中的某變異Mr與疾病存在關聯，其保留頻率為p1；若與疾病不存在關聯，其保留頻率為p2。則重復抽取次數B滿足[14]：

(3)

式(3)中，無關聯變異保留頻率p2≈1/K，若Mr無主效應或弱主效應，但存在交互作用，則Mr保留頻率p1為：

(4)

實際應用中，重復值B比理論上小得多，因為Mr可能與多個SNVs存在交互作用或邊際效應，將增加p1的概率。

(4)二階段BGTA篩選，并構建遺傳交互網絡及驗證

考慮到運算速度及探索低階交互作用更具解釋性，本文采用二階段BGTA算法篩選SNVs。第一階段定義k=10，選擇GTD得分前100的關聯子集變異。第二階段在第一階段基礎上，定義k=2，篩選一階交互作用變異，再根據各變異對子集的GTD得分進行置換檢驗，并采用FDR(false discovery rate)進行校正，檢驗水準為0.10。用最終篩選出的SNV對構建遺傳變異交互網絡，并映射成基因交互網絡。本研究對篩選的一階交互作用采用logistic回歸驗證(檢驗水準為0.05)，并與隨機森林(random forest，RF)聯合決策樹方法比較。

3.隨機森林聯合決策樹方法

本研究同時采用RF聯合決策樹方法，根據RF重要性評分和袋外估計誤差進行第一階段篩選。重要性評分基于MDG(mean decrease Gini)和MDA(mean decrease accuracy)兩個指標[6]；接著基于CHISQ分類規則構建決策樹進行第二階段低階交互篩選，檢驗水準為0.10，并應用logistic回歸模型進行驗證(檢驗水準為0.05)。

4.統計分析軟件

在SAS 9.4軟件中編程實現數據清洗和統計檢驗，以及基因交互作用信息計算和置換檢驗。使用BGTA軟件包(http：//statgene.stat.columbia.edu/)完成BGTA算法。變異對映射遺傳交互網絡采用R軟件的ggplot2包和igraph包[15]。RF采用R軟件中的Randomforest包，決策樹采用SPSS 17.0實現。

結 果

1.兩階段BGTA算法篩選變異及構建遺傳交互網絡

(1)BGTA第一階段篩選

第一階段設置k=10、B=150萬，采用BGTA算法對248個候選SNVs進行篩選，結果顯示GTD得分前100位的關聯子集共包含76個遺傳變異(含61個低頻變異)，分屬于36個基因，編碼22個蛋白。計算各SNV的保留頻次，其中大于分位數閾值Q3+1.8×(Q3～Q1)的有6個(圖1)。χ2檢驗主效應結果顯示其中4個(26，49，101和184)SNVs的P<0.05，FDR校正后均無統計學關聯(表2)，推測這些SNVs的作用可能為與其他SNVs的交互作用。

圖1 BGTA算法第一階段SNVs保留頻次

SNVs基因MAFP?FDR26ATP1A40.230.0140.0748PIK3CD0.210.0350.1049PIK3R30.350.1040.0791PIK3CA0.470.0330.09101CAMK2A0.490.0750.07184ANPEP0.420.0470.07

(2)第二階段交互作用分析

第二階段設置k=2、B=50萬，對第一階段篩選的76個變異分析，結果保留了1102對不可約子集。隨后的1000次置換檢驗和FDR校正后共得到82對SNVs，包含56個遺傳變異(含44個低頻變異)。進一步構建遺傳交互網絡，可以看到，26(ATP1A4)、49(PIK3R3)、52(REN)、184(ANPEP)和247(THOP1)與其他SNVs存在較多的連線(圖2)，作為可能的樞紐變異，值得重點關注。進一步將SNVs交互網絡映射到基因交互網絡見圖3。

圖2 高血壓遺傳變異交互網絡

圖3 遺傳變異網絡映射到基因交互網絡

(3)logistic驗證

采用logistic回歸驗證82對SNVs的相乘和相加交互作用。如表3所示，共有16對(21個)SNVs存在相乘或相加交互，其中相乘交互12對，相加交互10對；21個SNVs中含低頻變異15個，其中26、49和169與其他SNVs的交互作用對分別為6、3和3對。

表3 BGTA的logistic交互作用驗證結果

2.隨機森林聯合決策樹分析

(1)第一階段RF篩選

利用RF基于MDG共篩選出35個SNVs(無低頻變異)。χ2檢驗顯示其中7個(20、83、49、101、184、186、75)與高血壓關聯有統計學意義(P<0.05)。基于MDA共篩選出69個SNVs(含30個低頻變異)。χ2檢驗顯示9個(83、92、75、192、101、154、247、148、183)與高血壓關聯有統計學意義(P<0.05)，但經FDR校正后均無統計學意義(見表4)。

表4 采用RF第一階段篩選卡方檢驗有統計學意義SNVs信息

(2)第二階段決策樹分析低階交互

采用CHISQ分類規則構建決策樹。基于MDG篩選結果，修剪后保留3層，篩選出5個SNVs(186、20、206、101和123)，其中206和123無主效應，186和206，20、206和101可能存在交互作用。基于MDA篩選結果，修剪后保留3層，篩選出7個SNVs(83、101、48、177、206、144和1)，除83和101外其他SNVs無主效應，83和48，48、144和1，101和177可能存在交互作用。

(3)logistic驗證結果

采用logistic回歸驗證決策樹結果。基于MDG，有2對變異186和206、20和123存在無主效應的相乘交互，同時也可能存在相加交互。基于MDA，有2對相乘交互和3對相加交互(見表5)。

表5 決策樹的logistic交互作用驗證結果

3.兩種方法結果比較

篩選變異的第一階段，BGTA比RF篩選出較多低頻變異和編碼蛋白信息，見表6。在低階交互作用分析的第二階段，BGTA仍篩選出較多低頻變異。盡管主要的交互作用都可在logistic回歸中得到驗證，但BGTA比決策樹識別的無主效應交互作用更多，見表7。

表6 BGTA與RF篩選交互作用變異數比較

表7 BGTA算法與決策樹識別的交互作用對結果比較

討 論

研究顯示，基因-基因間的交互作用(上位效應)是解釋復雜性狀“遺傳缺失”的重要原因之一，尋找廣泛的基因-基因交互作用有助于了解復雜疾病的生物學機制[16-17]。很多情況下遺傳變異對疾病沒有主效應，而是通過基因-基因間的交互效應對疾病造成影響，而傳統方法幾乎很難直接檢測到這種無主效應的交互效應[18-19]。BGTA算法計算病例對照中多變異聯合基因型的差異，并采用逐一向后剔除法來評價各變異的重要性。當無主效應或弱主效應的SNVs仍有較高保留頻次時，其重要性表現為該SNV與其他SNVs的交互作用；特別是如與多個其他SNVs交互，則代表疾病關聯通路上潛在的樞紐基因或關鍵變異。本研究將BGTA算法用于GAW19數據，分析高血壓相關通路基因交互作用。結果顯示和隨機森林相比，BGTA算法能有效地識別無主效應的交互作用，并對低頻變異有較高的效能，借此構建的遺傳交互網絡有更強的生物解釋性。

已有的研究顯示RAAS系統的激活與高血壓發生存在明顯的關聯[20]。作為一種常見復雜性狀，高血壓的遺傳力為20%～60%[21-22]。但是，即便在已確認的高血壓關聯生物通路中，由GWAS發現的遺傳變異對遺傳度的解釋也十分有限[23]。如近期的一項高血壓GWAS僅識別到RAAS通路中少量遺傳變異的主效應，包括前期已識別的基因REN、ACE、AGT、CYP11B2等[24]，因此，探索遺傳變異中更廣泛的交互作用是填補現有“遺傳缺失”的重要途徑。

本文通過兩階段BGTA算法篩選RAAS候選通路中可能存在交互作用的遺傳變異并構建遺傳交互網絡，顯示該通路中存在大量的無主效應、包括很多低頻變異間的交互作用，其中26(ATP1A4)、247(THOP1)、49(PIK3R3)、183(ANPEP)、52(REN)與其他SNVs存在多項交互作用，顯示在遺傳網絡中可能起著關鍵的樞紐作用。logistic回歸驗證也發現交互作用主要集中在少數變異如26、49、169與其他變異間。進一步將遺傳變異交互網絡映射到基因交互網絡，可直觀顯示基因REN、ACE、ATP1A4、PIK3R3、THOP1、ANPEP、DAGLB等間的高度連接，揭示了RAAS通路內基因的廣泛交互作用。特別是PIK3R3、ATP1A4等作為潛在的關鍵樞紐基因在后續生物學機制研究中更值得關注。這些結果與KEGG分子信號通路中的基因關聯也較為一致。

盡管如此，本研究也存在一些局限性。在研究設計上，由于GAW19的遺傳數據僅提供奇數染色體的基因，所以本研究遺傳交互網絡不能代表完整的RAAS生物通路。在算法原理上，BGTA算法基于聯合效應，可很好地識別無主效應的交互，但當變異有較強主效應時可能會高估變異間交互作用；此外，對于MAF更低的稀有變異，BGTA算法仍面臨挑戰。

總之，結合通路分析策略，BGTA算法可探索復雜疾病全基因組關聯研究中的重要遺傳變異，識別復雜疾病的遺傳交互作用，尤其是無主效應時低頻變異交互及關鍵樞紐基因，并通過構建遺傳交互網絡可視化遺傳結構，為理解復雜性狀的遺傳生物學機制提供有益參考。