基于R語言的互信息網絡模型在乳腺癌易感基因檢測分析中的應用①

王淑棟,張善強,賀思程

(中國石油大學(華東)計算機與通信工程學院,青島 266580)

隨著人類基因組測序工作的逐步完成,大量的數據為全基因組關聯分析提供了豐富的素材,也涌現出許多數據分析方法[1-4].人類基因組計劃得出人類所有的基因共由39 000多個已經編碼蛋白的基因序列以及30億堿基組成.而國際單體型圖計劃[5]得到了SNP的300萬個位點.兩個計劃的實施給生物學領域帶來了眾多的數據信息,為全基因組研究中提供了方便.GWAS因其優勢得到了很多的應用.大量研究成果顯示關聯研究具有很多的優勢[6].

Ghoussaini等[7]在2012年針對乳腺癌相關基因進行研究,共得到了3個致病相關的位點,rs10771399不僅在乳腺癌的發展中起著關鍵作用,在骨轉移中也有著同樣的重要性.2013年,維爾漢姆等[8]關于躁郁癥的數據進行分析,得出與躁郁癥相關的SNP位點及致病基因.2014年,廣川等[9,10]針對心肌梗塞病設計了病例對照實驗,從實驗中得到了有關疾病的致病基因和SNP,使心肌梗塞病得到了合理的解釋.

GWAS能夠幫助人們更好的解釋復雜疾病成因,但是它也有不足.一方面,復雜疾病多種多樣,其中的影響因素也很多,如何確切地得到與特定的功能相聯系的位點是個不小的難題;另一方面,對于GWAS結果,它在不同群體中的影響程度并不一樣;目前的大部分研究主要針對簡單疾病,沒有涉及到基因間的相互作用.

而針對基因間的相互作用,可以通過互信息建立網絡進行表達.GWAS網絡方法將GWAS數據進行網絡建模,通過比較疾病數據與對照數據得出的網絡的不同,進行后續的相關統計量的分析及解釋.

本文試圖通過互信息表示SNP之間的相互作用關系,進而建立SNP與SNP之間的網絡.在此基礎上,進行全基因組關聯研究,找到結構性關鍵SNPs.

1 互信息網絡建模

隨著生物網絡的研究深入發展,研究者對元素之間的相關性的描述越來越準確,互信息作為兩個元素之間的相關信息度量,具有很多的優勢,其中最具優勢的就是它的熵表示,不僅是對元素出現概率的表示,更是體現了元素之間的離散程度及相互之間的關系,對于給定的兩個SNP表達序列,他們之間的數據存在著差異,而利用互信息可以充分表達SNP之間的差異性及依賴性,互信息越大,說明兩個SNP之間的關聯程度越緊密;反之,則說明聯系越小,從而找到跟所有的SNP聯系較大的節點,即是關鍵SNP.本文通過互信息建立相互作用網絡,從而分析網絡結構的差異性.設是兩個SNP的基因型數據在個體之間表達形成的向量,是X和Y的聯合概率分布,是他們之間的聯合熵,定義為:

對于兩個隨機變量之間存在的關系,H(X)表示隨機變量X蘊含的不確定性,而條件熵則是已知條件Y時隨機變量X所余下的不確定性,那樣,就表示已知條件Y后X包含的信息量.進而還可以證明這個值關于X和Y是對稱的,即且都等于由此X和Y之間的互信息可以計算,互信息記為

因為SNP數據是每個SNP仿真1000組得到的數據,每三個數據代表一個個體,首先需要對數據進行處理使得數據能夠表示基因型,我們確定使用0,1,2三個數來表示每個個體內表達的基因型,再根據公式(2)計算得到所有的SNP之間的互信息.具體計算過程如下:

(1)我們首先得到每個SNP的基因型可能性序列數據,假設共有N個個體,則每一行包含2N個SNP堿基可能性數據,0代表出現,1代表不出現.

例如:假定兩個個體關于5個SNPs的基因型數據如下:

SNP 1:AA AA

SNP 2:GG GT

SNP 3:CC CT

SNP 4:CT CT

SNP 5:AG GG

輸出的正確仿真數據如下所示:

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出現的可能性序列分別是100和010.

(2)每個SNP的基因型表達數據作為一個向量,x,y表示來自SNP集合I中的其中的兩個SNP向量.

(3)根據每個SNP的基因型表達量的分布,計算得到每兩個SNP之間存在的互信息值.所有SNP之間的互信息構成互信息矩陣,記作矩陣中的每行代表一個SNP,每一列代表此SNP與另一個SNP之間的互信息.

