曹一鳴 胡 曼 徐永利*

(1.北京化工大學 數理學院， 北京 100029； 2.北京兒童醫(yī)院 眼科， 北京 100045)

引 言

基因的多位點組合與表型的關聯分析是生物信息學的研究熱點，近年來的研究表明，在對復雜性遺傳疾病的分析上，相比于單基因位點，多基因的多位點組合與表型的關聯更為顯著。然而，現有的方法僅適用于分析多位點組合與離散型表型的關聯，不能分析其與連續(xù)型表型的關聯。

原發(fā)性青光眼屬于復雜遺傳病[1]，是全球第二位的不可逆性致盲眼病[2-3]，在我國，原發(fā)性青光眼(不包括先天性青光眼)的致盲率約為10%[4]。原發(fā)性開角型青光眼(primary open-angle glaucoma，POAG)的單眼致盲率在北京市農村和市區(qū)分別為15.4%和10.90%[5]。在我國，患青光眼的人群中原發(fā)性以及先天性青光眼患者的比例占80%以上[6]。針對青光眼進行基因層面的進一步研究有助于對該疾病的早期診斷。病理性眼壓升高是大部分青光眼患者都會出現的臨床癥狀，在一項大型多隊列研究中，發(fā)現眼壓(intraocular pressure， IOP)與原發(fā)性開角型青光眼之間存在著強烈的關聯[7]。因此研究與眼壓相關的變異有助于發(fā)現更多影響青光眼發(fā)病的基因。

測序技術的進步使研究人員可對大規(guī)模樣本數據進行測序，同時可縱向比較一大批樣本的基因組。全基因組關聯分析(genome-wide association study，GWAS)是一種在全基因組層面上的無假設分析(即關于特定基因或基因座位無現有的假設，但也不存在這些基因之間關聯的零假設)，其對跨基因組和疾病的標記(單核苷酸多態(tài)性，即SNP)之間的關聯進行測試，通常涉及300 000個或更多合理多態(tài)且在基因組中平均分布的標記[8]。近年來，GWAS試圖進一步研究多個SNP對個體表型差異的綜合作用，從而衡量它們對疾病發(fā)病的影響大小。

傳統(tǒng)GWAS一般通過對Case-Control數據的分析來篩選出能夠影響復雜遺傳病的SNP位點。2012年，一項基于單位點對青光眼作用的GWAS研究[9]給出了結果，即經過T檢驗后具有較高顯著性的兩個單位點：位于17號染色體基因GAS7上的rs11656696和位于1號染色體基因TMCO1上的rs7555523。然而對于復雜遺傳病，它在基因層面的致病機理不能僅通過幾個SNP位點單獨的作用得到解釋，考慮多個位點的共同作用和上位性相互作用等一系列基因間的相互作用關系非常必要[10]。

近年來，一些學者通過對Case-Control數據的分析，在多個SNP位點組合與疾病關聯性的研究方面已經取得了一些成果。Sapin等[10]利用蟻群算法對Ⅰ、Ⅱ型糖尿病，腸炎和類風濕性關節(jié)炎4種疾病進行了研究，在4個數據集上均發(fā)現了對疾病影響統(tǒng)計顯著性較高的雙SNP位點組合。Liu等[11]利用Hi- seeker模型在對乳腺癌和小腸吸收不良癥的Case-Control數據研究中，發(fā)現了對疾病影響統(tǒng)計顯著性較高的三SNP位點組合。以上兩種方法都是針對Case-Control數據(每個病例的狀態(tài)是患病或健康)的分析，而與POAG之間存在強烈關聯的眼壓數據是連續(xù)變化型數據，每個病例的眼壓值可以是確定范圍內的任意實數，而不是有限的幾個狀態(tài)。據我們所知，針對連續(xù)型數據和SNP位點的關聯性分析目前還沒有相關研究。

本文設計了一種蟻群算法，用來從數以億計(約1.767×1011對)的雙位點組合里篩選與眼壓最為相關的組合。利用本文設計的蟻群算法，對歐洲生物信息學研究所(European Bioinformatics Institute，EBI)所轄數據庫中一項記錄了樣本人群眼壓和基因的原發(fā)性青光眼數據進行分析，并從中發(fā)現了一些具有較高統(tǒng)計顯著性的雙位點組合；特別是在2號染色體上發(fā)現了多個SNP位點與同一個基因具有較高的顯著性關聯，這意味著在此數據集上可能存在一些特定位點，通過上位性相互作用影響其他基因的表達量。

1 算法部分

1.1 蟻群算法

蟻群算法是一種啟發(fā)式算法，能夠在節(jié)省計算資源的前提下從大量數據中發(fā)掘出有效信息。相比于遍歷方法的計算復雜度O(n2) (n為總位點數)，蟻群算法的計算復雜度為O(m) ，其中

