一种基于深度自动编码器的基因表达数据聚类方法

1.本发明属于基因表达数据分析领域，具体涉及一种基于深度自动编码器的基因表达数据聚类方法。


背景技术：

2.随着基因芯片检测技术日趋成熟，产生了大量的基因表达数据，人们可从公共数据集中获取基因表达数据。利用数据挖掘方法挖掘基因表达数据重要信息，可从分子层面探讨疾病的机理和生物标志物。聚类分析是数据挖掘的重要技术之一，在基因表达数据研究中具有很重要的作用。基因表达数据的聚类分析研究可分为3类：基因聚类、样本聚类和双向聚类。其中，基因聚类以基因为对象，以样本表达值为特征，能够用于判断差异基因在不同情况下的表达模式，将表达模式相同或相近的基因聚集，这些同类基因可能具有相似的功能，共同参与同一代谢过程或存在于同一细胞通路中。
3.基因聚类可将表达模式相近的基因聚在一起，以发现未知基因功能以及基因之间的调控关系。现有技术中有基于划分的聚类算法、基于层次的聚类算法、基于密度的聚类算法以及基于图的聚类算法等经典的聚类算法，这些算法都存在着各自的优缺点。目前，应用较多的基因表达数据的基因聚类算法包括基于划分的k-means、基于层次的凝聚层次和birch、混合高斯模型(gmm)以及基于图的谱聚类，也有一些研究针对聚类算法的缺点进行改进的用于基因表达数据的聚类算法，但几乎都需指定聚类簇数，因此使用的都是模拟数据集或已知聚类簇数的数据集。k-means是一种经典的基于划分的无监督静态聚类方法，相比凝集层次聚类算法更快速更简单，且适用于多种数据集并有不错的聚类效果，但也需要给定聚类簇数，而现实中基因表达数据集基因聚类簇数往往是未知的，这成为了聚类分析的难点。r语言的mclust包专用于确定聚类簇数，内有14种算法用于选取聚类簇数。然而，将r语言的mclust包用于基因表达数据集基因聚类时却常因基因表达数据集的维数高计算量大而不能得到结果。
4.综上所述，亟需一种不需要给定聚类簇数且计算量较少的基因表达数据聚类方法。


技术实现要素：