假定存在一個集合的SNP基因型數據D,其中所擁有的SNP的集合我們記作I,可由互信息計算公式(2)得到一個互信息矩陣.定義一個建立在關于SNP基因型數據D的互信息網絡.是邊賦權圖,其中V表示點集合、每個網絡中的節點i∈V表示一個SNP,基因i和j之間的互信息計算值wij定義為每條邊的權重.在下面的表述中,我們將基因i∈I以及頂點i∈V等同起來看待.

2 基于網絡統計量的關鍵基因選取

利用上述方法得到的SNP相關網絡中各節點(SNP)的網絡結構參數來描述特定生物過程中基因的重要性.首先給出幾個重要的能夠反映網絡結構特點的網絡統計量的相關定義[11].

(1)度(K):在網絡中,度指的是與該點相連接的邊數目.節點度可以表示該點的重要程度,節點度越大,表示該點在網絡中越重要.而網絡的平均度可以通過計算所有的點的度,后取平均數計算得到.

(2)平均路徑長度(L):定義為網絡中所有的點之間兩兩求得的距離的平均數,網絡中的任意兩點i,j的距離即邊的條數,則兩點之間的平均路徑長度表示為所有的點之間的平均距離,記作:其中N表示網絡中的節點數目.

(3)聚類系數(C):網絡中節點i有Ki個邊與之連接,那么與該點可能連接的最大邊數為若這Ki個節點之間真實邊為Ei,則它與總的所有情況下的邊比例,計算得到節點i的聚類系數很顯然,0≤C≤1.C=0代表網絡中的點為孤立點;C=1表示網絡中的所有點之間都是互相連接的,視為全局耦合網絡.

(4)介數(B):網絡中介數的概念可以分為兩類,一類是點介數,另一類是邊介數.節點k的介數定義為,其中,C(i,j)代表i與j間最短路徑總數,表示中間點為k時,i與j間的所有路徑總數.介數反映了節點k在i和j之間的流通量和重要程度.網絡中某個節點的介數越大,說明該點在網絡中信息傳播的信息量就越大,越容易在該點造成網絡堵塞.假設兩組連接度很高的網絡中間只有少數點連接,那么這幾個少數點介數就會很大,即很多的信息在流通的過程中經過這幾個點,很容易造成堵塞,從而造成數據信息丟失.因此,最大介數的增大會降低網絡同步能力.

(5)模塊度(Q):模塊度也稱作模塊化度量值,是用來衡量網絡強度的統計量.最早是Newman提出的,它用來描述網絡社團以及劃分的好壞.假定網絡共分為k個社團,代表一個k×k維的矩陣.故模塊度可以定義為:其中,是矩陣中的數值之和(行或列),eij用來表示社區i和社區j之間的邊的數量.模塊度可以區分社區劃分的好壞.若是劃分的好,則社區內部節點相似度較大,而在社區外邊相似度較低.Q越大,越接近1,代表社區擁有一個很好的劃分結構,使得社區的劃分合理化.通常設定的值是在0.3與0.7之間.

本文中我們主要選擇5個參數進行分析比較,對于給定的參數進行最終的分析,從而找到影響網絡的重要因素,依據此類統計量進行歸納分析,得出相應的參數.

我們對由SNP數據設定不同的互信息閾值而形成網絡,針對其中大于閾值的邊,做去掉處理,而針對小于閾值的邊進行保留操作,從網絡圖可以分析出統計量所對應的參數變化,得到有益信息量.

根據網絡中SNP之間互信息計算的值,選擇閾值范圍為0.1到0.63.共設置63個閾值,在每個閾值的條件下,統計計算相應的網絡結果,從而得到一致性網絡,根據網絡的相似性程度選擇對實驗組和對照組差別較大的統計量進行分析.我們最終選擇了度作為區分依據,并分析能夠區分實驗組和對照組的取值范圍,得出最佳的閾值,對于不同的數據,得到的互信息值也不同,所以需要根據數據得到的互信息范圍,由網絡統計量得到取值范圍,得到互信息取值的交集,能夠區分對照組和實驗組數據,從而確定最佳的互信息閾值.這樣就能夠保證所取的閾值不受樣本數量的大小影響,而是根據樣本的不同情況得到相應的閾值.對于節點i,我們定義,?d代表了這個節點的度差異值,在該公式中,分別代表了這個節點在實驗組與對照組網絡中節點的度.