m=genmax×n_ant

(1)

式中，genmax表示最大迭代次數，n_ant表示螞蟻個數。相比于遍歷計算平方階的時間復雜度O(n2)，蟻群算法的計算復雜度為線性階O(m)。

本文借鑒了Sapin等[10]用來分析Case-Control數據的算法，設計出一種通過研究連續(xù)型變量和雙SNP位點組合之間的關聯，來尋找可能具有致病性雙位點組合的蟻群算法。通常的蟻群算法是一種探索優(yōu)化路徑的概率型算法，螞蟻個數表示每一次算法迭代中可行解個數的參數。螞蟻從隨機起點開始，移動到與當前節(jié)點有連接的下一個節(jié)點上。在選擇下個節(jié)點的過程中，螞蟻會根據滿足條件的節(jié)點上信息素值的大小，對較大信息素值有偏向性地隨機決定下一步，某節(jié)點信息素值越高則螞蟻選擇該節(jié)點的概率越高。節(jié)點被選擇后，螞蟻在選擇的節(jié)點上按一定規(guī)則留下一定量的信息素。在本文設計的算法中，信息素用來刻畫每個SNP位點對于預測眼壓的重要程度。

提前為算法設置初始信息素值、螞蟻個數、競賽規(guī)模、禁忌頻率這4個參數。另外，由于沒有明確終止條件，需要設定最大迭代次數genmax。設置參數后進入算法的循環(huán)部分。在算法的每一次循環(huán)中，首先采用競賽選擇為每一只螞蟻找到一個初始的雙位點組合，即讓每只螞蟻隨機且具有偏向性地從數據集子集中選擇SNP組合；然后評估每只螞蟻攜帶的雙位點組合對受測試者眼壓的劃分情況，并計算適應度函數；循環(huán)的最后，對本次循環(huán)最優(yōu)的雙位點組合進行信息素更新，相應信息素值P(SNP1)+P(SNP2)疊加在雙位點組合的每個SNP上；最后，在所有SNP位點上乘一個小于1的常數r_loss作為信息素值的損失率，以此來蒸發(fā)一定量的信息素值，本文中使用了0.99這個乘數值。

又考慮到某些位點的信息素值可能過高，從而頻繁出現在迭代中，因此本文算法結合了禁忌表每隔一定的子循環(huán)數(1 000次)禁忌信息素最高的一個位點。從生物學意義上來說，由于基因間相互作用的影響，兩個單體作用都不大的SNP位點也可能組成對青光眼發(fā)病影響很大的位點組合，這樣的組合較難被發(fā)掘，未被解釋的青光眼遺傳機制極有可能在這些位點組合中。結合禁忌表可以找到更多諸如此類對青光眼發(fā)病綜合作用大、而組合中單個位點影響都很小的雙位點組合。

傳統(tǒng)的蟻群算法使用輪盤賭法進行路徑選擇，輪盤賭法又稱為比例選擇法，其基本思想是各位點被選中的概率與其信息素大小成正比。該方法能選擇既隨機又偏向于具有最大信息素的路徑。由于本文基因數據包含數量巨大的SNP(單位點個數接近600 000個)，當變量數很高時，直接使用輪盤賭法對計算量消耗會成倍增長。因此，本文在輪盤賭法之前加上競賽選擇步驟。競賽選擇即先隨機選擇一個指定大小的位點子集，再在這個子集上用輪盤賭法選擇2個SNP。

1.2 算法流程

基于1.1節(jié)分析設計的算法，設計偽代碼如下。

1)初始化參數(螞蟻個數n_ant和最大迭代次數genmax)，將每個SNP位點上的信息素初始化為0；

2)repeat以下步驟：

3)for 所有螞蟻 do：

4)通過競賽選擇為每一只螞蟻找到兩個SNP位點的雙位點組合；

5)計算每一只螞蟻的雙位點組合的適合度函數值；

6)end for；

7)更新具有最高適應度函數值的兩個SNP位點的信息素；

8)for 所有SNP位點 do：

9)乘以損失率r_loss來蒸發(fā)一定量的信息素；

10)end for；

11)每隔一定迭代次數禁忌信息素最高的一個位點；

12)until最大迭代次數后終止循環(huán)。

上面描述的算法偽代碼中，主要包括3個步驟：

步驟1初始化參數(設置螞蟻個數，最大迭代次數，競賽規(guī)模和禁忌頻率)，將每個SNP上的信息素初始化為0，然后對步驟2迭代genmax次；

