基于深度学习和进化计算的特征基因选择方法

摘要

一种基于深度学习和进化计算的特征基因查找方法，包括以下步骤：1)计算差异表达水平指数，根据二八准则建立一级基因池；2)根据每个基因的表达水平，计算其映射后的密度矩阵和距离矩阵，绘制决策图，利用多重线性回归分析拟合二元平面，自动确定聚类中心；3)构建深度基因表达预测网络计算一级基因池的基因‑基因敏感性信息GGSI，4)依据GGSI值剔除冗余基因，建立二级基因池；5)基于GGSI值对布谷鸟搜索算法进行二进制编码，选出最紧凑的基因集，建立三级基因池。本发明提供一种基于分级结构建立的特征基因选择框架，能够较好的提取出关键基因，采用多重线性回归分析结合深度学习算法和优化算法选择最紧凑特征基因集的方法。

说明书

技术领域

本发明属于生物信息领域，具体涉及一种特征基因的选择方法。

背景技术

随着基因测序技术的进步，高通量测序技术提供了大量基因表达数据。基因表达谱作为疾病发生、遗传扰动和药物治疗后的细胞表达模式的捕获工具而被广泛应用。面对大量高维的基因表达数据，如何从中挖掘有用的信息成为生物信息领域的研究重点。

由于基因表达数据的高维性质，挖掘其中具有关键信息的特征基因对后续的基因数据处理和样本表现型分析具有十分重要的意义。目前，有研究人员结合二进制粒子群优化算法(BPSO)和过滤法寻找到最佳的基因子集，从而简化了基因选择的过程并获得了相对遗传算法(GA)而言具有更高的分类精度。更进一步的研究，通过结合遗传算法、粒子群优化算法和极限学习机(ELM)设计了一种新的基因选择方法。然而，这些混合方法在起到简化基因集的作用和达到处理高效性的同时，缺乏可解释性。

因此，针对基因表达间的高度非线性相关性和庞大数据量，结合深度学习算法和二进制布谷鸟搜索(BCS)算法，建立神经网络和分层模型，简化特征基因选择与样本分类过程，增强其可解释性，对生物信息领域的基因表达分析具有重大意义，能够指导后续的表现型关联性分析、癌症靶向基因标志物检测等临床应用。

发明内容

为了克服已有特征基因的选择方法的不足，本发明提供一种能够实现聚类中心自动确定、样本分类更加精确的基于深度学习和进化计算的特征基因选择方法。

本发明解决其技术问题所采用的技术方案是：

一种基于深度学习和进化计算的特征基因选择方法，所述方法包括以下步骤：

1)选择差异表达的基因，建立一级基因池，过程如下：

1.1)计算原始基因池中各个基因的差异表达水平指数，即IIC-FC指数：

式(1)中，c表示原始基因池中的基因个数，和分别表示基因i和基因j的表达水平均值，和分别表示基因i和基因j的表达水平标准差，函数max{·,·}和min{·,·}分别表示取最大值和最小值，ln(·)为对数函数；若IIC-FC指数越大，表明该基因包含的样本分类信息丰度越高，用于样本分类获得高精度的可能性越大；

1.2)根据二八准则，按照步骤1.1)中得到的IIC-FC值从高到低对原始基因池中的所有基因进行排序，然后从中选择前20％的基因作为一级基因池；

2)对一级基因池中的基因进行自动聚类，过程如下：

2.1)基于密度聚类算法计算一级基因池中每个基因的距离值δ和密度值ρ；定义欧拉距离作为基因i和基因j之间的相似度距离，计算公式如下：

式(2)中，N_sam表示样本个数，gE(i,k)表示基因i在样本k中的表达值；

首先计算各个基因的密度值ρ，得到密度矩阵计算公式如下：

其中

式(3)中，ρ_i表示基因gene_i的密度值，一级基因池为相应的指标集为I_P＝{1,2,...,N_pool1}；其中定义表达值累加和最大处的基因gene_max与最小处的gene_min之间的相似度距离的2％作为d_c值，计算公式如下：

d_c＝0.02*d_gene(gene_max,gene_min)>

