随着人类基因组测序工作的逐步完成, 大量的数据为全基因组关联分析提供了丰富的素材, 也涌现出许多数据分析方法^[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. 本文通过互信息建立相互作用网络, 从而分析网络结构的差异性. 设 $X = ({x_1},{x_2},...,{x_n})$ , $Y = ({y_1},{y_2},...,{y_m})$ 是两个SNP的基因型数据在个体之间表达形成的向量, $ p({x_i},{y_j})(i = 1,2,...,n;$ $j = 1,2,...,m)$ 是X和Y的联合概率分布, $H(X,Y)$ 是他们之间的联合熵, 定义为:

$H(X,Y) = \sum\nolimits_{i = 1}^n {\sum\nolimits_{j = 1}^m {p({x_i},{y_j})} } {\log _2}\left[ {\frac{1}{{p({x_i},{y_j})}}} \right]$

(1)

对于两个随机变量之间存在的关系, H(X)表示随机变量X蕴含的不确定性, 而条件熵 $H(X,Y)$ 则是已知条件Y时随机变量X所余下的不确定性, 那样, $H(X) - H(X|Y)$ 就表示已知条件Y后X包含的信息量. 进而还可以证明这个值关于X和Y是对称的, 即 $H(X) - H(X|Y) = H(Y) - H(Y|X)$ , 且都等于 $H(X) + H(Y) - $ $ H(X,Y)$ . 由此X和Y之间的互信息可以计算, 互信息记为 $MI(X,Y)$ :

$MI(X,Y) = H(X) + H(Y) - H(X,Y)$

(2)

$MI(X,Y)$ 越大, 表明X和Y的相关程度越大; $MI(X,Y) = 0$ 表示X和Y相互独立. 因此, 互信息表达了两个SNP之间相互依赖的程度.

因为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之间的互信息构成互信息矩阵, 记作 $M = {\{ {W_{ij}}\} _{m*n}}(|I| = n)$ , 矩阵中的每行代表一个SNP, 每一列代表此SNP与另一个SNP之间的互信息.

假定存在一个集合的SNP基因型数据D, 其中所拥有的SNP的集合我们记作I, $M = {\{ {W_{ij}}\} _{m*n}}(|I| = n)$ 可由互信息计算公式(2)得到一个互信息矩阵. 定义一个建立在关于SNP基因型数据D的互信息网络. $G[D] = (V,E;w)$ 是边赋权图, 其中V表示点集合、每个网络中的节点i∈V表示一个SNP, ${\forall _{i,j}} \in V$ , 基因i和j之间的互信息计算值w_ij定义为每条边 $(i,j) \in E$ 的权重. 在下面的表述中, 我们将基因i∈I以及顶点i∈V等同起来看待.

2 基于网络统计量的关键基因选取

利用上述方法得到的SNP相关网络中各节点(SNP)的网络结构参数来描述特定生物过程中基因的重要性. 首先给出几个重要的能够反映网络结构特点的网络统计量的相关定义^[11].

(1) 度(K): 在网络中, 度指的是与该点相连接的边数目. 节点度可以表示该点的重要程度, 节点度越大, 表示该点在网络中越重要. 而网络的平均度可以通过计算所有的点的度, 后取平均数计算得到.

(2) 平均路径长度(L): 定义为网络中所有的点之间两两求得的距离的平均数, 网络中的任意两点i, j的距离即边的条数, 则两点之间的平均路径长度表示为所有的点之间的平均距离, 记作: $L = \displaystyle\frac{2}{{N(N - 1)}}\sum\nolimits_{i > j} {{d_{ij}}} $ , 其中N表示网络中的节点数目.

(3) 聚类系数(C): 网络中节点i有K_i个边与之连接, 那么与该点可能连接的最大边数为 ${K_i}({K_i} - 1)/2$ , 若这K_i个节点之间真实边为E_i, 则它与总的所有情况下的边比例, 计算得到节点i的聚类系数C_i: ${C_i} = 2{E_i}/$ ${K_i}({K_i} - 1)$ . 很显然, 0≤C≤1. C=0代表网络中的点为孤立点; C=1表示网络中的所有点之间都是互相连接的, 视为全局耦合网络.

(4) 介数(B): 网络中介数的概念可以分为两类, 一类是点介数, 另一类是边介数. 节点k的介数定义为 ${B_k} = \sum\nolimits_{i \ne j} {{C_k}} (i,j)/C(i,j)$ , 其中, C(i, j)代表i与j间最短路径总数, ${C_k}(i,j)$ 表示中间点为k时, i与j间的所有路径总数. 介数反映了节点k在i和j之间的流通量和重要程度. 网络中某个节点的介数越大, 说明该点在网络中信息传播的信息量就越大, 越容易在该点造成网络堵塞. 假设两组连接度很高的网络中间只有少数点连接, 那么这几个少数点介数就会很大, 即很多的信息在流通的过程中经过这几个点, 很容易造成堵塞, 从而造成数据信息丢失. 因此, 最大介数的增大会降低网络同步能力.

(5) 模块度(Q): 模块度也称作模块化度量值, 是用来衡量网络强度的统计量. 最早是Newman提出的, 它用来描述网络社团以及划分的好坏. 假定网络共分为k个社团, $E = ({e_{ij}})$ 代表一个k×k维的矩阵. 故模块度可以定义为: $Q = \sum\nolimits_i {{Q_i}} = \sum\nolimits_i {({e_{ij}} - {a_i}^2)} $ , 其中, ${a_i} = \sum\nolimits_i {{e_{ij}}} $ 是矩阵中的数值之和(行或列), e_ij用来表示社区i和社区j之间的边的数量. 模块度可以区分社区划分的好坏. 若是划分的好, 则社区内部节点相似度较大, 而在社区外边相似度较低. Q越大, 越接近1, 代表社区拥有一个很好的划分结构, 使得社区的划分合理化. 通常设定的值是在0.3与0.7之间.

本文中我们主要选择5个参数进行分析比较, 对于给定的参数进行最终的分析, 从而找到影响网络的重要因素, 依据此类统计量进行归纳分析, 得出相应的参数.

我们对由SNP数据设定不同的互信息阈值而形成网络, 针对其中大于阈值的边, 做去掉处理, 而针对小于阈值的边进行保留操作, 从网络图可以分析出统计量所对应的参数变化, 得到有益信息量.

根据网络中SNP之间互信息计算的值, 选择阈值范围为0.1到0.63. 共设置63个阈值, 在每个阈值的条件下, 统计计算相应的网络结果, 从而得到一致性网络, 根据网络的相似性程度选择对实验组和对照组差别较大的统计量进行分析. 我们最终选择了度作为区分依据, 并分析能够区分实验组和对照组的取值范围, 得出最佳的阈值, 对于不同的数据, 得到的互信息值也不同, 所以需要根据数据得到的互信息范围, 由网络统计量得到取值范围, 得到互信息取值的交集, 能够区分对照组和实验组数据, 从而确定最佳的互信息阈值. 这样就能够保证所取的阈值不受样本数量的大小影响, 而是根据样本的不同情况得到相应的阈值. 对于节点i, 我们定义, $\Delta d = {d_i}^d - {d_i}^c$ , ∆d代表了这个节点的度差异值, 在该公式中, ${d_i}^d$ 和 ${d_i}^c$ 分别代表了这个节点在实验组与对照组网络中节点的度.

我们都知道, 在复杂网络中, 节点度能够代表节点的作用和影响力. 本文从网络结构差异的角度去衡量各个统计量^[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中, 纵坐标表示相应的统计量, 横坐标代表阈值, 虚线表示对照组数据显示效果, 实线表示实验组数据显示效果.

实验发现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可以看出, 每一个统计量都有不同的阈值范围, 平均度K的范围较大, 0.08<K<0.35; 其他的统计量阈值范围相差不大, 基本在0.2到0.3之间. 结合图1, 选择0.28为阈值构建网络.

依据图2, 实验组和对照组的图像是有很大差异的. 在对照组, 节点之间联系较弱且存在更多的孤立点. 但是在实验组中, 很多的孤立点不再是独立的, 并且拥有了更多的联系. 对照组中存在36个连接点和9个孤立点, 而实验组中存在39个连接点喝6个孤立点. 这表明我们选取的阈值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.

4.3 参数评估

在查找“结构性关键SNPs”时, 我们需要从网络平均度出发, 对选取网络中的关键SNPs设置合适的差值参数. 如果选取的差值参数比较小, 对SNPs选取限制比较宽泛, 一些不相关的SNPs也会选取到SNPs集合内, 从而导致假阳性. 反之, 如果选取过于严苛, 反而会遗漏一些比较重要的节点, 导致假阴性.

我们选取基因BRCA2, 得到它在阈值为0.28时候的网络, 如图2所示. 选择不同的差值参数, 得到一系列不同的结构性关键SNPs, 如表3所示.

当互信息阈值设定为0.28时, 网络中度的最大变化量是5. 当r≥5时, 所得的关键SNPs只有节点39, 对网络影响较大的节点25却被忽略. 当r≤2时, 所得的关键SNPs只有13个, 这里面也包括了其中的非零点.

5 结论与分析

本文通过国际项目HapMap3中以及Hapgen2软件生成的13号染色体上BRCA2基因生成仿真数据. 利用互信息表示SNPs间的相互作用. 构建实验组和对照组的网络, 根据阈值及差值参数筛选出关键SNPs. 最后, 对我们所选择的参数进行了评估, 证明我们所选定的参数能够反映结构的变化, 能够较好地选择出预设的关键SNPs. 通过数值实验发现: 样本数目会影响互信息的大小, 样本数较小时, 互信息较高, 样本数较大时, 互信息逐渐降低, 本文认为, 样本数偏少, 则特异性个体数目不完备, 样本数过多, 又会造成冗余, 增加了计算复杂度. 目前, 确定合适的上下界仍然是一个具有挑战的问题.