步驟2通過競賽選擇為每一只螞蟻找到包含兩個SNP位點的組合，計算每一只螞蟻位點組合的適合度函數，更新具有最高適應度函數值的兩個SNP的信息素，最后對所有SNP位點的信息素值乘以損失率r_loss來蒸發(fā)一定量的信息素；

步驟3每隔一定迭代次數(本文為每1 000次)禁忌信息素最高的一個SNP位點，當迭代次數達到最大迭代次數genmax后，終止算法循環(huán)。

1.3 基因型數據的數字表示，雙位點組合的邏輯關系及適應度函數的計算

在1.2節(jié)所述的算法流程中，適應度函數用來刻畫基于某一雙位點組合狀態(tài)將全體病例分組的組間眼壓差異。雙位點組合把所有對象劃分為陽性組(即該組眼壓各項統(tǒng)計參數偏高)和陰性組，這兩組的眼壓平均值之差越大表明這個雙位點組合對眼壓預測的意義越顯著。

對于確定的每個雙位點組合都有9種基因型及4種布爾邏輯運算關系。每一個SNP位點的堿基構成可以分為3類，即常見純合子(AA)、罕見純合子(aa)以及雜合子(Aa)，故兩個位點的組合有9種可能的基因型?？紤]了兩個位點之間的4種布爾邏輯運算[10]，即：

AND 當且僅當第一個SNP取特定值并且第二個SNP取特定值時，個體是陽性的；

OR 當且僅當第一個SNP或第二個SNP取特定值時，個體是陽性的；

AND NOT 當且僅當第一個SNP取特定值，且第二個SNP不取特定值時，個體是陽性的；

XOR 當且僅當兩個SNP中的一個具有特定值而另一個SNP不采用特定值時，個體是陽性的。

當給定了基因型組合和邏輯關系后，可以把樣本眼壓劃分為兩組，從而計算適應度函數值。9種基因型組成和4種邏輯關系形成36種組合方式，故可得到36個適應度函數值，取這36個值中的最大值作為該雙位點組合最終的適應度函數值。本文算法使用的適應度函數為

(2)

式中，Mp表示陽性一組的人數；Pp表示陽性一組的眼壓平均值；Mn表示陰性一組的人數；Pn表示陰性一組的眼壓平均值；Pm表示所有受測試者的眼壓平均值。

2 實驗部分

2.1 實驗數據及處理

利用本文的蟻群算法分析了來自EBI的澳洲人群青光眼患病情況調查數據。該數據包含2 627個合格受測試者樣本(未測到眼壓值的不合格樣本有134個)，在594 398個SNP位點上的基因數據和眼壓值?；驍祿雌渌诓煌旧w位置劃分為26個子數據(1～22、XY、X、Y、MT，XY表示X染色體和Y染色體的重疊部分，MT表示線粒體)，數據表每行記錄一個SNP位點的信息，前5列記錄該位點編號、所在位置和堿基成分，之后每3列對應一個受測試者是顯性純合子、雜合子或隱性純合子的概率。最終的輸入數據是一個594 398×2 627的基因數據矩陣和一個2 627×1的眼壓數據矩陣。

2.2 蟻群算法的參數設置

首先用包含8 713個位點的21號染色體進行預實驗，來選取算法的核心參數(螞蟻數)。計算21號染色體上任意雙位點組合的適應度函數值，并按該值降序排列這些組合，前750個組合具有統(tǒng)計顯著性。接下來，在其他參數相同的條件下只改變螞蟻數，把計算結果中屬于這750個位點對的個數作為評價算法計算效率的標準。在預實驗中，將螞蟻數分別設置為10、20和30，競賽規(guī)模都設置為20，迭代次數分別設置為300 000、150 000、100 000次(確保算法總的計算次數相同)。圖1展示了在不同螞蟻數下，隨著計算次數增加蟻群算法在預實驗子集上能搜索到的前750位組合個數。

從圖1中可以看到，螞蟻數目設置為20時，蟻群算法搜索到的關鍵位點對個數最多。所以在接下來的全位點實驗中將螞蟻數恒定地取為20，同時將競賽規(guī)模設置為20。另外，每1 000次迭代禁忌一次位點。

2.3 全部位點蟻群實驗結果

在全位點實驗中使用的計算機包含參數為Intel i7- i7- 6850K的CPU，2 400 MHz的32G內存以及Ubantu 14.04的操作系統(tǒng)。根據計算機性能將算法的最大迭代次數設置為800 000次。完成本文設計的蟻群算法的800 000次迭代總共用時220 h。表1列出了秩和檢驗P值最小的10個雙位點組合。在所有結果中，有6對雙位點組合的秩和檢驗P值小于1×10-4，秩和檢驗P值小于1×10-2的一共有137對結果。