然后计算各个基因的距离值，得到距离矩阵每个基因gene_i的距离值定义为δ_i，首先查找比gene_i密度大的基因集，得到集合P'＝{gene_j}，然后查找P'中与gene_i的距离最近的基因，则得到δ_i＝d_gene(i,j')；

定义函数f_γ是关于变量ρ和δ的二元离散函数，对应于三维空间中的坐标值是(ρ,δ,f_γ)，则得到双变量离散函数为：

式(5)中，f_γ取密度值和距离值的乘积的对数值作为函数值；表示大于零的较小正数(一般取)，其目的是为了在密度为零时，即没有点落在以d_c为半径的圈内时式子仍能够成立，并没有实际物理意义；

2.2)根据步骤2.1)中得到的双变量离散函数γ＝f_γ(ρ,δ)，利用最小二乘法进行二元斜面的拟合，得到拟合平面为z_γ＝b₁+b₂ρ+b₃δ，计算各个数据点的残差值为ε_γi＝y_γi(ρ,δ)-γ_i(ρ,δ)，绘制残差直方图ε_γi-h，并根据钟型曲线的正态拟合得到方差值为σ_γ，利用3σ原则自动确定处在置信区间外的η个奇异点作为聚类中心，记为c_γ；

3)构建深度基因表达预测网络(DGEPN)，计算一级基因池的“基因-基因”敏感性信息(GGSI)，过程如下：

3.1)区分样本的训练集、验证集和测试集，搭建多层神经网络作为深度基因表达预测网络；首先假定有N_sam个样本，每个样本在原始基因池中有N_gene维基因表达，其中M_in个基因作为神经网络的输入层，即将一级基因池中的所有基因作为输入；M_out个基因作为神经网络的输出层，即将原始基因池中对应一级基因池的补集作为输出；满足M_in＝N_pool1，M_in+M_out＝N_gene，则每个样本可以表示为然后随机挑选N_train个样本组成训练集，N_val个样本组成验证集，N_test个样本作为测试集，其中N_train:N_val:N_test＝4:1:1；

在神经网络的输入层和输出层间具有l_hidden层隐藏层(一般取l_hidden＝3)，隐藏层的神经元之间进行全连接；第i层隐藏层的神经元个数NumNeu_i为：

NumNeu_i＝M_in+i*Z{abs(M_in-M_out)/l_hidden+1}>

式(6)中，函数abs(·)表示取绝对值，Z{·}表示取整函数；为了减少训练参数从而加快训练，对神经网络中的隐藏层设置dropout比率(一般取dropout＝0.25)；为了克服梯度弥散问题，激活函数采用ReLU函数；

3.2)计算基因间的敏感性信息，基于步骤3.1)中搭建的深度基因表达预测网络，定义第i个输入基因到第j个输出基因的敏感性信息为：

式(7)中，表示第j个输出神经元对第i个输入神经元的导数；敏感性信息指数越大，说明第i个基因对第j个基因的非线性相关性越高；

3.3)根据步骤3.2)中得到的基因间敏感性信息，计算神经网络输入层中每个输入基因归一化后的敏感性指数：

3.4)训练神经网络，定义输出值与期望值的平均平方误差作为损失函数，即：

当训练代数达到预设的最大值或者损失值小于预设的阈值，结束训练，并用测试集的平均绝对误差(MAE)作为神经网络模型的评价指标：

误差越小，说明计算GGSI的网络模型性能越好；

4)建立二级基因池，过程如下：

4.1)计算步骤2)中得到的每个基因簇的最大GGSI值和平均其中K是基因簇的个数；然后依据簇间剔除和簇内剔除策略构建二级基因池；

4.2)首先进行簇间剔除，若第j个簇的GGSI最大值远小于其它簇的GGSI最大值，即则剔除该基因簇中的所有基因；

4.3)然后进行簇内剔除，若第i个簇内的第j个基因的GGSI值小于该簇的GGSI平均值，即则剔除该基因，其中NumK_i表示基因簇i中包含的基因个数；

5)基于二进制布谷鸟搜索算法进行优化搜索最紧凑的基因集合，建立三级基因池，过程如下：

5.1)进行鸟巢的编码初始化，首先将二级基因池中的每个基因根据其GGSI值从高到低进行排序，然后对第i个鸟巢的第j维，即第i个解的第j个基因进行初始化编码：

