基于條件互信息的乳腺癌易感基因網(wǎng)絡的構建和分析?

王淑棟 范曉丹 王新贈

（1.中國石油大學（華東）計算機與通信工程學院 青島 266580）（2.山東科技大學數(shù)學與系統(tǒng)科學學院 青島 266590）

1 引言

隨著基因芯片技術和高通量測序技術的發(fā)展，產(chǎn)生的大量數(shù)據(jù)為全基因組關聯(lián)研究（GWAS）提供了豐富的素材，期間也出現(xiàn)了許多數(shù)據(jù)處理方法［1～4］。近年來，大量研究成果顯示GWAS 具有很多優(yōu)勢：2014 年，Hirokawa 等［5］利用病例組和對照組數(shù)據(jù)對心肌梗塞疾病做了全基因組關聯(lián)研究，并確定了兩個新的與心肌梗塞發(fā)病機理相關易感性位點：PLCL2 和AP3D1-DOT1L-SF3A2。2016 年，Direk 等［6］通過薈萃分析先前兩個GWAS 研究的結果發(fā)現(xiàn)，位于FHIT 內(nèi)含子區(qū)域的一個新的抑郁癥狀相關的位點（rs9825823，P=1.0*10-9）。

從SNP 數(shù)據(jù)出發(fā)，度量SNP 間的相關性，并構建SNP-SNP 相互作用網(wǎng)絡，可有效挖掘SNP 間的關系，進而從生物分子網(wǎng)絡的角度認識生命現(xiàn)象并揭示生命活動的基本規(guī)律，有助于預測未知SNP功能、認識疾病發(fā)病機理、加速藥物開發(fā)等。隨著對生物網(wǎng)絡［7～9］研究的深入，對元素間相關性的度量方法也越來越多，傳統(tǒng)上主要有皮爾遜相關系數(shù)、斯皮爾曼相關系數(shù)等，被廣泛用于測量變量間的線性關系，但無法區(qū)分間接關聯(lián)和直接關聯(lián)。偏相關性（PC）由于可以檢測變量間的直接關聯(lián)，被廣泛使用，Barzel等［10］應用PC指標構建了一種動態(tài)相關性基因調(diào)控網(wǎng)絡，消除了基因間的間接影響，能有效區(qū)分基因間的直接調(diào)控和間接調(diào)控。然而，基于PC 的方法忽略了非線性系統(tǒng)（如生物分子網(wǎng)絡）中起重要作用的非線性相關性，因此近年來，互信息（MI）和條件互信息被廣泛應用于線性和非線性關聯(lián)的量化中。但MI 不能檢測直接關聯(lián)或依賴關系，且具有高估問題。CMI可以量化變量間的非線性直接依賴關系，優(yōu)于PC 和MI，因此被廣泛應用于許多領域［11～12］進行網(wǎng)絡直接依賴的推斷。

目前很大一部分GWAS主要針對簡單疾病，且很少涉及SNP間非線性直接依賴關系，如何準確定位疾病相關的SNPs仍是個不小的難題。本文針對基于MI 構建SNP-SNP 相互作用網(wǎng)絡假陽性邊偏高的問題，通過CMI 表示SNP 間的相關性，將乳腺癌相關的SNP數(shù)據(jù)進行網(wǎng)絡建模，進行全基因組關聯(lián)研究及節(jié)點網(wǎng)絡中心性的分析解釋，最終找到可能的致病SNPs。

2 數(shù)據(jù)來源與處理

本文使用了HapMap3 中位于13 號染色體上的包含88 個SNPs 的BRCA2 基因數(shù)據(jù)，包含.leg 文件、.hap文件及.map文件。

