基于貝葉斯網絡的基因變異間的因果關系發現與驗證

蔡瑞初 甄啟祺 陳 薇 郝志峰,2

1(廣東工業大學計算機學院 廣東 廣州 510006)2(佛山科學技術學院數學與大數據學院 廣東 佛山 528000)

0 引 言

全基因組關聯研究(Genome-Wide Association Study，GWAS)當前主要側重于單核苷酸多態性(Single Nucleotide Polymorphisms，SNPs)與人類疾病的關聯研究[1]。研究表明，人類常見疾病的復雜性狀涉及多個基因變異之間的相互影響，而基因變異間存在連鎖不平衡(Linkage Disequilibrium，LD)現象[2]。由于基因按一定次序在染色體上呈線性排列，關于基因變異間的相互作用大多是從基因的相對物理位置來探討，然而在基因測序數據中，一些位于不同染色體上的基因，其測序數據都反映出了非偶然性的連鎖不平衡現象[3]，單純研究相鄰基因之間變異的相互作用會很受局限，故多基因變異相互影響是GWAS中不可回避的關鍵問題。

在多基因變異間的相互作用研究中，Yang等[4]提出了多變量全基因組關聯測試的方法，對多變量定量表現型數據進行分析。這是一種線性混合模型，實驗表明，該方法較好地把第一類錯誤控制在合理范圍內，并且比不同族群結構和相關性的邊緣檢驗更有效。然而該方法的適用條件為變量之間滿足線性關系，該假設較為苛刻。

事實上，由于多個基因變異之間的相互作用效果并不是簡單地線性加權，傳統的線性模型遠不能滿足這種多基因模型[5]。Monneret等[6]著重于邊緣因果估計方法，提出了基于高斯結構方程的高斯因果模型，然而該模型是一種線性模型，與真實情況中基因相互影響的非線性并不符合。

立足于非線性角度，Botta等[7]對隨機森林進行擴展，得到T樹模型。通過實驗與線性模型相比，T樹的預測能力表現出明顯優勢，并表明了多個SNPs之間的非線性效應。在進一步的基因座識別研究中，T樹模型不僅復現了已發現的多數基因座，并且發現了與克羅恩病相關的兩個新的符合生物學的易感位點。而在更復雜的非線性關系如網絡關系的研究中，T樹的效果有待驗證。Lim等[8]提出了BTR(BoolTraineR)算法，對布爾模型進行了優化，使用單細胞數據重構了基因調控網絡，然而沒有進一步考慮基因調控和基因變異之間的關系。

在基因變異和基因調控方面，現有的研究主要關注具體基因與疾病之間的關系。在基因變異上，Pan等[9]闡明了基因TNNT2其新變體R205Q的因果作用，Elsaadany等[10]通過全外顯子組測序研究了WWOX基因中的新純合突變，并提出了WWOX基因變異與人類癲癇性腦疾病之間的關系。基因調控方面，Wang等[11]研究了N6-甲基腺苷對脂肪的生成代謝起到了關鍵的調控作用，Jin等[12]研究了microRNA-192對膀胱癌細胞的調控，表明microRNA-192可能可以通過調控細胞周期，從而抑制膀胱癌細胞。Liu等[13]在Sia等[14]的基礎上構建了CNV-ICC-TRN，結合KEGG信號通路進行分析，發現基因變異的影響會發生在通路組件上。Liu等認為，整合基因變異和調控網絡二者的信息可更深層次地研究肝癌致病機理，然而目前基因變異和基因調控這兩者關系的相關成果相對較少，是當前亟需進行的研究。

鑒于上述研究狀況，針對部分不足，本文對基因變異間的因果關系進行研究，發現在同一信號通路中相互調控的基因，其變異間具有較強的因果性。本文通過改進的Parents-Children(Improved Parent-Children，IPC)算法，檢測出具有直接相互作用的基因對，并通過結構識別的方法判斷出因果方向。實驗結果表明，相互調控的基因變異之間具有較強因果性。本文還發現了一批具有較強因果性的基因變異可供相關研究參考。

1 相關知識

1.1 SNP與基因

單核苷酸多態性(Single-Nucleotide Polymorphism,SNP)指在基因組中某特定位點上的單核苷酸所發生的變異，該變異在人類群體中占有一定的規模比例[15]。在多代遺傳中，SNPs一般都會穩定地遺傳下去，故在GWAS中，一些跟人類疾病相關聯基因的變異，可以用SNPs作為穩定標記。

基因是一段DNA序列，由多個脫氧核糖核苷酸組成，是生物體的遺傳分子。人體的性狀由基因和生活環境控制調節，其中基因起著決定性作用。不同基因的DNA序列長度并不相同，而SNP是對單核苷酸狀態的描述，由此容易得知，基因與SNP之間是屬于一對多的關系。