式(11)中，d表示解的维度，即二级基因池中的基因个数；式(11)表示GGSI排名靠前20％的基因和排名靠后20％的基因分别编码为1和0，剩余基因进行0或者1的随机编码；

5.2)进行糟糕解的替换，为了保证新产生的解都是由优秀的个体组成，基于概率p_a(一般取p_a∈[0,1])淘汰不适应环境的糟糕解，并将它们替换为新的解决方案，替换公式为：

式(12)中，d表示解的维度；式(12)表示对某个糟糕解，将其GGSI排名靠前10％的基因和排名靠后10％的基因分别编码为1和0，剩余基因的编码不变；

5.3)进行鸟巢位置的更新，更新公式为：

式(13)中，Sig(·)表示sigmoid函数，Avg(·)表示二级基因池中所有基因的GGSI的平均值，其中

Levy(λ)～u＝s^-λ，λ∈(1,3]，α＞0表示步长尺度缩放因子；

5.4)布谷鸟搜索优化算法的适应度函数定义为：

fitness(net_i)＝MAE(net_i)+ξ*GeneNumber(net_i)>

式(14)中，MAE(net_i)表示第i个解在二级基因池中所选择的基因集合在DGEPN中预测剩余基因表达的平均绝对误差，GeneNumber(net_i)表示第i个解所选的基因个数，ξ∈{10^-2,10^-3,10^-4}是一个尺度规模参数，ξ的取值根据GeneNumber(net_i)的值确定；

5.5)进行优化搜索迭代，当达到预设的迭代次数，或者适应度函数值小于预设的阈值，则停止迭代，将此时的最优解进行解码，得到最优三级基因池。

进一步，在所述步骤1.1)中，通过计算差异表达水平指数，体现基因包含的样本分类信息丰度，并且适用于多分类数据集的基因差异表达水平计算。

再进一步，在所述步骤2)中，通过计算每个基因的密度值和距离值，确定二元离散函数，并对其进行多重线性回归分析得到二元拟合平面，依据拟合平面和预设的置信区间自动确定聚类中心。

更进一步，在所述步骤3)中，通过构建深度神经网络学习基因间的敏感性信息，并将其作为步骤4)中一级基因池的冗余基因剔除依据和步骤5)中布谷鸟搜索算法的二进制编码依据。

在所述步骤1)、步骤4)和步骤5)中，通过构建三级基因池查找最紧凑的特征基因，能够尽可能的剔除冗余基因，同时尽可能的保留具有高分类信息丰度的特征基因。

本发明的技术构思为：基于深度基因表达预测网络和二进制布谷鸟搜索算法，通过分级筛选结构，逐步选出具有高分类信息丰度的特征基因。首先根据二八准则挑选出一级基因池，再计算差异表达水平指数，然后基于各个基因的距离值和基因密度值，绘制决策图；接着运用多重线性回归分析方法，自动确定聚类中心；同时构建深度基因表达预测网络，计算GGSI值(基因间敏感性指数)，剔除冗余基因，建立二级基因池；最后基于基因间敏感性指数，根据GGSI值进行布谷鸟搜索算法的二进制编码，选出最紧凑的特征基因，建立三级基因池。

本发明的有益效果主要表现在：设计的差异表达水平指数的计算方法，适用于多分类数据集的基因差异表达水平计算；采用多重线性回归分析拟合二元平面，可以自动确定聚类中心；基于深度神经网络学习基因间的敏感性信息，能有效剔除冗余基因，并根据二进制布谷鸟搜索算法查找最紧凑的特征基因集合，保留具有高分类信息丰度的特征基因。

附图说明

图1是基于深度学习和进化计算方法的特征基因选择方法框图。

图2是基于深度学习神经网络的深度基因表达预测网络整体框架模型示意图。

图3是基于基因间敏感性信息的二进制布谷鸟寻优算法流程图。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图3，一种基于深度学习和进化计算的特征基因选择方法，包括以下步骤：

1)选择差异表达的基因，建立一级基因池，过程如下：

1.1)计算原始基因池中各个基因的差异表达水平指数，即IIC-FC指数：

式(1)适用于多分类数据集的基因差异表达水平计算，其中c表示原始基因池中的基因个数，和分别表示基因i和基因j的表达水平均值，和分别表示基因i和基因j的表达水平标准差，函数max{·,·}和min{·,·}分别表示取最大值和最小值，ln(·)为对数函数；若IIC-FC指数越大，表明该基因包含的样本分类信息丰度越高，用于样本分类获得高精度的可能性越大；

