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

王淑棟 范曉丹 王新贈

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

1 引言

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

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

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

2 數據來源與處理

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

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

3 基于CMI 的SNP-SNP 相互作用網絡的構建

假設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 時的聯合概率值。MI是根據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 條件下的聯合概率，p( x，y，z )表示變量X 、Y 和Z 的聯合概率。CMI是根據變量X 和Y 在變量Z 下的條件獨立性評估的，定義如下：

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

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

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

4 最佳CMI閾值的選取

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

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

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

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

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

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

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

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

5 致病SNPs的定位

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

我們需要設置合適的差異值參數，如果差異值參數較小，對致病SNPs的選取限制較少，一些不相關的SNPs 也會被選到可能的致病SNPs 集合內，導致假陽性。反之，如果選取過于嚴苛，會遺漏致病SNPs，導致假陰性。圖2 中最佳CMI 閾值0.2 的情況下，兩組網絡的平均度大致相差2.5，因此認為SNP 節點的度中心性差異值大于2.5 時才有研究的必要，選取Δd ≥3 的SNP，得到包含11 個SNPs的集合S1。此外，兩網絡平均的接近中心性、介數中心性及特征向量中心性大致相差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 個集合中表現都不錯。我們又在相同的仿真數據下，將本文方法找到的可能的致病SNPs 集合與Wang 等［15］利用基于MI 的參數取值方法選擇出的結構性關鍵SNPs 集合進行了比較，實驗結果顯示兩集合中均包含了預設的致病SNPs，但本文集合普遍較小，避免了互信息存在有偏估計，導致錯誤率偏高的問題。

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

6 結語

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