1.2 信號通路

在多細胞生物體內，細胞外的細胞因子、激素等分子信號特異性地與細胞膜結合或直接穿過細胞膜，經過級聯轉導后多種信號在細胞內相互識別作用，共同調控內部的生化反應。細胞內多種信號的作用途徑交織成各種錯綜復雜的調控網絡，進行生化功能的調控，這樣的網絡稱為信號通路(Signal Pathway)[16]。京都基因與基因組百科全書通路(Kyoto Encyclopedia of Genes and Genomes Pathway，KEGG Pathway)是一個手工繪制出這些信號轉導網絡圖的集合。通過KEGG Pathway，本文可以把處于同一個信號通路的基因聚集起來，驗證處于同一信號通路中，相互調控的基因變異間因果關系的存在性。

1.3 貝葉斯網

貝葉斯網是一類概率圖模型。該模型借助一個有向無環圖，表示一系列的隨機變量X1,X2,…,Xn及其相互之間的依賴關系，其中每個節點代表一個隨機變量，并對應一個概率分布表。圖模型直觀地反映了各個隨機變量之間的依賴關系，而聯合概率分布表則具體量化了隨機變量之間的概率依賴強度。貝葉斯網如圖1所示，全部隨機變量的聯合概率分布可通過這些隨機變量的邊緣概率分布和條件概率分布相乘而得到，記為：

(1)

式中：Pa(Xi)為Xi的父節點集。

圖1 貝葉斯網模型

2 數據處理

2.1 數據描述

實驗所用的SNPs數據由英國WTCCC(Wellcome Trust Case Control Consortium)機構提供，本文所使用的數據是多種復雜疾病整合到一起的數據文件，共有16 179個樣本，每個樣本含有394 747個SNPs的變異情況。表1選取了其中兩個SNPs的染色體位點和主次等位基因等基本屬性的信息，部分SNPs變異情況如表2所示。

表1 每個SNP的染色體位點及其基本屬性

表2 部分SNPs變異數據

表2列出了數據集中3個樣本的3個SNPs情況，每一列為一個SNP的變異情況，如表2中的rs2980300_T，這是一個命名為rs2980300的SNPs，其后綴“_T”代表該SNP位點突變為次等位基因T；0、1、2三個取值分別代表樣本1在該SNP位點上沒有發生突變，樣本2有1條染色體上的該SNP位點突變為此等位基因T，樣本3則2條染色體全都發生了突變。

2.2 SNPs數據編碼為基因數據

通過NCBI可以查詢到rs命名的SNPs是否位于基因區域、屬于哪一個基因。根據基因和SNPs為一對多的關系，SNPs數據可以編碼成基因粒度，編碼方式如表3所示，394 747個SNPs中有3個位于基因51150區域。

表3 基因51150的編碼情況表

3 基于IPC算法的基因變異因果關系發現

3.1 IPC算法框架

基于所得的離散型基因變異數據構建貝葉斯網模型，如果檢測到一個基因變量的父子(Parents and Children，PC)節點，則可以推斷出調控網絡中直接相互作用的基因。通過條件獨立性檢驗構造局部因果結構，可以檢測出貝葉斯網中基因變量的PC節點[17-19]。受貝葉斯半監督方法(Bayesian Semi-Supervised method，BASSUM)[17]中采用PC算法檢測變量間因果關系思想的啟發，本節提出一種改進的PC算法，即IPC算法(Improved-PC)，檢測同一信號通路中基因變異間相互作用的因果關系。

記PC(T)為目標基因T的PC節點集。由貝葉斯網性質可知，PC(T)是導致基因T發生變異的直接原因或由基因T變異導致的直接結果的集合，即PC(T)中所有基因變量跟T都是不獨立的。記vi為跟基因T直接相互作用的基因，由此可得vi∈PC(T)的一個必要條件：vi跟T不滿足相互獨立。根據此條件，在檢測目標基因T的PC節點時，遍歷PC(T)的候選基因集合C(T)，如果屬于C(T)的基因ci與T不相互獨立，則更新PC(T)，使得PC(T)=PC(T)∪ci，這里ci∈C(T)。由于上述條件是vi∈PC(T)的必要條件，據此檢測所得的結果很可能會引入一些非直接相互作用的基因，得到的是PC(T)的一個超集，故檢測結束后需剔除掉非直接相互作用的基因。

先對超集中的非直接相互作用的基因分析。圖2所示的PC(T)超集包含了非直接相互作用的基因fi，基因fi發生變異是目標基因T發生變異的間接原因，檢測時如果基因fi被誤加入PC(T)中，那么從貝葉斯網絡結構上看，節點fi就不可能跟節點T直接連接。同時，由于根據必要條件檢測所得的集合是真實PC(T)的超集，故fi跟T必然被PC(T)所d-分隔，從而在給定PC(T)的條件下，fi跟T條件獨立，其中PC(T)可以通過遍歷超集的冪集搜索獲得。此分析過程對目標基因T變異的間接結果同樣適用。

