基于通路分析的遺傳交互網(wǎng)絡(luò)方法：BGTA方法及應(yīng)用*

任志強(qiáng) 皮路程 劉桂炎 劉 晴 程金群 郜艷暉△

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

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

資料和方法

1.數(shù)據(jù)來(lái)源和預(yù)處理

(1)GAW19數(shù)據(jù)庫(kù)簡(jiǎn)介

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

(2)候選SNVs篩選

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

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

2.BGTA 算法原理

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

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

(1)

計(jì)算k個(gè)變異集{M1，M2，…，Mk}的GTDk得分與無(wú)Mr的k-1個(gè)變異集{M1，M2，…，Mr-1，Mr+1，…，Mk}的GTDk-1得分差異，可反映單變異Mr與疾病關(guān)聯(lián)的程度，即GTA得分(式(2))。

(2)

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

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

①?gòu)腒個(gè)SNVs的數(shù)據(jù)集中隨機(jī)選擇k個(gè)變異，一般k為2～10。

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

③重復(fù)②，直至子集中GTA得分均為負(fù)，定義為關(guān)聯(lián)子集，并計(jì)算關(guān)聯(lián)子集的GTD得分。

④重復(fù)①到③B次。

⑤計(jì)算B個(gè)隨機(jī)子集中保留SNVs的出現(xiàn)頻次及關(guān)聯(lián)子集的GTD得分并排序。

可以看到，BGTA算法考慮隨機(jī)子集的聯(lián)合效應(yīng)，無(wú)主效應(yīng)時(shí)聯(lián)合效應(yīng)考慮為交互效應(yīng)。對(duì)較為重要的SNVs(有主效應(yīng)或交互效應(yīng))，BGTA算法有較高的保留頻次。

(3)隨機(jī)子集變異數(shù)k和重復(fù)次數(shù)B的確定

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

(3)

式(3)中，無(wú)關(guān)聯(lián)變異保留頻率p2≈1/K，若Mr無(wú)主效應(yīng)或弱主效應(yīng)，但存在交互作用，則Mr保留頻率p1為：

(4)

實(shí)際應(yīng)用中，重復(fù)值B比理論上小得多，因?yàn)镸r可能與多個(gè)SNVs存在交互作用或邊際效應(yīng)，將增加p1的概率。

(4)二階段BGTA篩選，并構(gòu)建遺傳交互網(wǎng)絡(luò)及驗(yàn)證

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

3.隨機(jī)森林聯(lián)合決策樹方法

本研究同時(shí)采用RF聯(lián)合決策樹方法，根據(jù)RF重要性評(píng)分和袋外估計(jì)誤差進(jìn)行第一階段篩選。重要性評(píng)分基于MDG(mean decrease Gini)和MDA(mean decrease accuracy)兩個(gè)指標(biāo)[6]；接著基于CHISQ分類規(guī)則構(gòu)建決策樹進(jìn)行第二階段低階交互篩選，檢驗(yàn)水準(zhǔn)為0.10，并應(yīng)用logistic回歸模型進(jìn)行驗(yàn)證(檢驗(yàn)水準(zhǔn)為0.05)。

4.統(tǒng)計(jì)分析軟件

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

結(jié) 果

1.兩階段BGTA算法篩選變異及構(gòu)建遺傳交互網(wǎng)絡(luò)

(1)BGTA第一階段篩選

第一階段設(shè)置k=10、B=150萬(wàn)，采用BGTA算法對(duì)248個(gè)候選SNVs進(jìn)行篩選，結(jié)果顯示GTD得分前100位的關(guān)聯(lián)子集共包含76個(gè)遺傳變異(含61個(gè)低頻變異)，分屬于36個(gè)基因，編碼22個(gè)蛋白。計(jì)算各SNV的保留頻次，其中大于分位數(shù)閾值Q3+1.8×(Q3～Q1)的有6個(gè)(圖1)。χ2檢驗(yàn)主效應(yīng)結(jié)果顯示其中4個(gè)(26，49，101和184)SNVs的P<0.05，F(xiàn)DR校正后均無(wú)統(tǒng)計(jì)學(xué)關(guān)聯(lián)(表2)，推測(cè)這些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)第二階段交互作用分析

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

