引言

大量研究表明，人体的许多表型性状差异以及对药物和疾病的易感性都可能与某些基因位点相关联，故选取恰当的方法来找寻出相关致病位点，对于疾病的治疗和预防具有重要意义。

目前关于遗传位点与患病信息关联性的研究成果较多，其中全基因组关联分析是研究患病信息和患病位点之间关联性的重要方法，这种方法通过对试验个体的健康状况和位点编码的统计关联分析来确定致病位点，从而发现遗传病或性状的遗传机理^[1]。皮尔逊卡方检验^[2]和Logistic回归模型^[3-4]也是致病位点关联性分析中常用的方法，目的是剔除与患病信息无关或影响作用非常小的基因位点。为了避免模型设定的错误，部分学者提出利用神经网络模型等非参数的方法拟合数据并进行位点筛查^[5-6]，但神经网络筛选方法是根据对每个位点在神经网络中的权重贡献率大小来挑选的，并没有考虑到位点之间是否存在交互作用，因此有时不能准确反映位点与遗传疾病的关联强度。

实际上，如癌症、糖尿病等复杂疾病并不是由单个位点独立引起的，而是与多个遗传位点的联合作用有关，即位点之间存在交互作用。交互作用在疾病的认识和发展中扮演着极其重要的角色，如果没有考虑位点间的交互作用，就无法真实准确地描述位点与患病的效应^[7-8]。探讨位点间交互作用对于提高复杂疾病的遗传解释度、构建复杂疾病的遗传风险评估模型、开发疾病诊疗个性化药物靶点并最终降低复杂疾病负担等方面均具有重要的理论和现实意义^[9-10]。文献[11-12]通过Logistic回归模型研究了位点及位点交互作用对遗传疾病的影响，得到对遗传疾病影响较大的位点组合。但是Logistic回归模型对位点及位点交互作用作了较强的参数形式的假定，这在处理实际问题时可能会产生模型设定误差，从而导致错误的推断结果。

因此，针对位点之间可能存在交互作用的情形，本文提出新的方法筛选与遗传疾病关联性最强的位点组合。针对每个可能的候选位点组合，利用神经网络方法训练数据，并以预报准确率作为评价标准，利用粒子群算法(PSO)通过迭代计算逐步逼近与遗传疾病相关性最强的位点组合，该方法缩短了计算时间，避免了Logistic回归模型等方法对位点作用较强的参数形式的假定。

1 基于粒子群算法的位点组合选取 1.1 基于粒子群算法的迭代计算

本文将预报准确率作为评价标准，将全部给定的候选位点集合看作样本空间，把样本空间的所有非空子集作为候选的位点组合，根据给定的评价标准找出关联性最强的位点组合。其中预报准确率根据以下方法得到：首先将所有样本数据分为训练样本和预测样本，然后以相应位点的取值为输入、遗传疾病信息为输出，分别对某个指定的候选位点组合利用神经网络方法训练数据，最后基于预测样本进行预报并计算出相应的预测精度。

但实际应用中由于候选的全部基因位点数目巨大，导致寻找最优位点组合的计算强度大时间过长，因此本文利用粒子群算法，通过迭代逼近最终筛查出关联性最强的位点组合。粒子群中的每一个粒子均对应一个候选的位点组合，其包含的位点在全部候选位点中的位置与粒子中分量为1的位置相对应。以每个粒子对应的预测精度为适应度函数值，根据粒子群算法对粒子进行更新迭代直至收敛，从而得到给定评价标准下的最优粒子，即最优位点组合。

粒子群算法的原理为：在迭代搜寻的过程中，粒子通过向个体最优位置和群体最优位置学习的经验来搜索最优解，在搜寻过程中通过不断修正粒子的适应度、速度和位置来学习。正是由于每个粒子都在不断调整和其他较好位置粒子的差异，不断向前期找到的较好位置靠拢，才使得它们的收敛速度很快^[13-14]。吕思晨^[15]提出了将粒子群算法和遗传算法结合起来研究位点与疾病关联分析的算法，并用来寻找与疾病关联性较强的单个位点。