表1 秩和檢驗P值最優(yōu)的前十組雙位點組合Table 1 Top ten optimal combinations of SNPs of rank sum test P-value

括號內數字為滿足條件的人數占總人數的比例。

在之前關于青光眼數據的眼壓與單位點研究中，得到與該疾病顯著相關的一個單位點是位于17號染色體基因GAS7上的rs11656696[9]，在本文所分析的EBI數據集上該位點的秩和檢驗P值為5.03×10-2。本文所發(fā)現的雙位點組合對于眼壓預測的秩和檢驗P值低于該單位點的秩和檢驗P值，也即我們發(fā)現的雙位點組合對眼壓的影響更具有統(tǒng)計顯著性。

本文實驗中，秩和檢驗P值最小的組合是位點rs661221與rs2654703，該組合的秩和檢驗P值為1.84×10-5，且這個位點對劃分得到的陽性組和陰性組的眼壓樣本均滿足正態(tài)分布。

進一步計算了位點rs661221與位點rs2654703的T檢驗P值，為6.9×10-12，顯著優(yōu)于文獻[9]。位點rs661221位于1號染色體的基因RIMS3上，位點rs2654703位于3號染色體的基因FGF12上。另外，從表1可以看到，P值最小的10組雙位點組合劃分出的陽性組中，眼壓高于2.793 kPa的病例比例明顯高于普通人群中的比例，而眼壓高于2.793 kPa在臨床上是青光眼發(fā)病的重要指征。這表明本文算法發(fā)現的雙位點組合與青光眼發(fā)病因素顯著相關。

圖2(a)展示了秩和檢驗P值小于1×10-2的結果在染色體上分布的熱圖，圖中色塊的灰度深淺代表位點對的多少。橫縱坐標為位點所在染色體位置編號(1～21)，組成位點分別位于橫縱坐標標示的染色體上。圖2(b)是將圖2(a)的內容反映在點和邊上的分布圖，圖中節(jié)點度數代表兩位點在同一染色體的結果個數，有2對雙位點組合的組成位點都在2號染色體上；3號和11號染色體各有1對雙位點組合的組成位點都在同一染色體上；其余染色體不包含雙位點組合的組成位點都位于同一染色體的情況。圖2(b)中連接各節(jié)點的邊的粗細代表組成雙位點組合的各位點分別位于兩不同染色體的組合個數。如圖2(b)所示，最粗的一條邊(染色體之間雙位點組合最多)在1號染色體和3號染色體之間，包含8對雙位點組合；其次為2號染色體和6號染色體之間的邊，包含6對雙位點組合。

圖2表明，137對顯著雙位點組合涉及2號染色體的個數最多，因此對涉及到2號染色體位點的所有結果進行了分析，發(fā)現2號染色體基因LOC107985962上的位點與其他11個染色體上的位點組成了顯著性較高的組合。其中包括位點rs7742590，該位點與基因FOXC1都位于6號染色體，且距離基因FOXC1大約1 600萬個堿基。據文獻[12]報道，6號染色體上基因FOXC1被懷疑與POAG發(fā)病高度相關。

如表2所示，11個雙位點組合都包含2號染色體基因LOC107985962上的位點。以上位點中有5個位于明確的基因上，在這5個基因中，有兩個基因的功能已被確認(分別為ERICH6- AS1和PDE8B)，其余3個基因的功能尚未被確認。ERICH6- AS1屬于RNA反轉錄基因，PDE8B屬于磷酸二酯酶調控基因。根據美國國家生物技術信息中心(National Center for Biotechnology Information，NCBI)數據庫上收錄的常規(guī)組織表達量數據，以上兩個基因都在胎盤組織具有一定表達量。因此可以猜測該基因通過基因間相互作用影響其他基因的表達，從而引起胎兒的先天性青光眼。本文所提蟻群算法選取的雙位點組合與生物醫(yī)學已經證實的結果具有較高的相合性。本文結果的生物學意義有待通過進一步的生物實驗來驗證。

表2 與位點rs13414786相關的雙位點組合Table 2 Site pair combinations involving SNP rs13414786

3 結論

本文的研究提供了一種通過探索連續(xù)型變量與多位點組合的關聯性來發(fā)現與疾病相關的SNP組合的新思路。通過10 000位點子集的實驗找到了合理的參數，提高了所設計算法的計算效率，同時也增強了算法對邊際效應較高雙位點組合的探索能力。另外，通過全部位點的實驗，在2號染色體上找到了一系列未被報道的顯著性水平較高的位點組合，這些位點組合很有可能通過醫(yī)學實驗得到進一步證實，從而為探索與眼壓相關的基因功能和青光眼的發(fā)病原理提供新的線索。