我們都知道,在復雜網絡中,節點度能夠代表節點的作用和影響力.本文從網絡結構差異的角度去衡量各個統計量[12],進而對應到其中的節點,找到“結構性關鍵SNPs”.這種差異性貢獻分為正、負貢獻兩個方面.我們用r代表度的變化閾值.正貢獻SNP代表了該節點在病例組、對照組兩個網絡中度的貢獻?d≥r的SNP;同理,負貢獻SNP代表了該節點在以上兩個網絡中度的貢獻?d≤-r的SNP.

本文對基因BRCA2仿真數據建立病例組與對照組建立相互作用網絡進行數據實驗.對SNP互信息設置一個閾值范圍,分析產生的病例組和對照組SNPs互信息網絡的統計量:平均路徑長度、聚類系數、平均度、模塊度、平均介數隨閾值在其變化范圍內的增加而變化的情況.根據計算的網絡中SNP之間互信息的值,我們取互信息閾值的范圍為0至0.63,步長0.01,分析對應病例組與對照組的SNP相互作用網絡的上述網絡結構參數隨變化而變化的情況.

3 數據來源與處理

HapMap給出了人類基因組單核苷酸多態性(SNPs)和拷貝數多態性(CNPs)的分布情況.本文使用HapMap提供的三個文件進行實驗,包含了關于BRCA2的88個SNPs.下面是對三個文件的說明.

.hap文件是已知的單體型數據,其中行代表SNP,列表示單體型.每一個.hap文件都需要一個相應的legend文件,所有的等位基因都以0,1作為標記.

.legend文件是SNP標記位點數據,四列數據分別表示SNP的ID、堿基位置、堿基的0,1表示.

.map文件包含了小規模的重組率,共三列分別表示每個SNP的物理位置,距離左標記點的位置和距離右標記點的位置.

在這數據中,必須去掉全部為0或者全部為1的數據,因為這些數據對構建網絡結構沒有任何幫助.去掉這些多余的數據,共得到45條SNP數據.把3個文件放到一起,執行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組實驗組和對照組數據.隨機選定2個SNPs作為致病SNPs.它們的信息如下:rs206081和rs9534318,選取雜合子變異率分別是2.5和1.5,純合子變異率分別為2和4.5,上述樣本數據都包含SNP編號,SNP位置及0,1表達數據.

本文中,我們使用.gen文件,刪除前五列后把數據轉換成一個矩陣,其中每行表示一個向量,每三個數字代表一個個體,我們轉換成0,1,2表示.

4 試驗方法和結果

4.1 網絡統計量的比較

根據得到的互信息矩陣,大于閾值的向量之間表示相互關系較強,選定這些SNP作為節點建立網絡.分析比較網絡的6個特性.每個結構參數都反映著網絡的特性,進而可以顯示SNP間的互信息的變化,取0.01為步長,從0到0.63之間求得每一個閾值下的網絡結構特性值,得到圖1.圖1中,縱坐標表示相應的統計量,橫坐標代表閾值,虛線表示對照組數據顯示效果,實線表示實驗組數據顯示效果.

圖1 4個網絡結構的統計量隨閾值的增加的變化情況

實驗發現5個結構特性中,平均聚類系數B交織在一起,不能區分實驗組和對照組.

觀察圖1(a),當0＜t＜0.21時,網絡的平均介數B在在兩組中的變化趨勢走向大體相似.當0.21＜t＜0.63時,網絡的平均介數B逐漸減小.從圖中可以明顯的看出,病例組的平均介數要比對照組的平均介數高.于是,我們得到,隨著互信息閾值的增大,節點的介數也在不斷減小,網絡中邊越來越稀疏.

觀察圖1(b),當0.2＜t＜0.43時,實驗組與對照組的網絡有相對明顯的差異.于是我們可以得到,在這個變化區間內,平均路徑長度可以很好的區分病例組和對照組,而當t＞0.43時,網絡的邊越來越少,平均路徑長度趨近于0.

從模塊度Q隨閾值的變化圖1(c)看出,當閾值0＜t＜0.2 或 0.43＜t＜0.63 時,兩組中的模塊度Q逐步上升,但變化大致相同,而當 0.2＜t＜0.43 時,實驗組模塊度與對照組有較大區別.

觀察圖1(d),可以發現,在很長的一段閾值范圍內,病例組與對照組的網絡平均度有很大的區別,而隨著網絡的閾值增加,網絡的平均度越來越小,這與網絡的孤立點越來越多也是相對應的.