假定粒子群中有n个粒子，每个粒子均是分量为0或1的D维向量(D为可选的位点个数)，即S=(S₁, S₂, …, S_n)为粒子种群，其中S_i为第i个粒子。为了计算每个粒子的适应度，分别以粒子中所有取值为1的分量所对应的位点取值作为BP神经网络的输入值，则对预测样本的预报正确率可作为对应粒子的适应度函数值。在下一次迭代时，基于每个粒子的适应度值以及种群中最大的适应度值来调整其分量的飞行速度，然后基于速度修正粒子的选取概率，并重新生成一个新的具有n个粒子的种群。重复以上迭代步骤进行计算，直至收敛。

具体来说，粒子的速度决定粒子移动的方向和距离，在每次迭代过程中，第i个粒子的第j个分量根据个体适应度极值和群体适应度极值更新自身的速度，即：

$ v_{ij}^{k + 1} = \omega v_{ij}^k + {c_1}{r_1}\left( {P_{ij}^k-s_{ij}^k} \right) + {c_2}{r_2}\left( {P_{gj}^k-s_{ij}^k} \right) $

(1)

其中，ω为惯性权重；k为当前迭代次数；c₁和c₂是非负常数，称为加速度因子；r₁和r₂是[0, 1]上的随机数；P_i=(P_i1, P_i2, …, P_iD)^T为第i个个体适应度最优时对应的取值，表示迭代到第k次时第i个粒子的适应度最大取值(比较的是同一个粒子在整个迭代中的适应度值)；P_g=(P_g1, P_g2, …, P_gD)^T为全局适应度最大取值，表示迭代到第k次时所有粒子的适应度最大值(比较的是所有粒子在整个迭代中的适应度值)。

根据第i个粒子中第j个分量的速度v_ij更新对应位置的选取概率f(v_ij)，基于选取概率来判定第i个粒子对应的第j个分量的取值S_ij

