基于条件互信息的乳腺癌易感基因网络的构建和分析∗

王淑栋 范晓丹 王新赠

（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还需要进一步的研究。