為了保證構建的SNP-SNP 相互作用網(wǎng)絡更具代表性，需要刪除意義不大的數(shù)據(jù)，去掉.hap 文件中全部為0 或全部為1 的數(shù)據(jù)，得到45 條SNP 數(shù)據(jù)。利用以上3 個文件，使用HAPGEN2 進行數(shù)據(jù)仿真：隨機選定rs9534318 和rs9943876 作為致病SNPs，設定對照組和病例組的雜合子變異率分別是1.5 和2，純合子變異率分別是2.25 和4，分別仿真1000 組病例組和對照組數(shù)據(jù)。接下來刪除仿真產(chǎn)生的.gen 文件中的SNP 的ID、名稱、堿基位置及等位基因信息，并把剩余數(shù)據(jù)轉(zhuǎn)換成45 行3000 列的矩陣，每行表示一個SNP 向量，每3 個數(shù)字代表一個個體。為了后續(xù)操作方便，按照100轉(zhuǎn)換為0，010 轉(zhuǎn)換為1，001 轉(zhuǎn)換為2 的規(guī)律處理該矩陣，分別得到新的1000 個個體的病例組和對照組SNP 基因型數(shù)據(jù)D1和D2。

3 基于CMI 的SNP-SNP 相互作用網(wǎng)絡的構建

假設X 和Y 是兩個隨機變量，互信息代表使用Y 編碼X 時所需的信息，反之亦然，即變量X和Y 間的相關性可用MI( )X;Y 度量。MI 是在KL距離D［13］的基礎上定義的：

式中，p( x )表示變量X 為x 時的概率值，p( y )表示變量Y 為y 時的概率值，p(x，y)表示變量X 和Y 分別為x 和y 時的聯(lián)合概率值。MI是根據(jù)X 和Y 之間的相互獨立評估的，定義如下：

條件互信息表示兩個變量在第3 個變量下的條件依賴性，能夠量化變量間的非線性直接關系，變量 X 和Y 在變量 Z 下的條件互信息CMI(X;Y|Z)定義如下［14］：

式中，p( z )表示變量Z 為z 時的概率值，p( x|z )和p( y|z )分別表示變量X 和Y 在Z 條件下的概率，p(x，y|z) 表示變量X 和Y 在Z 條件下的聯(lián)合概率，p( x，y，z )表示變量X 、Y 和Z 的聯(lián)合概率。CMI是根據(jù)變量X 和Y 在變量Z 下的條件獨立性評估的，定義如下：

如果變量X 和Y 在變量Z 條件下相互獨立，則CMI(X;Y|Z) 為零；CMI(X;Y|Z) 越大，表明X和Y 的相關程度越大。本文基于CMI 構建SNP-SNP相互作用網(wǎng)絡時，CMI(X;Y|Z)表達了兩個SNPs 在第三個SNP 下的相互依賴程度，CMI(X;Y|Z)越大，說明X 和Y 兩個SNP間的關聯(lián)程度越緊密。

對于SNP 基因型數(shù)據(jù)為D ，我們假定其SNP集合為I={1 ，2，…i，…n} ，根據(jù)CMI 式（4）可得CMI 矩陣CONM={C MIij}n*n(|I | =n )。并定義關于D 的CMI 網(wǎng)絡為G[ D] =(V ，E;w )，G 是邊賦權圖，其中V 表示點集合，E 表示邊集合，節(jié)點i ∈V 表示SNP i，對于?i，j ∈V ，節(jié)點i 和j 間的CMI 計算值wij定義為網(wǎng)絡中的邊( i，j )∈E 的權重。

對于數(shù)據(jù)處理后得到的病例組SNP 基因型數(shù)據(jù)D1，我們將其擁有的SNP 基因型表達數(shù)據(jù)的集合記作I1。計算每兩個SNPs間的CMI值wij，得到關于D1的CMI 矩陣CONM1，每行代表一個SNP，每列代表此SNP 與另一個SNP 間的CMI 值，將CONM1的對角線及下三角元素設為0，并構建基于CMI 的病例組SNP-SNP 相互作用網(wǎng)絡G[ D1] 。對對照組SNP 基因型數(shù)據(jù)D2進行相同處理，得到CMI 矩陣CONM2及對照SNP-SNP 相互作用網(wǎng)絡G[ D2] 。