圖2 高血壓遺傳變異交互網(wǎng)絡(luò)

圖3 遺傳變異網(wǎng)絡(luò)映射到基因交互網(wǎng)絡(luò)

(3)logistic驗(yàn)證

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

表3 BGTA的logistic交互作用驗(yàn)證結(jié)果

2.隨機(jī)森林聯(lián)合決策樹分析

(1)第一階段RF篩選

利用RF基于MDG共篩選出35個(gè)SNVs(無(wú)低頻變異)。χ2檢驗(yàn)顯示其中7個(gè)(20、83、49、101、184、186、75)與高血壓關(guān)聯(lián)有統(tǒng)計(jì)學(xué)意義(P<0.05)?；贛DA共篩選出69個(gè)SNVs(含30個(gè)低頻變異)。χ2檢驗(yàn)顯示9個(gè)(83、92、75、192、101、154、247、148、183)與高血壓關(guān)聯(lián)有統(tǒng)計(jì)學(xué)意義(P<0.05)，但經(jīng)FDR校正后均無(wú)統(tǒng)計(jì)學(xué)意義(見表4)。

表4 采用RF第一階段篩選卡方檢驗(yàn)有統(tǒng)計(jì)學(xué)意義SNVs信息

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

采用CHISQ分類規(guī)則構(gòu)建決策樹?；贛DG篩選結(jié)果，修剪后保留3層，篩選出5個(gè)SNVs(186、20、206、101和123)，其中206和123無(wú)主效應(yīng)，186和206，20、206和101可能存在交互作用?；贛DA篩選結(jié)果，修剪后保留3層，篩選出7個(gè)SNVs(83、101、48、177、206、144和1)，除83和101外其他SNVs無(wú)主效應(yīng)，83和48，48、144和1，101和177可能存在交互作用。

(3)logistic驗(yàn)證結(jié)果

采用logistic回歸驗(yàn)證決策樹結(jié)果。基于MDG，有2對(duì)變異186和206、20和123存在無(wú)主效應(yīng)的相乘交互，同時(shí)也可能存在相加交互?；贛DA，有2對(duì)相乘交互和3對(duì)相加交互(見表5)。

表5 決策樹的logistic交互作用驗(yàn)證結(jié)果

3.兩種方法結(jié)果比較

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

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

表7 BGTA算法與決策樹識(shí)別的交互作用對(duì)結(jié)果比較

討 論

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

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

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

盡管如此，本研究也存在一些局限性。在研究設(shè)計(jì)上，由于GAW19的遺傳數(shù)據(jù)僅提供奇數(shù)染色體的基因，所以本研究遺傳交互網(wǎng)絡(luò)不能代表完整的RAAS生物通路。在算法原理上，BGTA算法基于聯(lián)合效應(yīng)，可很好地識(shí)別無(wú)主效應(yīng)的交互，但當(dāng)變異有較強(qiáng)主效應(yīng)時(shí)可能會(huì)高估變異間交互作用；此外，對(duì)于MAF更低的稀有變異，BGTA算法仍面臨挑戰(zhàn)。

總之，結(jié)合通路分析策略，BGTA算法可探索復(fù)雜疾病全基因組關(guān)聯(lián)研究中的重要遺傳變異，識(shí)別復(fù)雜疾病的遺傳交互作用，尤其是無(wú)主效應(yīng)時(shí)低頻變異交互及關(guān)鍵樞紐基因，并通過(guò)構(gòu)建遺傳交互網(wǎng)絡(luò)可視化遺傳結(jié)構(gòu)，為理解復(fù)雜性狀的遺傳生物學(xué)機(jī)制提供有益參考。