$ {S_{ij}} = \left\{ \begin{array}{l} 1\;\;\;\;\;\;\;R < f\left( {{v_{ij}}} \right)\\ 0\;\;\;\;\;\;\;{\rm{others}} \end{array} \right. $

(2)

其中R是服从[0, 1]上均匀分布的随机数，函数$f\left( {{v_{ij}}} \right) = \frac{{\exp \left( {{v_{ij}}} \right)}}{{1 + \exp \left( {{v_{ij}}} \right)}} $，v_ij不同，粒子的更新位置也不同，v_ij越大，粒子在该位置取1的概率越大。为了能够迭代计算，在初次迭代时，每个粒子的每个分量的取值可利用等概率来选择。模拟结果显示，算法的结果不受初值的影响。

1.2 基于神经网络的适应度计算

假定输入某l个位点构成的位点组合为

$ {\mathit{\boldsymbol{X}}_p} = {\left( {x_1^p, x_2^p, \cdots, x_l^p} \right)^{\rm{T}}} $

其中p=1, 2, …, W(W为训练样本总数)。神经网络第p个样本隐层第j个结点的输出值为

$ h_j^p = f\left( {{a_{0j}} + \sum\limits_{i = 1}^l {{\omega _{ij}}x_i^p} } \right), j = 1, 2, \cdots, L $

式中，f(·)为激活函数，其形式同式(2)中f(v_ij)，ω_ij为输入层和隐层之间的连接权值，a_0j为常数项，j为隐层结点数。

BP神经网络第p个样本输出层第k个结点值为

$ y_k^p = \sum\limits_{j = 1}^L {{\omega _{jk}}h_j^p, } k = 1, 2, \cdots, N $

(3)

其中，ω_jk为隐层到输出层的权值，h_j^p为第j个隐层结点的输出值，N为输出层结点数。

得到训练的网络结构后，输入预测样本并比较预测结果和真实结果，以位点组合对预测样本的预测准确率作为对应的粒子在粒子群算法中的适应度函数值。

1.3 算法构建

位点组合选取迭代计算步骤如下。

(1) 利用两点分布B(1, 0.5)生成n个粒子的每个分量，建立BP神经网络；对每个粒子，以粒子中所有位置为1的位点作为输入值，以预测样本的正确率作为适应度函数评估各粒子，记录第i个粒子的个体最优值并作为当前粒子适应度，全局最优值为适应度值最大的粒子。

(2) 根据PSO算法的公式(1)和(2)更新n个粒子的速度和位置，产生新的一组粒子，建立BP神经网络；将每个粒子以粒子中所有位置为1的位点作为输入值，以预测样本的正确率为适应度函数来评估各粒子。比较当前第i个粒子和个体最优的适应度函数值，将其中具有较大适应度函数值的粒子作为第i个粒子的个体最优；比较所有粒子和全局最优的适应度函数值，将具有最大适应度函数值的粒子作为全局最优。

(3) 判断是否满足停止准则，若满足，则将全局最优输出，结束；若不满足，则返回步骤(2)。

(4) 重复步骤(1)~(3)，可以得到多个位点组合。将出现在位点组合中次数最多的部分位点作为具有较好预测效果的位点组合。

2 实验及结果分析 2.1 数据来源

选取2016年研究生数学建模竞赛B题数据(http://gmcm.seu.edu.cn/01/1d/c12a285/page.htm)，样本由1000位试验者的患病信息和9445个位点信息遗传信息构成。样本分为患病者和健康者，两者各占50%，用1表示患病者、0表示健康者。每位试验者对应9445个位点，位点信息由碱基A，T，C，G的不同组合来表示，用两个碱基的组合表示一个位点的信息，一个位点有3种不同编码。

因为样本所给的位点有9445个，直接实施本文算法会导致计算量非常大，计算时间过长，所以需要首先对与患病信息无关的基因位点初步筛除。本文选择皮尔逊卡方检验的p值和Logistic回归中单变量的t检验的p值作为统计相关性度量指标。筛选掉与遗传疾病信息完全独立的基因位点和本身作用很小的基因位点，然后从剩余的相关性较强的位点集合中分别利用本文方法和神经网络权重分析法筛查关联性最强的部分位点，并根据预报准确率对两种方法进行对比。

2.2 初步筛选步骤 2.2.1 皮尔逊卡方检验

皮尔逊卡方统计量^[2]是用于检验实际分布与理论分布拟合优度指标，可以用于两个指标的独立性检验。以位点rs2273298位置为例，该位点取值为AA、AG和GG，并分别赋值为1, 2, 3；疾病信息取值为0、1。令t_ij为位点取值j而疾病取值为i的个体的数量，则位点rs2273298与疾病信息的独立性检验的卡方值为

$ {\chi ^2} = \sum\limits_{i, j} {\frac{{{{\left( {{t_{ij}}-{f_i}{q_j}/1000} \right)}^2}}}{{{f_i}{q_j}/1000}}} $

(4)

其中，${f_i} = \sum\limits_{j = 1}^3 {{t_{ij}}, i = 0, 1} $；${q_j} = \sum\limits_{i = 1}^2 {{t_{ij}}}, j = 1, 2, 3 $。

当位点与疾病独立时，式(4)中统计量服从自由度为2的卡方分布，因此可计算出样本卡方值对应的p值。若p值很小，说明位点与疾病之间存在相关性；若p值大于给定的临界值，说明两者之间相互独立。

利用MATLAB计算9445个位点与疾病之间的卡方统计量值和相应的p值，以卡方检验的p值小于0.01为标准，从中筛选出与疾病存在相关性的73个位点，结果如表 1所示。

2.2.2 Logistic回归

本文中位点的取值为碱基对，每一个位点编码方式均为3种。考虑到每个位点有3个取值，将位点拆分为两个0-1型变量。以第i个位点rs2273298为例，样本中位点rs2273298有AA、AG和GG3种不同编码方式，将此位点拆分为两个示性变量

$ {x_{2i-1}} = \left\{ \begin{array}{l} 1\;\;\;AA\\ 0\;\;\;{\rm{other}} \end{array} \right. $

$ {x_{2i}} = \left\{ \begin{array}{l} 1\;\;\;\;GG\\ 0\;\;\;\;{\rm{other}} \end{array} \right. $

根据此方法将通过卡方检验得到的73个位点拆分成146个示性变量。

观测指标变量为(X, y)，其中X=(x₁, x₂, …, x_2m)是m个(m=73)位点的信息，y∈{0, 1}是一个二分类的属性变量(y=1表示第i个样本是患病者，y=0表示第i个样本是健康者)。用Logistic回归^[3]建立患病识别准则模型：

$ p = \frac{{\exp \left( {{\beta _0} + {\beta _1}{x_1} + {\beta _2}{x_2} + \cdots + {\beta _{2m}}{x_{2m}}} \right)}}{{1 + \exp \left( {{\beta _0} + {\beta _1}{x_1} + {\beta _2}{x_2} + \cdots + {\beta _{2m}}{x_{2m}}} \right)}} $

(5)

其中，β_i(i=1, 2, …, 2m)为变量系数，Logistic回归值p为个体样本是否患病的概率。

根据回归系数的t检验来判断位点与疾病之间的相关性是否显著。由于样本量较大，t检验统计量近似服从标准正态分布

$ T = \frac{{{{\hat \beta }_i}-{\beta _i}}}{{{\sigma _{{{\hat \beta }_i}}}}} \sim N\left( {0, 1} \right) $

(6)

由式(6)结合标准正态分布表可计算出t检验统计量的p值

$ p = P\left\{ {\left| T \right| > t} \right\} $

以146个示性变量为自变量进行Logistic回归，以t检验的p值来衡量位点是否对患病概率存在影响。以p值小于0.05作为筛选标准进一步筛选出变量所对应的位点，得到显著相关水平较高的55个位点如表 2所示。

2.3 最优位点组合求解

根据初步筛选得到55个与某遗传疾病存在关联性的位点及位点组合的预报准确率，可以找出对遗传疾病发生可能性影响最大的位点组合。但是可能的位点组合数目有2⁵⁵-1个，要逐一计算对比寻找最优位点组合将产生非常大的计算量。为此本文利用粒子群算法通过迭代逼近来求出最优位点组合，在采用神经网络方法针对候选位点训练拟合数据时将样本分为训练样本和预测样本，其中训练样本900个，预测样本100个。将BP神经网络对预测样本的预测正确率作为粒子群算法中的适应度函数值，然后进行迭代求解。

2.4 结果对比及分析

根据在多个最优解中出现的次数选出重要程度最高，即与疾病关联性最强的11个位点如表 3所示。

将上述11个位点作为输入序列，患病信息为输出量，从真实数据和预测数据的比较(图 1)可以看出模型具有较好的预测效果。

为了比较本文所提出的基于神经网络和粒子群算法的位点组合的筛选方法(方法1)与单一神经网络权重分析法(方法2)的准确度差异，以55个初选得到的位点取值为输入值直接训练数据，以每个位点在神经网络中的权重贡献率大小排序，再从中选取权重贡献率最大的11个位点作为与疾病关联性最强的位点。

利用混淆矩阵和受试者工作特征曲线ROC下方面积(AUC)检验这两个方法预测效果，具体结果如表 4所示。

由表 4结果可知，基于粒子群算法所筛选的位点组合预测效果比单一神经网络权重法更好，尤其是针对患病者的预测精度高达86%，说明本文方法分类正确率更高。

3 结束语

针对寻找与疾病存在关联性的位点问题，本文提出以神经网络的预报准确率为评价标准，基于粒子群算法筛选与某遗传疾病相关的重要位点组合。建立的位点筛选方法具有较强的适应性，包含了位点之间存在交互作用的情形，与只考虑选取单个重要位点的单一神经网络权重法相比，由本文方法得到的位点组合的预测效果更好，适用范围也更广泛，为疾病诊断和遗传疾病分析提供了一种新方法。