5.针对现有技术存在的不足，本发明提出了一种基于深度自动编码器的基因表达数据聚类方法，该方法包括：
6.s1：获取基因表达数据；
7.s2：对基因表达数据进行预处理，得到mrna数据；
8.s3：将mrna数据输入到训练好的深度自动编码器中，得到降维后的mrna数据；
9.s4：采用mclust包确定降维后的mrna数据的聚类簇数；
10.s5：根据聚类簇数，采用k-means聚类算法对降维后的mrna数据进行聚类，得到聚类结果。
11.优选的，对获取的基因表达数据进行预处理包括：从基因表达数据中提取mrna数据；对mrna数据进行基因id转换、缺失值处理和删除重复基因，得到得到初步处理好的mrna数据；对初步处理好的mrna数据进行剔除过低表达量基因处理；使用方差稳定变换对剔除过低表达量基因后的数据进行标准化处理，得到预处理好的mrna数据。
12.优选的，对深度自动编码器进行训练包括：深度自动编码器包括编码器和解码器；
13.s31：采用编码器对输入mrna数据进行特征学习，得到降维后的mrna数据；
14.s32：采用解码器对降维后的mrna数据进行重构，得到重构的mrna数据；
15.s33：根据重构的mrna数据和原mrna数据计算损失函数，返回步骤s31，当损失函数最小时，得到训练好的深度自动编码器。
16.进一步的，深度自动编码器的编码器处理数据的公式为：
17.y＝f
θ
(x)＝fk(f
k-1
(...f1(x)))＝σk(σ
k-1
(...σ1(wx+b)))
18.其中，x表示输入向量，y表示降维后的mrna数据，f
θ
表示编码器，fi表示深度网络中第i层的映射函数，σk表示第k层的激活函数，w表示权重矩阵值，b表示偏移矩阵。
19.进一步的，深度自动编码器的解码器处理数据的公式为：
20.z＝g
θ
′
(y)＝gk(g
k-1
(...g1(y)))＝σk(σ
k-1
(...σ1(w
′
y+b
′
)))
21.其中，y表示输入的经过降维的mrna数据，z表示重构基因向量，g
θ
′
＝{g1,g2,...,gk}表示解码器，σk表示激活函数，w
′
表示权重矩阵值，b
′
表示偏移矩阵。
22.进一步的，损失函数为：
23.l＝-logp(x|z)＝-logp(x|z＝g
θ
′
(f
θ
(x)))
24.其中，l表示损失函数，x表示输入的原mrna数据，z表示输入x的一个重构基因向量，f
θ
(x)表示经过降维的mrna数据，g
θ
′
表示解码器。
25.优选的，采用mclust包确定降维后的mrna数据的聚类簇数包括：设定初始聚类簇数的范围；在聚类簇数的范围内，mclust包中所有算法均根据不同的聚类簇数计算bic值；将使得bci值最大时的聚类簇数作为mrna数据的聚类簇数。
26.进一步的，bic值的计算公式为：
[0027][0028]
其中，m表示mclust包中第m个算法，g表示聚类簇数，x表示输入的经过降维的mrna数据，表示极大似然估计，n表示基因数目，v表示估计参数个数。
[0029]
优选的，采用k-means聚类算法对降维后的mrna数据进行聚类包括：
[0030]
s51：定义初始聚类中心；
[0031]
s52：根据聚类中心对mrna数据进行聚类，得到聚类的mrna数据；
[0032]
s53：根据所有mrna数据与聚类中心的距离计算目标函数；
[0033]
s54：更新聚类中心，返回步骤s52，当目标函数最小时，停止聚类。
[0034]
进一步的，目标函数为：
[0035][0036]
其中，si表示第i个簇，μi表示簇si的均值向量聚类中心，k表示聚类簇数，x表示mrna数据。
[0037]
本发明的有益效果为：本发明采用mclust包确定聚类簇数，避免了人为设置聚类簇数的局限性，聚类结果更准确；将深度自动编码器与聚类方法相结合，深度自动编码器通过特征学习降低数据维度，加快了聚类运算速度，且有效地提升了聚类的效果；本发明不需要给定聚类簇数且计算量较少，聚类结果精度高，高精度基因聚类结果有助于更好的研究基因的未知功能以及基因之间的调控关系，具有良好的应用前景。
附图说明
[0038]
图1为本发明中基于深度自动编码器的基因表达数据聚类方法流程图；
[0039]
图2为本发明中一种优选实施例的甲状腺癌基因表达数据集预处理后的效果示意图；
[0040]
图3为本发明中深度自动编码器的网络模型。
具体实施方式
[0041]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0042]
本发明提出了一种基于深度自动编码器的基因表达数据聚类方法，如图1所示，所述方法包括：
[0043]
s1：获取基因表达数据；
[0044]
s2：对基因表达数据进行预处理，得到mrna数据；
[0045]
s3：将mrna数据输入到训练好的深度自动编码器中，得到降维后的mrna数据；
[0046]
s4：采用mclust包确定降维后的mrna数据的聚类簇数；
[0047]
s5：根据聚类簇数，采用k-means聚类算法对降维后的mrna数据进行聚类，得到聚类结果。
[0048]
一种基于深度自动编码器的基因表达数据聚类方法的具体实施例如下：
[0049]
采用癌症基因组图谱(tcga)数据库中甲状腺癌的基因表达数据集，对数据集进行mrna数据提取，对提取到的mrna数据进行基因id转换、缺失值处理和删除重复基因，得到初步处理后的mrna数据，此时基因数量由56404减少为19193个。
[0050]
对初步处理后的mrna数据进行进一步处理，包括：采用cpm归一化后剔除过低表达量基因，保留的基因为至少在两个样本里cpm归一化后表达量大于0.5，并使用方差稳定变换(vst)对数据进行标准化处理，最后基因表达数据矩阵为16971x568，其中基因个数为16971，样本数为568，表达量数据范围为6.29-23.38；预处理后的结果如图2所示。
[0051]
如图3所示，深度自动编码器网络包括输入层、隐藏层以及输出层；隐藏层包括5个全连接层；深度自动编码器的编码器包括输入层和隐藏层的3个全连接层，深度自动编码器的解码器包括输出层和隐藏层的2个全连接层，编码器与编码器均使用relu函数作为激活函数；其中，编码器第一层全连接层神经元个数为128，第二层全连接层为64，第三层全连接层为25，第三层的全连接层数量实质上是将输入数据压缩到这个维度；解码器全连接层可实现对编码器提取的特征进行重建并与原始输入数据进行对比，解码器第一个全连接层神
经元个数为64，第二个全连接层为128；输出层神经元个数为568，与原始数据大小保持一致，用于与输入基因表达数据进行对比从而迭代优化模型。
[0052]
对深度自动编码器进行训练包括：深度自动编码器包括编码器和解码器；
[0053]
s31：采用编码器对输入mrna数据进行特征学习，得到降维后的mrna数据；
[0054]
s32：采用解码器对降维后的mrna数据进行重构，得到重构的mrna数据；
[0055]
s33：根据重构的mrna数据和原mrna数据计算损失函数，返回步骤s31，当损失函数最小时，得到训练好的深度自动编码器。
[0056]
深度自动编码器的编码器处理数据的公式为：
[0057]
y＝f
θ
(x)＝fk(f
k-1
(...f1(x)))＝σk(σ
k-1
(...σ1(wx+b)))
[0058]
其中，x为输入向量，表示一个基因在各个样本的表达量，y为此基因的特征表示，此处为降维后的mrna数据，f
θ
＝{f1,f2,...,fk}为编码器，fi为深度网络中第i层的映射函数，σk表示第k层的激活函数，w为权重矩阵值即该基因提取的各个特征占各个样本的权重值，b是偏移矩阵；在这个神经网络中，第一层的输出作为第二层的输入，k个单一的隐藏层堆叠起来作为一个深度神经网络来拟合这个映射函数f
θ
。
[0059]
深度自动编码器的解码器处理数据的公式为：
[0060]
z＝g
θ
′
(y)＝gk(g
k-1
(...g1(y)))＝σk(σ
k-1
(...σ1(w
′
y+b
′
)))
[0061]
其中，y表示输入的经过降维的mrna数据，z表示重构基因向量，g
θ
′
＝{g1,g2,...,gk}表示解码器，σk表示第k层的激活函数，w
′
表示权重矩阵值，b
′
表示偏移矩阵。
[0062]
损失函数为：
[0063]
l＝-logp(x|z)＝-logp(x|z＝g
θ
′
(f
θ
(x)))
[0064]
其中，l表示损失函数，x表示输入原mrna数据，z表示输入x的一个重构基因向量，f
θ
(x)表示经过降维的mrna数据，g
θ
′
表示解码器。
[0065]
采用训练好的深度自动编码器对预处理后的mrna数据提取特征，得到降维后的mrna数据。
[0066]
采用mclust包确定降维后的mrna数据的聚类簇数，包括：设定初始聚类簇数的范围；在聚类簇数的范围内，mclust包中所有算法均根据不同的聚类簇数计算bic值；将使得bci值最大时的聚类簇数作为mrna数据的聚类簇数。
[0067]
bic值的计算公式为：
[0068][0069]
其中，m表示mclust包中第m个算法，g表示聚类簇数，x表示输入的经过降维的mrna数据，表示极大似然估计，n表示基因数目，v表示估计参数个数。
[0070]
对于给定的基因表达数据集d，d＝(x1,x2,x3,...,xn)，每个基因数据对象xi为m维向量，采用k-means聚类算法对降维后的mrna数据进行聚类，包括：
[0071]
s51：定义初始聚类中心；
[0072]
s52：根据聚类中心对mrna数据进行聚类，将n个基因数据对象划分为k簇，即s＝{s1,s2,s3,...,sk}，得到聚类的mrna数据；
[0073]
s53：根据所有mrna数据与聚类中心的距离计算目标函数；
[0074]
s54：更新聚类中心，返回步骤s52，当目标函数最小时，停止聚类。
[0075]
更新聚类中心的公式为：
[0076][0077]
其中，si表示第i个簇，μi表示簇si的均值向量聚类中心，x表示mrna数据。
[0078]
目标函数为：
[0079][0080]
其中，e表示簇中数据点与聚类中心之间的距离平方差之和，si表示第i个簇，μi表示簇si的均值向量聚类中心，k表示聚类簇数，x表示mrna数据。
[0081]
当距离度量为欧氏距离时，在误差平方函数的基础上对目标函数e进行极小化，当目标函数e收敛时，得到最佳聚类结果；由于聚集在一起的基因在所有样本中具有相似的表达模式，这可能表明这些基因的共同调控；此外，来自同一簇的基因可能执行类似的细胞功能，有助于注释新发现的基因，用户可通过对聚类结果的分析，实现从分子层面探讨疾病的机理和生物标志物，帮助用户获得更多的生物学信息。
[0082]
采用戴维森堡丁指数dbi、calinski-harabasz指标(ch指标)以及轮廓系数作为评价指标对本发明进行评价。戴维森堡丁指数dbi、ch指标以及轮廓系数为聚类算法的内部评价方法，dbi计算任意两类别的类内平均距离之和除以两聚类中心距离求最大值，然后求平均；ch指标通过计算类中各点与类中心的距离平方和来度量类内的紧密度，通过计算各类中心点与数据集中心点距离平方和来度量数据集的分离度，最终ch指标由分离度与紧密度的比值得到；轮廓系数根据欧几里得距离量化了聚类的紧密度和分离度。
[0083]
dbi表示为：
[0084][0085]
其中，n是类别个数，ci是第i个类别的中心，σi是类别i中所有的点到簇质心的平均距离；d(ci,cj)是两类别中心点之间的距离，dbi的值越小，意味着类内距离越小同时类间聚类越大，因此聚类效果越好。
[0086]
ch指标表示为：
[0087][0088]
其中，n为聚类数目，k为当前的类，是类间距离差矩阵的迹，是类内距离差矩阵的迹，z是整个数据集的均值，zj是第j个簇cj的均值，ch(k)值越大，聚类效果越好。
[0089]
轮廓系数表示为：
[0090][0091]
其中，s(i)表示第i个基因的轮廓系数；a(i)为第一簇内不相似度，表示第i基因到同簇内其它不同基因距离的平均值，体现紧密度；b(i)为第二簇间不相似度，表示第i个基
因到其它簇的所有基因平均距离的最小值，体现分离度。对于一个数据集合，它的轮廓系数是所有数据点(基因)的平均轮廓系数。轮廓系数的取值在[-1,1]之间，值越大聚类效果越好。
[0092]
将dbi、轮廓系数以及ch指数作为聚类方法的综合评价指标；常见的主成分分析(pca)与k-means结合的dbi、轮廓系数以及ch指数分别为1.26、23059.76、0.24；本发明的深度自编码器与k-means结合的dbi、轮廓系数以及ch指数分别为0.64、71963.67、0.46；与本发明相比，本发明的聚类效果更好。
[0093]
本发明采用mclust包确定聚类簇数，避免了人为设置聚类簇数的局限性，聚类结果更准确；将深度自动编码器与聚类方法相结合，深度自动编码器通过特征学习降低数据维度，加快了聚类运算速度，且有效地提升了聚类的效果；本发明不需要给定聚类簇数且计算量较少，聚类结果精度高，高精度基因聚类结果有助于更好的研究基因的未知功能以及基因之间的调控关系，具有良好的应用前景。
[0094]
以上所举实施例，对本发明的目的、技术方案和优点进行了进一步的详细说明，所应理解的是，以上所举实施例仅为本发明的优选实施方式而已，并不用以限制本发明，凡在本发明的精神和原则之内对本发明所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。