當t＞0.62時,病例、對照組中都只有一個包含四個節點的全耦合子網,聚類系數C、平均路徑長度L兩者相等,且都為1.當t＞0.63時,平均路徑長度L、聚類系數C是缺失的,平均介數B以及其他三個統計量值均為0.

總之,平均聚類系數C不能區分兩組數據,平均路徑長度L和平均介數B能夠區分但是閾值具有一定局限性.平均度可以在很大的范圍內把實驗組和對照組分別出來,我們選擇平均度作為區分的依據.

從圖1中我們得到每個統計量能夠區分兩組的閾值范圍,如表1.

表1 各統計量能夠區分實驗組和對照組的閾值范圍

從表1可以看出,每一個統計量都有不同的閾值范圍,平均度K的范圍較大,0.08＜K＜0.35;其他的統計量閾值范圍相差不大,基本在0.2到0.3之間.結合圖1,選擇0.28為閾值構建網絡.

依據圖2,實驗組和對照組的圖像是有很大差異的.在對照組,節點之間聯系較弱且存在更多的孤立點.但是在實驗組中,很多的孤立點不再是獨立的,并且擁有了更多的聯系.對照組中存在36個連接點和9個孤立點,而實驗組中存在39個連接點喝6個孤立點.這表明我們選取的閾值0.28是合適的.經過多次仿真數據試驗,對于結合數據互信息得到閾值范圍,而后確定互信息閾值的方法都是有效的.

圖2 閾值為0.28的條件下,實驗組和對照組互信息網絡

4.2 獲得“結構性關鍵SNP”

結構決定功能,而結構的差異決定了功能的差異,本文將這種差異細化到每個節點上,而平均度可以很好的區分病例組和病例組,所以我們選擇每個SNP位點的平均度來刻畫SNP在病例組和對照組的差異,計算每個網絡的每個節點的節點度差異,當節點的度在病例、對照組中的變化差異比較大時,說明這兩個組的網絡結構差異較大.從兩組網絡的數據分析來說,節點度的增量有正有負,所以,節點在病例組中的度也有增減之分,即存在正、負貢獻SNPs.度變化量增加最大的是節點39,增加值的大小是5,同理,減少量最大的是16,41,減少值的大小是2.

當閾值為0.28時,對照組網絡中的平均度大致等于2,從而可以得到,當病例、對照組網絡中節點度的變化值大于等于3時,其對網絡結構影響較大.故可設?d=3,由此,我們可以獲得對網絡結構有顯著影響4個SNPs,如表2,其中rs206081,rs9534318為預設致病SNPs.

表2 給定參數為3的條件下,部分結構性關鍵SNPs的信息及度的變化量

4.3 參數評估

在查找“結構性關鍵SNPs”時,我們需要從網絡平均度出發,對選取網絡中的關鍵SNPs設置合適的差值參數.如果選取的差值參數比較小,對SNPs選取限制比較寬泛,一些不相關的SNPs也會選取到SNPs集合內,從而導致假陽性.反之,如果選取過于嚴苛,反而會遺漏一些比較重要的節點,導致假陰性.

我們選取基因BRCA2,得到它在閾值為0.28時候的網絡,如圖2所示.選擇不同的差值參數,得到一系列不同的結構性關鍵SNPs,如表3所示.

表3 不同參數r的取值下關鍵SNPs個數

當互信息閾值設定為0.28時,網絡中度的最大變化量是5.當r≥5時,所得的關鍵SNPs只有節點39,對網絡影響較大的節點25卻被忽略.當r≤2時,所得的關鍵SNPs只有13個,這里面也包括了其中的非零點.

5 結論與分析

本文通過國際項目HapMap3中以及Hapgen2軟件生成的13號染色體上BRCA2基因生成仿真數據.利用互信息表示SNPs間的相互作用.構建實驗組和對照組的網絡,根據閾值及差值參數篩選出關鍵SNPs.最后,對我們所選擇的參數進行了評估,證明我們所選定的參數能夠反映結構的變化,能夠較好地選擇出預設的關鍵SNPs.通過數值實驗發現:樣本數目會影響互信息的大小,樣本數較小時,互信息較高,樣本數較大時,互信息逐漸降低,本文認為,樣本數偏少,則特異性個體數目不完備,樣本數過多,又會造成冗余,增加了計算復雜度.目前,確定合適的上下界仍然是一個具有挑戰的問題.

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 汪小帆,李翔,陳關榮.復雜網絡理論及其應用.北京:清華大學出版社,2006:35–38.

12 賈華仟.復雜網絡分析方法在全基因組關聯研究中的應用[碩士學位論文].青島:山東科技大學,2015.