圖2 根據必要條件檢測到的PC(T)超集

類似地，如果基因vi屬于真實PC(T)，則基因vi跟目標基因T直接相互作用，故與fi情況相反，基因vi和目標基因T不會被任何集合所d-分隔，即不存在條件集合使vi跟T滿足條件獨立。根據vi和fi在條件獨立性上的差異，可以把超集中的非直接相互作用的基因fi和屬于真實的PC(T)基因vi分辨出來。

根據以上分析，可以得到IPC算法步驟如算法1所示。

算法1IPC算法

輸入：包含所有變量的集合V、目標變量T

輸出：PC(T)

1. 初始化：PC(T)=?，C(T)=V-{T}；

2.FORvi∈C(T)-PC(T)DO；

//遍歷；

3.g(vi)=G2(vi,T)；

//構造基因vi和

//目標基因T的

//G2統計量

//如果G2統計量落在拒絕域中，

//說明基因vi和目標基因T

//不滿足相互獨立

5.PC(T)=PC(T)∪vi；

6.ENDIF

7.ENDFOR

//通過條件獨立性檢驗刪除冗余信息

8.FORfi∈PC(T)DO

9.IF?S?PC(T)-{fi}，fi⊥T|STHEN

10.PC(T)=PC(T)-{fi}；

//遍歷PC(T)減去fi后的冪集，以此

//作為條件進行條件獨立性檢驗

11.ENDIF

12.ENDFOR

IPC算法檢測到基因的PC節點后，需要通過因果結構的識別來判定PC基因之間的因果方向。由于直接確定兩個相互連接的PC節點間的方向比較困難，因此，本文通過對三個節點的基本有向無環結構進行識別，進而判斷其中的父節點(Parent)和子節點(Child)。三個節點構成的基本結構包括圖3的4種情況。

(a) 順連 (b) 順連

(c) 分連 (d) 匯連圖3 三個節點構成的基本結構

這4種情況可以從獨立性和條件獨立性上的差異來實現結構的識別。情況(a)-(c)中，節點vi和節點vj并不滿足相互獨立，而在給定節點T這一條件下，節點vi和節點vj滿足條件獨立。這三種結構在獨立性和條件獨立性的性質上等價，無法進行有效區分，但情況(d)會呈現出相反的性質。

情況(d)屬于匯連結構。該結構跟前面三種結構的區別在于，在節點T未知的情況下，節點vi和節點vj滿足相互獨立，但在給定節點T的條件下，節點vi和節點vj不滿足條件獨立。這跟前面三種結構的性質相反。因此，利用匯連結構在條件獨立性上的特性，本文可以通過算法2識別出三節點構成的匯連結構，推斷其中的父親節點和孩子節點。

記從節點vi出發，到節點T終止的有向邊為，具體的識別過程如算法2所示。

算法2匯連結構識別算法

輸入：目標變量T、T的PC節點集PC(T)

輸出：包含目標變量T的有向邊集E

1. 初始化：E=?；

2.FORvi,vj∈PC(T)DO

3.IFvi⊥vjANDNOTvi⊥vj|TTHEN

4.E=E∪{,}

5.ENDIF

6.ENDFOR

不同于一般的相關性分析，在目標變量T和PC節點集PC(T)中，算法2基于條件獨立性檢驗，先有效識別出變量中的匯連結構，進而推斷出該結構中的因果方向，保證了檢測得到的因果關系的可靠性。

3.2 獨立性檢驗

本文所采用的獨立性檢驗基于卡方檢驗實現。對一些較為復雜而龐大的數據進行獨立性檢驗時，平方差的計算會增加卡方檢驗的運算量，所以可以構造G2檢驗統計量[20]來作為卡方統計量的代替。

可借助列聯表來考察兩種屬性的獨立性。如表4所示，Omn表示樣本中屬性1取值為Am且屬性2取值為Bn的觀察計數。其中M為屬性1的取值數目，N為屬性2的取值數目。

表4 M×N列聯表

如果零假設“屬性1和屬性2相互獨立”成立，則兩者的聯合分布率等于邊緣分布率之積。故聯合分布率的極大似然估計為：

(2)

由零假設和式(1)可知，樣本數據中屬性1取值為m且屬性2取值為n的期望樣本量為:

(3)

而在備擇假設中，兩屬性的聯合分布率為:

(4)

類似地，在進行條件獨立性檢驗時，如果在給定屬性3取值為l的條件下，檢驗屬性1和屬性2是否相互獨立時，對數據中滿足給定屬性3的取值l的樣本構造屬性1和屬性2的列聯表，由以上推導可得期望樣本量為：

