王淑棟,張善強,賀思程
(中國石油大學(華東)計算機與通信工程學院,青島 266580)
隨著人類基因組測序工作的逐步完成,大量的數(shù)據(jù)為全基因組關聯(lián)分析提供了豐富的素材,也涌現(xiàn)出許多數(shù)據(jù)分析方法[1-4].人類基因組計劃得出人類所有的基因共由39 000多個已經(jīng)編碼蛋白的基因序列以及30億堿基組成.而國際單體型圖計劃[5]得到了SNP的300萬個位點.兩個計劃的實施給生物學領域帶來了眾多的數(shù)據(jù)信息,為全基因組研究中提供了方便.GWAS因其優(yōu)勢得到了很多的應用.大量研究成果顯示關聯(lián)研究具有很多的優(yōu)勢[6].
Ghoussaini等[7]在2012年針對乳腺癌相關基因進行研究,共得到了3個致病相關的位點,rs10771399不僅在乳腺癌的發(fā)展中起著關鍵作用,在骨轉(zhuǎn)移中也有著同樣的重要性.2013年,維爾漢姆等[8]關于躁郁癥的數(shù)據(jù)進行分析,得出與躁郁癥相關的SNP位點及致病基因.2014年,廣川等[9,10]針對心肌梗塞病設計了病例對照實驗,從實驗中得到了有關疾病的致病基因和SNP,使心肌梗塞病得到了合理的解釋.
GWAS能夠幫助人們更好的解釋復雜疾病成因,但是它也有不足.一方面,復雜疾病多種多樣,其中的影響因素也很多,如何確切地得到與特定的功能相聯(lián)系的位點是個不小的難題;另一方面,對于GWAS結果,它在不同群體中的影響程度并不一樣;目前的大部分研究主要針對簡單疾病,沒有涉及到基因間的相互作用.
而針對基因間的相互作用,可以通過互信息建立網(wǎng)絡進行表達.GWAS網(wǎng)絡方法將GWAS數(shù)據(jù)進行網(wǎng)絡建模,通過比較疾病數(shù)據(jù)與對照數(shù)據(jù)得出的網(wǎng)絡的不同,進行后續(xù)的相關統(tǒng)計量的分析及解釋.
本文試圖通過互信息表示SNP之間的相互作用關系,進而建立SNP與SNP之間的網(wǎng)絡.在此基礎上,進行全基因組關聯(lián)研究,找到結構性關鍵SNPs.
隨著生物網(wǎng)絡的研究深入發(fā)展,研究者對元素之間的相關性的描述越來越準確,互信息作為兩個元素之間的相關信息度量,具有很多的優(yōu)勢,其中最具優(yōu)勢的就是它的熵表示,不僅是對元素出現(xiàn)概率的表示,更是體現(xiàn)了元素之間的離散程度及相互之間的關系,對于給定的兩個SNP表達序列,他們之間的數(shù)據(jù)存在著差異,而利用互信息可以充分表達SNP之間的差異性及依賴性,互信息越大,說明兩個SNP之間的關聯(lián)程度越緊密;反之,則說明聯(lián)系越小,從而找到跟所有的SNP聯(lián)系較大的節(jié)點,即是關鍵SNP.本文通過互信息建立相互作用網(wǎng)絡,從而分析網(wǎng)絡結構的差異性.設是兩個SNP的基因型數(shù)據(jù)在個體之間表達形成的向量,是X和Y的聯(lián)合概率分布,是他們之間的聯(lián)合熵,定義為:
對于兩個隨機變量之間存在的關系,H(X)表示隨機變量X蘊含的不確定性,而條件熵則是已知條件Y時隨機變量X所余下的不確定性,那樣,就表示已知條件Y后X包含的信息量.進而還可以證明這個值關于X和Y是對稱的,即且都等于由此X和Y之間的互信息可以計算,互信息記為
因為SNP數(shù)據(jù)是每個SNP仿真1000組得到的數(shù)據(jù),每三個數(shù)據(jù)代表一個個體,首先需要對數(shù)據(jù)進行處理使得數(shù)據(jù)能夠表示基因型,我們確定使用0,1,2三個數(shù)來表示每個個體內(nèi)表達的基因型,再根據(jù)公式(2)計算得到所有的SNP之間的互信息.具體計算過程如下:
(1)我們首先得到每個SNP的基因型可能性序列數(shù)據(jù),假設共有N個個體,則每一行包含2N個SNP堿基可能性數(shù)據(jù),0代表出現(xiàn),1代表不出現(xiàn).
例如:假定兩個個體關于5個SNPs的基因型數(shù)據(jù)如下:
SNP 1:AA AA
SNP 2:GG GT
SNP 3:CC CT
SNP 4:CT CT
SNP 5:AG GG
輸出的正確仿真數(shù)據(jù)如下所示:
SNP1 rs1 1000 A C 1 0 0 1 0 0
SNP2 rs2 2000 G T 1 0 0 0 1 0
SNP3 rs3 3000 C T 1 0 0 0 1 0
SNP4 rs4 4000 C T 0 1 0 0 1 0
SNP5 rs5 5000 A G 0 1 0 0 0 1
所以,在SNP3上,兩個等位基因上堿基分別為C和T,所以每個個體與之相對應的堿基組合CC,CT,TT出現(xiàn)的可能性序列分別是100和010.
(2)每個SNP的基因型表達數(shù)據(jù)作為一個向量,x,y表示來自SNP集合I中的其中的兩個SNP向量.
(3)根據(jù)每個SNP的基因型表達量的分布,計算得到每兩個SNP之間存在的互信息值.所有SNP之間的互信息構成互信息矩陣,記作矩陣中的每行代表一個SNP,每一列代表此SNP與另一個SNP之間的互信息.
假定存在一個集合的SNP基因型數(shù)據(jù)D,其中所擁有的SNP的集合我們記作I,可由互信息計算公式(2)得到一個互信息矩陣.定義一個建立在關于SNP基因型數(shù)據(jù)D的互信息網(wǎng)絡.是邊賦權圖,其中V表示點集合、每個網(wǎng)絡中的節(jié)點i∈V表示一個SNP,基因i和j之間的互信息計算值wij定義為每條邊的權重.在下面的表述中,我們將基因i∈I以及頂點i∈V等同起來看待.
利用上述方法得到的SNP相關網(wǎng)絡中各節(jié)點(SNP)的網(wǎng)絡結構參數(shù)來描述特定生物過程中基因的重要性.首先給出幾個重要的能夠反映網(wǎng)絡結構特點的網(wǎng)絡統(tǒng)計量的相關定義[11].
(1)度(K):在網(wǎng)絡中,度指的是與該點相連接的邊數(shù)目.節(jié)點度可以表示該點的重要程度,節(jié)點度越大,表示該點在網(wǎng)絡中越重要.而網(wǎng)絡的平均度可以通過計算所有的點的度,后取平均數(shù)計算得到.
(2)平均路徑長度(L):定義為網(wǎng)絡中所有的點之間兩兩求得的距離的平均數(shù),網(wǎng)絡中的任意兩點i,j的距離即邊的條數(shù),則兩點之間的平均路徑長度表示為所有的點之間的平均距離,記作:其中N表示網(wǎng)絡中的節(jié)點數(shù)目.
(3)聚類系數(shù)(C):網(wǎng)絡中節(jié)點i有Ki個邊與之連接,那么與該點可能連接的最大邊數(shù)為若這Ki個節(jié)點之間真實邊為Ei,則它與總的所有情況下的邊比例,計算得到節(jié)點i的聚類系數(shù)很顯然,0≤C≤1.C=0代表網(wǎng)絡中的點為孤立點;C=1表示網(wǎng)絡中的所有點之間都是互相連接的,視為全局耦合網(wǎng)絡.
(4)介數(shù)(B):網(wǎng)絡中介數(shù)的概念可以分為兩類,一類是點介數(shù),另一類是邊介數(shù).節(jié)點k的介數(shù)定義為,其中,C(i,j)代表i與j間最短路徑總數(shù),表示中間點為k時,i與j間的所有路徑總數(shù).介數(shù)反映了節(jié)點k在i和j之間的流通量和重要程度.網(wǎng)絡中某個節(jié)點的介數(shù)越大,說明該點在網(wǎng)絡中信息傳播的信息量就越大,越容易在該點造成網(wǎng)絡堵塞.假設兩組連接度很高的網(wǎng)絡中間只有少數(shù)點連接,那么這幾個少數(shù)點介數(shù)就會很大,即很多的信息在流通的過程中經(jīng)過這幾個點,很容易造成堵塞,從而造成數(shù)據(jù)信息丟失.因此,最大介數(shù)的增大會降低網(wǎng)絡同步能力.
(5)模塊度(Q):模塊度也稱作模塊化度量值,是用來衡量網(wǎng)絡強度的統(tǒng)計量.最早是Newman提出的,它用來描述網(wǎng)絡社團以及劃分的好壞.假定網(wǎng)絡共分為k個社團,代表一個k×k維的矩陣.故模塊度可以定義為:其中,是矩陣中的數(shù)值之和(行或列),eij用來表示社區(qū)i和社區(qū)j之間的邊的數(shù)量.模塊度可以區(qū)分社區(qū)劃分的好壞.若是劃分的好,則社區(qū)內(nèi)部節(jié)點相似度較大,而在社區(qū)外邊相似度較低.Q越大,越接近1,代表社區(qū)擁有一個很好的劃分結構,使得社區(qū)的劃分合理化.通常設定的值是在0.3與0.7之間.
本文中我們主要選擇5個參數(shù)進行分析比較,對于給定的參數(shù)進行最終的分析,從而找到影響網(wǎng)絡的重要因素,依據(jù)此類統(tǒng)計量進行歸納分析,得出相應的參數(shù).
我們對由SNP數(shù)據(jù)設定不同的互信息閾值而形成網(wǎng)絡,針對其中大于閾值的邊,做去掉處理,而針對小于閾值的邊進行保留操作,從網(wǎng)絡圖可以分析出統(tǒng)計量所對應的參數(shù)變化,得到有益信息量.
根據(jù)網(wǎng)絡中SNP之間互信息計算的值,選擇閾值范圍為0.1到0.63.共設置63個閾值,在每個閾值的條件下,統(tǒng)計計算相應的網(wǎng)絡結果,從而得到一致性網(wǎng)絡,根據(jù)網(wǎng)絡的相似性程度選擇對實驗組和對照組差別較大的統(tǒng)計量進行分析.我們最終選擇了度作為區(qū)分依據(jù),并分析能夠區(qū)分實驗組和對照組的取值范圍,得出最佳的閾值,對于不同的數(shù)據(jù),得到的互信息值也不同,所以需要根據(jù)數(shù)據(jù)得到的互信息范圍,由網(wǎng)絡統(tǒng)計量得到取值范圍,得到互信息取值的交集,能夠區(qū)分對照組和實驗組數(shù)據(jù),從而確定最佳的互信息閾值.這樣就能夠保證所取的閾值不受樣本數(shù)量的大小影響,而是根據(jù)樣本的不同情況得到相應的閾值.對于節(jié)點i,我們定義,?d代表了這個節(jié)點的度差異值,在該公式中,分別代表了這個節(jié)點在實驗組與對照組網(wǎng)絡中節(jié)點的度.
我們都知道,在復雜網(wǎng)絡中,節(jié)點度能夠代表節(jié)點的作用和影響力.本文從網(wǎng)絡結構差異的角度去衡量各個統(tǒng)計量[12],進而對應到其中的節(jié)點,找到“結構性關鍵SNPs”.這種差異性貢獻分為正、負貢獻兩個方面.我們用r代表度的變化閾值.正貢獻SNP代表了該節(jié)點在病例組、對照組兩個網(wǎng)絡中度的貢獻?d≥r的SNP;同理,負貢獻SNP代表了該節(jié)點在以上兩個網(wǎng)絡中度的貢獻?d≤-r的SNP.
本文對基因BRCA2仿真數(shù)據(jù)建立病例組與對照組建立相互作用網(wǎng)絡進行數(shù)據(jù)實驗.對SNP互信息設置一個閾值范圍,分析產(chǎn)生的病例組和對照組SNPs互信息網(wǎng)絡的統(tǒng)計量:平均路徑長度、聚類系數(shù)、平均度、模塊度、平均介數(shù)隨閾值在其變化范圍內(nèi)的增加而變化的情況.根據(jù)計算的網(wǎng)絡中SNP之間互信息的值,我們?nèi)』バ畔㈤撝档姆秶鸀?至0.63,步長0.01,分析對應病例組與對照組的SNP相互作用網(wǎng)絡的上述網(wǎng)絡結構參數(shù)隨變化而變化的情況.
HapMap給出了人類基因組單核苷酸多態(tài)性(SNPs)和拷貝數(shù)多態(tài)性(CNPs)的分布情況.本文使用HapMap提供的三個文件進行實驗,包含了關于BRCA2的88個SNPs.下面是對三個文件的說明.
.hap文件是已知的單體型數(shù)據(jù),其中行代表SNP,列表示單體型.每一個.hap文件都需要一個相應的legend文件,所有的等位基因都以0,1作為標記.
.legend文件是SNP標記位點數(shù)據(jù),四列數(shù)據(jù)分別表示SNP的ID、堿基位置、堿基的0,1表示.
.map文件包含了小規(guī)模的重組率,共三列分別表示每個SNP的物理位置,距離左標記點的位置和距離右標記點的位置.
在這數(shù)據(jù)中,必須去掉全部為0或者全部為1的數(shù)據(jù),因為這些數(shù)據(jù)對構建網(wǎng)絡結構沒有任何幫助.去掉這些多余的數(shù)據(jù),共得到45條SNP數(shù)據(jù).把3個文件放到一起,執(zhí)行Hapgen2軟件,代碼如下:
./hapgen2 -m BRCA2.map -l BRCA2.legend -h BRCA2.hap -o BRCA2.out -dl 31820136 1 2.5 2 31847382 0 1.5 4.5 -n 5000 5000.
分別仿真了5000組實驗組和對照組數(shù)據(jù).隨機選定2個SNPs作為致病SNPs.它們的信息如下:rs206081和rs9534318,選取雜合子變異率分別是2.5和1.5,純合子變異率分別為2和4.5,上述樣本數(shù)據(jù)都包含SNP編號,SNP位置及0,1表達數(shù)據(jù).
本文中,我們使用.gen文件,刪除前五列后把數(shù)據(jù)轉(zhuǎn)換成一個矩陣,其中每行表示一個向量,每三個數(shù)字代表一個個體,我們轉(zhuǎn)換成0,1,2表示.
根據(jù)得到的互信息矩陣,大于閾值的向量之間表示相互關系較強,選定這些SNP作為節(jié)點建立網(wǎng)絡.分析比較網(wǎng)絡的6個特性.每個結構參數(shù)都反映著網(wǎng)絡的特性,進而可以顯示SNP間的互信息的變化,取0.01為步長,從0到0.63之間求得每一個閾值下的網(wǎng)絡結構特性值,得到圖1.圖1中,縱坐標表示相應的統(tǒng)計量,橫坐標代表閾值,虛線表示對照組數(shù)據(jù)顯示效果,實線表示實驗組數(shù)據(jù)顯示效果.
圖1 4個網(wǎng)絡結構的統(tǒng)計量隨閾值的增加的變化情況
實驗發(fā)現(xiàn)5個結構特性中,平均聚類系數(shù)B交織在一起,不能區(qū)分實驗組和對照組.
觀察圖1(a),當0<t<0.21時,網(wǎng)絡的平均介數(shù)B在在兩組中的變化趨勢走向大體相似.當0.21<t<0.63時,網(wǎng)絡的平均介數(shù)B逐漸減小.從圖中可以明顯的看出,病例組的平均介數(shù)要比對照組的平均介數(shù)高.于是,我們得到,隨著互信息閾值的增大,節(jié)點的介數(shù)也在不斷減小,網(wǎng)絡中邊越來越稀疏.
觀察圖1(b),當0.2<t<0.43時,實驗組與對照組的網(wǎng)絡有相對明顯的差異.于是我們可以得到,在這個變化區(qū)間內(nèi),平均路徑長度可以很好的區(qū)分病例組和對照組,而當t>0.43時,網(wǎng)絡的邊越來越少,平均路徑長度趨近于0.
從模塊度Q隨閾值的變化圖1(c)看出,當閾值0<t<0.2 或 0.43<t<0.63 時,兩組中的模塊度Q逐步上升,但變化大致相同,而當 0.2<t<0.43 時,實驗組模塊度與對照組有較大區(qū)別.
觀察圖1(d),可以發(fā)現(xiàn),在很長的一段閾值范圍內(nèi),病例組與對照組的網(wǎng)絡平均度有很大的區(qū)別,而隨著網(wǎng)絡的閾值增加,網(wǎng)絡的平均度越來越小,這與網(wǎng)絡的孤立點越來越多也是相對應的.
當t>0.62時,病例、對照組中都只有一個包含四個節(jié)點的全耦合子網(wǎng),聚類系數(shù)C、平均路徑長度L兩者相等,且都為1.當t>0.63時,平均路徑長度L、聚類系數(shù)C是缺失的,平均介數(shù)B以及其他三個統(tǒng)計量值均為0.
總之,平均聚類系數(shù)C不能區(qū)分兩組數(shù)據(jù),平均路徑長度L和平均介數(shù)B能夠區(qū)分但是閾值具有一定局限性.平均度可以在很大的范圍內(nèi)把實驗組和對照組分別出來,我們選擇平均度作為區(qū)分的依據(jù).
從圖1中我們得到每個統(tǒng)計量能夠區(qū)分兩組的閾值范圍,如表1.
表1 各統(tǒng)計量能夠區(qū)分實驗組和對照組的閾值范圍
從表1可以看出,每一個統(tǒng)計量都有不同的閾值范圍,平均度K的范圍較大,0.08<K<0.35;其他的統(tǒng)計量閾值范圍相差不大,基本在0.2到0.3之間.結合圖1,選擇0.28為閾值構建網(wǎng)絡.
依據(jù)圖2,實驗組和對照組的圖像是有很大差異的.在對照組,節(jié)點之間聯(lián)系較弱且存在更多的孤立點.但是在實驗組中,很多的孤立點不再是獨立的,并且擁有了更多的聯(lián)系.對照組中存在36個連接點和9個孤立點,而實驗組中存在39個連接點喝6個孤立點.這表明我們選取的閾值0.28是合適的.經(jīng)過多次仿真數(shù)據(jù)試驗,對于結合數(shù)據(jù)互信息得到閾值范圍,而后確定互信息閾值的方法都是有效的.
圖2 閾值為0.28的條件下,實驗組和對照組互信息網(wǎng)絡
結構決定功能,而結構的差異決定了功能的差異,本文將這種差異細化到每個節(jié)點上,而平均度可以很好的區(qū)分病例組和病例組,所以我們選擇每個SNP位點的平均度來刻畫SNP在病例組和對照組的差異,計算每個網(wǎng)絡的每個節(jié)點的節(jié)點度差異,當節(jié)點的度在病例、對照組中的變化差異比較大時,說明這兩個組的網(wǎng)絡結構差異較大.從兩組網(wǎng)絡的數(shù)據(jù)分析來說,節(jié)點度的增量有正有負,所以,節(jié)點在病例組中的度也有增減之分,即存在正、負貢獻SNPs.度變化量增加最大的是節(jié)點39,增加值的大小是5,同理,減少量最大的是16,41,減少值的大小是2.
當閾值為0.28時,對照組網(wǎng)絡中的平均度大致等于2,從而可以得到,當病例、對照組網(wǎng)絡中節(jié)點度的變化值大于等于3時,其對網(wǎng)絡結構影響較大.故可設?d=3,由此,我們可以獲得對網(wǎng)絡結構有顯著影響4個SNPs,如表2,其中rs206081,rs9534318為預設致病SNPs.
表2 給定參數(shù)為3的條件下,部分結構性關鍵SNPs的信息及度的變化量
在查找“結構性關鍵SNPs”時,我們需要從網(wǎng)絡平均度出發(fā),對選取網(wǎng)絡中的關鍵SNPs設置合適的差值參數(shù).如果選取的差值參數(shù)比較小,對SNPs選取限制比較寬泛,一些不相關的SNPs也會選取到SNPs集合內(nèi),從而導致假陽性.反之,如果選取過于嚴苛,反而會遺漏一些比較重要的節(jié)點,導致假陰性.
我們選取基因BRCA2,得到它在閾值為0.28時候的網(wǎng)絡,如圖2所示.選擇不同的差值參數(shù),得到一系列不同的結構性關鍵SNPs,如表3所示.
表3 不同參數(shù)r的取值下關鍵SNPs個數(shù)
當互信息閾值設定為0.28時,網(wǎng)絡中度的最大變化量是5.當r≥5時,所得的關鍵SNPs只有節(jié)點39,對網(wǎng)絡影響較大的節(jié)點25卻被忽略.當r≤2時,所得的關鍵SNPs只有13個,這里面也包括了其中的非零點.
本文通過國際項目HapMap3中以及Hapgen2軟件生成的13號染色體上BRCA2基因生成仿真數(shù)據(jù).利用互信息表示SNPs間的相互作用.構建實驗組和對照組的網(wǎng)絡,根據(jù)閾值及差值參數(shù)篩選出關鍵SNPs.最后,對我們所選擇的參數(shù)進行了評估,證明我們所選定的參數(shù)能夠反映結構的變化,能夠較好地選擇出預設的關鍵SNPs.通過數(shù)值實驗發(fā)現(xiàn):樣本數(shù)目會影響互信息的大小,樣本數(shù)較小時,互信息較高,樣本數(shù)較大時,互信息逐漸降低,本文認為,樣本數(shù)偏少,則特異性個體數(shù)目不完備,樣本數(shù)過多,又會造成冗余,增加了計算復雜度.目前,確定合適的上下界仍然是一個具有挑戰(zhàn)的問題.
1 Pharoah PDP,Tsai YY,Ramus SJ,et al.GWAS metaanalysis and replication identifies three new susceptibility loci for ovarian cancer.Nature Genetics,2013,45(4):362–370e2.[doi:10.1038/ng.2564]
2 Xu ZL,Taylor JA.SNPinfo:Integrating GWAS and candidate gene information into functional SNP selection for genetic association studies.Nucleic Acids Research,2009,37(S2):W600–W605.
3 Larsson M,Duffy DL,Zhu G,et al.GWAS findings for human iris patterns:Associations with variants in genes that influence normal neuronal pattern development.The American Journal of Human Genetics,2011,89(2):334–343.[doi:10.1016/j.ajhg.2011.07.011]
4 Jia PL,Zheng SY,Long JR.dmGWAS:Dense module searching for genome-wide association studies in proteinprotein interaction networks.Bioinformatics,2011,27(1):95–102.[doi:10.1093/bioinformatics/btq615]
5 Collins FS,Morgan M,Patrinos A.The human genome project:Lessons from large-scale biology.Science,2003,300(5617):286–290.[doi:10.1126/science.1084564]
6 Yong Y,He L.SHEsis,a powerful software platform for analyses of linkage disequilibrium,haplotype construction,and genetic association at polymorphism loci.Cell Research,2005,15(2):97–98.[doi:10.1038/sj.cr.7290272]
7 Ghoussaini M,Fletcher O,Michailidou K.Genome-wide association analysis identifies three new breast cancer susceptibility loci.Nature Genetics,2012,44(3):312–318.[doi:10.1038/ng.1049]
8 Winham SJ,Cuellar-Barboza AB,Oliveros A.Genome-wide association study of bipolar disorder accounting for effect of body mass index identifies a new risk allele in TCF7L2.Molecular Psychiatry,2014,19(9):1010–1016.[doi:10.1038/mp.2013.159]
9 Hirokawa M,Morita H,Tajima T.A genome-wide association study identifies PLCL2 and AP3D1-DOT1LSF3A2 as new susceptibility loci for myocardial infarction in Japanese.European Journal of Human Genetics,2015,23(3):374–380.[doi:10.1038/ejhg.2014.110]
10 Goh KI,Cusick ME,Valle D.The human disease network.Proceedings of the National Academy of Sciences of the United States of America,2007,104(21):8685–8690.[doi:10.1073/pnas.0701361104]
11 汪小帆,李翔,陳關榮.復雜網(wǎng)絡理論及其應用.北京:清華大學出版社,2006:35–38.
12 賈華仟.復雜網(wǎng)絡分析方法在全基因組關聯(lián)研究中的應用[碩士學位論文].青島:山東科技大學,2015.