1.2)根据二八准则，按照步骤1.1)中得到的IIC-FC值从高到低对原始基因池中的所有基因进行排序，依据排序选择前20％的基因作为一级基因池；

2)对一级基因池中的基因进行自动聚类，过程如下：

2.1)基于密度聚类算法计算一级基因池中每个基因的距离值δ和密度值ρ；定义欧拉距离作为基因i和基因j之间的相似度距离，计算公式如下：

式(2)中，N_sam表示样本个数，gE(i,k)表示基因i在样本k中的表达值；

首先计算各个基因的密度值ρ，得到密度矩阵计算公式如下：

其中

式(3)中，ρ_i表示基因gene_i的密度值，一级基因池为相应的指标集为I_P＝{1,2,...,N_pool1}；其中定义表达值累加和最大处的基因gene_max与最小处的gene_min之间的相似度距离的2％作为d_c值，计算公式如下：

d_c＝0.02*d_gene(gene_max,gene_min)>

然后计算各个基因的距离值，得到距离矩阵每个基因gene_i的距离值定义为δ_i，首先查找比gene_i密度大的基因集，得到集合P'＝{gene_j}，然后查找P'中与gene_i的距离最近的基因，则得到δ_i＝d_gene(i,j')；

定义函数f_γ是关于变量ρ和δ的二元离散函数，对应于三维空间中的坐标值是(ρ,δ,f_γ)，则得到双变量离散函数为：

式(5)中，f_γ取密度值和距离值的乘积的对数值作为函数值；表示大于零的较小正数(一般取)，其目的是为了在密度为零时，即没有点落在以d_c为半径的圈内时式子仍能够成立，并没有实际物理意义；

2.2)根据步骤2.1)中得到的双变量离散函数γ＝f_γ(ρ,δ)，利用最小二乘法进行二元斜面的拟合，得到拟合平面为z_γ＝b₁+b₂ρ+b₃δ，计算各个数据点的残差值为ε_γi＝y_γi(ρ,δ)-γ_i(ρ,δ)，绘制残差直方图ε_γi-h，并根据钟型曲线的正态拟合得到方差值为σ_γ，利用3σ原则自动把处在置信区间外的η个奇异点确定为聚类中心，并记为c_γ；

3)构建深度基因表达预测网络(DGEPN)，网络框架示意图如图2所示，主要包括输入层、隐藏层、输出层；然后计算一级基因池的“基因-基因”敏感性信息(GGSI)，过程如下：

3.1)区分样本的训练集、验证集和测试集，搭建多层神经网络作为深度基因表达预测网络；首先假定有N_sam个样本，每个样本在原始基因池中有N_gene维基因表达，其中M_in个基因作为神经网络的输入层，即将一级基因池中的所有基因作为输入；M_out个基因作为神经网络的输出层，即将原始基因池中对应一级基因池的补集作为输出；若满足M_in＝N_pool1，M_in+M_out＝N_gene，则每个样本可以表示为然后随机挑选N_train个样本组成训练集，N_val个样本组成验证集，N_test个样本作为测试集，其中N_train:N_val:N_test＝4:1:1；

在神经网络的输入层和输出层间设有l_hidden层隐藏层(一般取l_hidden＝3)，隐藏层的神经元之间进行全连接；第i层隐藏层的神经元个数NumNeu_i为：

NumNeu_i＝M_in+i*Z{abs(M_in-M_out)/l_hidden+1}>

式(6)中，函数abs(·)表示取绝对值，Z{·}表示取整函数；为了减少训练参数从而加快训练，对神经网络中的隐藏层设置dropout比率(一般取dropout＝0.25)；激活函数采用ReLU函数，解决梯度弥散问题；

3.2)计算基因间的敏感性信息，基于步骤3.1)中搭建的深度基因表达预测网络，定义第i个输入基因到第j个输出基因的敏感性信息为：

式(7)中，表示第j个输出神经元对第i个输入神经元的导数；敏感性信息指数越大，说明第i个基因对第j个基因的非线性相关性越高；

3.3)根据步骤3.2)中得到的基因间敏感性信息，计算神经网络输入层中每个输入基因归一化后的敏感性指数：