(5)

式中：各項的右上標l代表給定條件屬性的取值為l。

在M×N列聯表中，任意單元格中的觀察樣本數都對應著唯一一個期望樣本數，為了便于記錄，這里把M×N個單元格的觀察樣本數記為O1，O2，…，OM×N，期望樣本數則對應為E1，E2，…，EM×N。于是可以得到G2檢驗統計量作為卡方檢驗的近似估計，表示為:

(6)

另一方面是自由度的計算。對每個列聯表的G2檢驗統計量其所對應的自由度df(degrees of freedom)為:

df=(M-1)×(N-1)

(7)

在給定條件屬性時，自由度為:

df=(M-1)×(N-1)×S(Con)

(8)

式中：S(Con)表示條件屬性值的數目。

4 實 驗

4.1 實驗方案

要驗證同一信號通路中相互調控的基因，其變異間具有因果關系這一假設，需要通過IPC算法檢測出具有因果關系的基因對，并將該實驗結果分別跟真實信號通路中相互調控的基因對進行匹配。本次實驗的卡方獨立性檢驗顯著性水平設為0.01。

4.2 實驗結果分析

本文使用了KEGG Pathway的真實信號通路數據來評估實驗結果。為了評估實驗結果的匹配程度和IPC算法的性能，本文使用召回率R(Recall)、準確率P(Precision)和F1值(F1-sorce)來評價實驗效果。三者的定義如下：

(9)

(10)

(11)

式中：TP代表IPC算法檢測出來的基因變異對和真實信號通路網絡結構中相互調控的基因對匹配的數目；FN是真實信號通路網絡結構中含有而IPC算法沒有檢測出的基因變異對的數目；FP為IPC算法檢測出而真實信號通路網絡結構中沒有出現的基因變異對的數目。

根據上述指標，最后得到的實驗結果與真實信號通路的對比效果如圖4所示。

(a) 召回率

(b) 準確率

(c) F1值圖4 實驗檢測結果跟真實信號通路網絡結構的匹配評分

可以看出，算法的部分檢測結果與真實信號通路中相互調控的基因對相匹配，準確率比較理想，其中某個包含161個基因的信號通路中，其準確率達到90%以上。該實驗結果表明，部分相互調控的基因之間的變異具有一定的因果關系。

為了進一步驗證這種匹配的非偶然性和IPC算法的有效性，本文進一步進行了隨機對照實驗作為對比實驗。實驗保持真實信號通路的拓撲結構不變，將同一信號通路中的所有節點所對應的基因標簽隨機混淆，信號通路中原有的調控信息破壞掉，并以隨機信息替換，得到了隨機對照的網絡結構。結果發現其中并沒有基因對可以與IPC算法結果相匹配，部分實驗結果如表5所示。這一結果驗證了這種匹配的非偶然性，基因變異間的因果關系跟信號通路中的基因調控關系存在某種關聯。

表5 本文IPC算法實驗、LD算法對比實驗和隨機對照實驗的指標對比

續表5

在圖4中一個具有161個基因的信號通路，其召回率、準確率和F1值都突然明顯升高，其中原因可能跟信號通路的類型有關。進一步分析，基因數目超過10個的信號通路共有181個，總體情況如表6所示。其中，細胞過程和環境信息處理這兩類信號通路的平均基因數較多，其平均準確率也偏低；遺傳信息處理類型的信號通路只有4個，代表性較弱；人類疾病和生物體系統這兩類的平均基因數都在48左右，然而平均準確率分別是0.11和0.18，相差較大，可能還有其他因素本實驗中沒有考慮到；實驗效果比較理想的是新陳代謝類型的信號通路，相較于人類疾病和生物體系統類型，雖然其平均基因數較少，但平均準確率增幅很大，達到0.38，且樣本量有40，具有一定的代表性。這也說明了對于某些信號通路，通過本文算法可以有效地檢測出其基因變異間所蘊含的因果關系。

表6 基因數超過10個的信號通路的情況

表7給出了部分跟真實信號通路網絡結構相匹配的IPC算法檢測結果，調控類型列中ECrel代表酶與酶的催化逐次反應，PPrel代表該基因對的相互作用表現為蛋白質與蛋白質的相互作用。

表7 檢測基因變異的因果關系所驗證到的基因的調控情況

5 結 語

本文研究了相互調控的基因之間的變異的因果關系。把WTCCC研究機構的SNPs數據編碼為基因粒度的數據，并使用因果關系發現IPC算法檢測出了數據中的因果關系。將檢測結果分別跟KEGG Pathway的真實信號通路網絡結構數據和隨機對照實驗網絡結構對比，發現因果結構在部分真實信號通路網絡中獲得了驗證，而跟基于連鎖不平衡的對比實驗則體現了IPC算法具備更良好的性能。