4 最佳CMI閾值的選取

本文中我們選擇平均度和平均介數(shù)兩個網(wǎng)絡統(tǒng)計量的參數(shù)進行分析比較，根據(jù)網(wǎng)絡的相似程度，確定能夠有效區(qū)分病例組和對照組SNP-SNP相互作用網(wǎng)絡的最佳CMI 閾值。首先，根據(jù)SNPs間的CMI 值，選擇CMI 閾值T 的范圍為0.01～0.58，以0.01 為步長設置58 個閾值。然后，在每個閾值下，對網(wǎng)絡G[ D1] 和G[ D2] 中權值小于閾值的邊進行刪除，權值大于閾值的邊進行保留，分別得到新的58 個病例組和58 個對照組網(wǎng)絡。當T＞0.58 時，病例組和對照組網(wǎng)絡中的孤立點所占比例非常大，邊特別稀疏，這也證實了我們初步確定的閾值范圍是合理的。最后，對比58 個閾值下病例組和對照組網(wǎng)絡的兩個統(tǒng)計量的參數(shù)，并分析統(tǒng)計量參數(shù)能夠區(qū)分兩個網(wǎng)絡的T 的取值范圍，從而確定最佳的CMI閾值。

本文得到的病例組和對照組網(wǎng)絡的統(tǒng)計量隨閾值增加的變化情況如圖1，其中縱坐標表示相應的統(tǒng)計量，橫坐標表示CMI 閾值T，實線表示對照組的情況，虛線表示病例組的情況。

當T＞0.58 時，網(wǎng)絡的平均度和平均介數(shù)趨于0，證實了沒有研究的必要。圖1（a）中，當0.14＜T＜0.3時，病例組與對照組網(wǎng)絡的平均度區(qū)別較大，隨T 的增加，平均度越來越小，這與網(wǎng)絡中孤立點越來越多是對應的。圖1（b）中，當0.17＜T＜0.25時，網(wǎng)絡的平均介數(shù)在兩組中區(qū)別較大，當T 大于一定值時，平均介數(shù)減小，這說明隨著T 增大，網(wǎng)絡中的邊越來越稀疏。根據(jù)兩個網(wǎng)絡統(tǒng)計量確定的T 的范圍，結合圖1，最終選擇最佳CMI 閾值為0.2。在此閾值下，病例組和對照組SNP-SNP 相互作用網(wǎng)絡如圖2所示。

圖1 網(wǎng)絡統(tǒng)計量隨閾值增加的變化情況

圖2 CMI閾值為0.2時，病例組和對照組SNP-SNP相互作用網(wǎng)絡

相同方法可得最佳MI 閾值為0.21，相應的病例組和對照組網(wǎng)絡如圖3所示。

圖3 MI閾值為0.21時，病例組和對照組SNP-SNP相互作用網(wǎng)絡

經(jīng)過多次實驗，得到的最佳CMI閾值均為0.2，證實了本文方法是有效的，也表明最佳閾值為0.2是合理的。圖2 中病例組和對照組網(wǎng)絡有很大差異，對照組網(wǎng)絡節(jié)點間聯(lián)系較弱且存在15 個孤立點；而病例組網(wǎng)絡中，很多孤立點不再獨立且具有了較多聯(lián)系。對比圖2 和圖3，發(fā)現(xiàn)圖3 對照組和病例組網(wǎng)絡中分別有7 個和5 個孤立點，圖2 中一些沒有聯(lián)系的SNPs 節(jié)點，在圖3 較高的MI 閾值0.21 下的網(wǎng)絡中卻存在聯(lián)系，且圖3 中網(wǎng)絡的邊聯(lián)系更加密切，這均證實了MI 具有高估變量間相關性，導致網(wǎng)絡有較高的假陽性邊的問題。

5 致病SNPs的定位