3.4)训练神经网络，定义输出值与期望值的平均平方误差作为损失函数，即：

当训练代数达到预设的最大值或者损失值小于预设的阈值，结束训练，并用测试集的平均绝对误差(MAE)作为神经网络模型的评价指标：

误差越小，说明计算GGSI的网络模型性能越好；

4)建立二级基因池，过程如下：

4.1)计算步骤2)中得到的每个基因簇的最大GGSI值和平均其中K是基因簇的个数；然后依据簇间剔除和簇内剔除策略构建二级基因池；

4.2)首先进行簇间剔除，若第j个簇的GGSI最大值远小于其它簇的GGSI最大值，即则剔除该基因簇中的所有基因，保留其他基因簇；

4.3)然后进行簇内剔除，若第i个簇内的第j个基因的GGSI值小于该簇的GGSI平均值，即则剔除该基因，保留该簇内其他基因，其中NumK_i表示基因簇i中包含的基因个数；

5)基于二进制布谷鸟搜索算法进行优化搜索最紧凑的基因集合，建立三级基因池，过程如下：

5.1)进行鸟巢的编码初始化，首先将二级基因池中的每个基因根据其GGSI值从高到低进行排序，然后对第i个鸟巢的第j维，即第i个解的第j个基因进行初始化编码：

式(11)中，d表示解的维度，即二级基因池中的基因个数；式(11)表示GGSI排名靠前20％的基因和排名靠后20％的基因分别编码为1和0，剩余基因进行0或者1的随机编码；

5.2)进行糟糕解的替换，为了保证新产生的解都是由优秀的个体组成，基于概率p_a(一般取p_a∈[0,1])淘汰不适应环境的糟糕解，并将它们替换为新的解决方案，替换公式为：

式(12)中，d表示解的维度；式(12)表示对某个糟糕解，将其GGSI排名靠前10％的基因和排名靠后10％的基因分别编码为1和0，剩余基因的编码不变；

5.3)进行鸟巢位置的更新，更新公式为：

式(13)中，Sig(·)表示sigmoid函数，Avg(·)表示二级基因池中所有基因的GGSI的平均值，其中Levy(λ)～u＝s^-λ，λ∈(1,3]，α＞0表示步长尺度缩放因子；

5.4)布谷鸟搜索优化算法的适应度函数定义为：

fitness(net_i)＝MAE(net_i)+ξ*GeneNumber(net_i)>

式(14)中，MAE(net_i)表示第i个解在二级基因池中所选择的基因集合在DGEPN中预测剩余基因表达的平均绝对误差，GeneNumber(net_i)表示第i个解所选的基因个数，ξ∈{10^-2,10^-3,10^-4}是一个尺度规模参数，ξ的取值根据GeneNumber(net_i)的值确定；

5.5)进行优化搜索迭代，当达到预设的迭代次数，或者适应度函数值小于预设的阈值，则停止迭代，将此时的最优解进行解码，得到最优三级基因池。寻优算法的流程图如图3所示，具体搜索算法的步骤如下：

参数定义:鸟巢个数m，每个鸟巢维度d，每个鸟巢的适应度值全局适应度值globalfit，迭代次数T，每个解的平均绝对误差基因池尺度缩放参数ξ，步长尺度缩放因子α，每一代的最小适应度值minfit以及对应的索引minIndex，全局最优解替换概率p_a.

开始

对于每个鸟巢do

对于每个维度j(j＝1,2,...,d)，do

根据式(11)初始化鸟巢

fit(i)＝+∞

globalfit＝+∞

对于每次迭代t(t＝1,2,...,T)，do

对于每个鸟巢do

训练DGEPN并将平均绝对误差保存为MAE(i)

Ifthen

对于每个维度j(j＝1,2,...,d)，do

[minfit,minIndex]←min{fit(i)|i∈(1,n)}

If(minfit<globalfit)，then

globalfit＝minfit

对于每个维度j(j＝1,2,...,d)，do

对于每个鸟巢do

对于每个维度j(j＝1,2,...,d)，do

依据概率p_a∈[0,1]选择糟糕解

根据式(12)替换糟糕解

对于每个鸟巢do

对于每个维度j(j＝1,2,...,d)，do

根据式(13)更新鸟巢位置结束。