通過節(jié)點中心性可以了解節(jié)點在網(wǎng)絡中的重要性，分為度中心性、接近中心性、介數(shù)中心性以及特征向量中心性。節(jié)點在病例組和對照組網(wǎng)絡中的中心性差異一定程度上決定了節(jié)點的功能。對于SNP節(jié)點i，我們定義Δdi= ||代表該節(jié)點的度中心性差異值，其中和分別代表該節(jié)點在病例組和對照組網(wǎng)絡中的度中心性；同理，定義Δci= ||、Δbi= ||及Δei= ||分別代表該節(jié)點的接近中心性、介數(shù)中心性及特征向量中心性差異值。以上差異值越大，說明節(jié)點越能區(qū)分病例組與對照組網(wǎng)絡。本文將以上差異值細化到了每個SNP節(jié)點，來刻畫其在兩組網(wǎng)絡中的差異。

我們需要設置合適的差異值參數(shù)，如果差異值參數(shù)較小，對致病SNPs的選取限制較少，一些不相關的SNPs 也會被選到可能的致病SNPs 集合內(nèi)，導致假陽性。反之，如果選取過于嚴苛，會遺漏致病SNPs，導致假陰性。圖2 中最佳CMI 閾值0.2 的情況下，兩組網(wǎng)絡的平均度大致相差2.5，因此認為SNP 節(jié)點的度中心性差異值大于2.5 時才有研究的必要，選取Δd ≥3 的SNP，得到包含11 個SNPs的集合S1。此外，兩網(wǎng)絡平均的接近中心性、介數(shù)中心性及特征向量中心性大致相差2e-04、4.5及0.13，同樣可得包含23個Δc ≥2e-04 的SNPs的集合S2，包含20 個Δb ≥4.5 SNPs 的集合S3，及包含16 個Δe ≥0.13 的SNPs 的集合S4。計算S1，S2，S3 及S4的交集，最終得到集合S，包含4 個可能的致病SNPs，如表1。

表1 可能的致病SNPs的信息

本文方法找到的可能的致病SNPs 只有4 個，且其中rs9534318，rs9943876 為預設的致病SNPs。將集合S1、S2、S3 和S4 中的SNPs 按差異值Δd 、Δc 、Δb 及Δe 從大到小排序，rs9534318 的Δd 為4，在S2、S3 和S4 中 分 別 排 第8、第9 及 第1；rs9943876 的Δd 為3，在S2、S3 和S4 中分別排第、第1 及第6，在4 個集合中表現(xiàn)都不錯。我們又在相同的仿真數(shù)據(jù)下，將本文方法找到的可能的致病SNPs 集合與Wang 等［15］利用基于MI 的參數(shù)取值方法選擇出的結構性關鍵SNPs 集合進行了比較，實驗結果顯示兩集合中均包含了預設的致病SNPs，但本文集合普遍較小，避免了互信息存在有偏估計，導致錯誤率偏高的問題。

根據(jù)上述定位致病SNPs 的方法，我們分別針對有1個、2個和3個預設致病SNPs的情況，對病例組和對照組含500、1000、3000、5000 個個體的情況進行了多次實驗，基本每次實驗得到的可能的致病SNPs集合S中都包含預設的致病SNPs，且S大小合適，這說明本文定位致病SNPs 的方法比較準確有效。

6 結語

本文采用HapMap3 計劃13 號染色體上乳腺癌相關的BRCA2 基因數(shù)據(jù)，利用Hapgen2 仿真病例組和對照組SNP 數(shù)據(jù)，基于CMI 計算SNPs 間的相關性，構建了病例組和對照組SNP-SNP 相互作用網(wǎng)絡，并在最佳CMI 閾值下，根據(jù)節(jié)點的網(wǎng)絡中心性差異值參數(shù)，篩選出了可能的致病SNPs。實驗結果表明本文能夠高效準確地選擇出預設的致病SNPs。但是綜合考慮多個乳腺癌易感基因數(shù)據(jù)，進而定位可能的致病SNPs還需要進一步的研究。