CN113889192A - 一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法 - Google Patents

Abstract

本发明公开了一种基于深层降噪自编码器的单细胞RNA‑seq数据聚类方法，首先调整单细胞RNA‑seq数据的批次效应并对数据进行标准化从而减少技术噪声造成的不利影响；其次使用基于零膨胀负二项分布的深层降噪自编码器有效挖掘单细胞RNA‑seq数据的特征信息；然后使用快速独立成分分析方法对单细胞RNA‑seq数据进行降维，提高方法模型的计算效率；最后通过基于期望最大化的高斯混合模型对细胞展开更加准确的聚类，并使用T分布随机近邻嵌入的方法对最终的单细胞RNA‑seq数据聚类结果进行可视化。本发明能够有效降低单细胞RNA‑seq数据维度高、噪声大等特点对数据聚类的干扰，准确学习单细胞RNA‑seq数据的基因表达信息进而聚类细胞，为基因网络构建、发现细胞类型、癌症早期发现与治疗提供帮助。

Description

一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法

技术领域

本发明属于生物信息学中的单细胞RNA-seq数据分析技术领域，具体涉及一种基于深层降噪自编码器的单细胞RNA-seq(RiboNucleicAcid-sequencing)数据聚类方法。

背景技术

随着测序技术的快速发展，研究人员获得了大量的单细胞RNA-seq数据。无监督聚类对单细胞RNA-seq数据的分析有着重要的作用，针对单细胞RNA-seq数据的聚类方法不仅可以识别未知的细胞类型，还能揭示细胞的异质性。通过对单细胞RNA-seq数据的聚类方法研究，研究人员可以更加准确地鉴定细胞状态、搭建细胞之间的网络结构、深入理解癌症细胞等的分化过程，为未来癌症的早期发现和治疗奠定基础。目前，层次聚类、谱聚类和具有噪声的基于密度的聚类方法等传统聚类方法已经被广泛使用，但是，单细胞RNA-seq数据具有其独特的特征，这使得传统的聚类方法无法对这些数据进行有效聚类。

发明内容

为了克服以上技术问题，本发明提供了一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，将自编码器和快速独立成分分析等方法结合起来以实现单细胞RNA-seq数据聚类过程中特征学习、降维等目的，最终使用高斯混合聚类对单细胞RNA-seq数据进行聚类，并通过引入零膨胀负二项分布重构数据来降低数据噪声对聚类效果的影响。

为了实现上述目的，本发明采用的技术方案是：

一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，包括以下步骤；

1)调整批次效应与数据标准化预处理：

选用从ArrayExpress和GEO数据库中下载的5个公开的真实单细胞RNA-seq数据集来对单细胞进行聚类，进一步发现细胞亚型，为相关癌症的早期发现和针对性治疗提供帮助，这5个公开数据集中的基因表达值取自于各种组织细胞，包括GSE60361、GSE65525、GSE72056、GSE76312和GSE103322，读取原始的单细胞RNA-seq数据并对其进行批次效应调整和标准化预处理，避免因为样本数据在不同的批次处理和测量时引入的与生物状态不相关的系统性技术偏差；

2)数据重构与降噪：

由于单细胞RNA-seq数据中存在着大量的零值，这些零值既可能表示某些细胞的部分基因确实没有表达，也可能是技术误差导致的结果，这些噪声将会极大地干扰细胞亚型的发现，将经过对数标准化处理后的单细胞RNA-seq数据输入到深层降噪自编码器中，深层降噪自编码器使用零膨胀负二项分布来重构数据，重构出的数据能够较好地保存生物的原始特征；

3)数据降维：

深层降噪自编码器重构出的单细胞RNA-seq数据仍然是高维的，高维的单细胞RNA-seq数据给细胞亚型的鉴定带来了非常大的困难，使用快速独立成分分析的方法降低样本数据的维度，消除数据中的冗余部分，进一步避免因为数据中的冗余部分干扰癌症的早期发现和相关治疗；

4)高斯混合聚类与数据可视化：

得到低维低噪声的单细胞RNA-seq数据之后，使用高斯混合模型对细胞进行聚类并确定细胞类型，得到的细胞类型即为发现的潜在细胞亚型，采用T分布随机近邻嵌入方法可视化最终的聚类结果，并结合现有的细胞及癌症数据库上对聚类结果进行分析，帮助医生开展癌症的早期发现。

所述步骤1)中对单细胞RNA-seq数据调整批次效应和标准化预处理的步骤包括：首先使用层次贝叶斯模型来调整单细胞RNA-seq数据的批次效应，同时解决测量灵敏度带来的不确定性问题；接着筛选出基因表达量正常的细胞；然后，采用对数标准化方法对数据的测序深度和基因长度进行标准化。

所述步骤2)中所使用的深层降噪自编码器通过零膨胀负二项分布来重构单细胞RNA-seq数据，整个自编码器有三个输出，分别学习零膨胀负二项分布的零膨胀因子、均值和方差；

所述的待分析单细胞RNA-seq数据使用X来表示，自编码器中编码阶段表示为h(X)＝σ_h(WX+b)，W表示编码过程中的权重矩阵，b则表示偏置项，而自编码器的解码阶段则和编码阶段相对应，将编码后的数据进行重构，自编码器的输入维度和用于训练的单细胞RNA-seq数据的维度一致，编码器和解码器各有五层网络，在负二项分布(NB)模型的基础上加入一个零膨胀因子，也可以理解为在零点处增加一个脉冲函数，即用零膨胀负二项分布(Zero-Inflated Negative Binomial)来建模单细胞RNA-seq数据，公式化表示为ZINB(X|π,μ,θ)＝πδ₀(X)+(1-π)ΝΒ(X|μ,θ)，如果Y＝σ_o(W′h(X)+b′)表示解码器的最后一个隐层，在其后添加三个独立的全连接层，也就是说整个自编码器会有三个输出，分别学习零膨胀负二项分布的零膨胀因子、均值和方差，降噪自编码器的降噪部分的损失函数表示为L_d＝-log(ZINB(X|π,μ,θ))。

所述步骤3)中使用快速独立成分分析的方法降低单细胞RNA-seq数据的维度，独立成分分析假设所有数据的各部分彼此独立，并认为所有成分同等重要，将原始的数据分解成统计意义相互独立的非高斯数据成分的线性组合；

假设重构后的单细胞RNA-seq数据服从模型X＝AS，其中的S是具有独立分量的未知源数据，A是未知混合矩阵，S中的每个独立成分和A中的每个混合系数均是未知的，独立成分分析方法只通过X中每个观测到的信号数据来预估混合系数与独立分量，该方法会先对原始的数据进行中心化以及白化预处理，经过上述预处理后，采用快速独立成分分析的方法来处理样本数据，首先是初始化向量w，定义W＝A^-1，w是W中的行向量。其次，令w⁺＝E{Xg(w^TX)}-E{g′(w^TX)}w，上式中的g是一个非线性标量函数，并令w＝w⁺/||w⁺||，如果上述过程没有收敛，则不断重复这一步骤，最后，使用快速独立成分分析方法对包含重要信息的几个独立分量进行估计，实现降低单细胞RNA-seq数据维度的目的。

所述步骤4)中使用高斯混合模型对细胞进行聚类并确定细胞类型，具体步骤包括：

首先初始化高斯混合分布的模型参数，然后基于期望最大化算法反复迭代优化模型的参数；期望最大化算法中的E迭代步骤：基于第i个高斯混合分量计算出第i个样本数据的后验概率γ_ji：

；期望最大化算法中的M迭代步骤：反复迭代优化模型的其他参数μ_i，∑_i和α_i，基于以下公式计算：

当实验过程中达到最大的迭代次数时候停止迭代，如果不满足上述条件则继续迭代更新参数，最后，样本数据x_j的聚类标签λ_j使用λ_j＝argmaxγ_ji计算得到，使用T分布随机近邻嵌入方法对最终的聚类结果进行可视化，使聚类结果显示到二维坐标上。

所述的使用高斯混合模型进行初始化过程中，采用k-means++来解决质心初始化的问题，其方法是从输入的数据点集合中随机选择一个点作为第一个聚类中心；对与数据集中的每一个对象，计算它与最近聚类中心的相似性；选择一个新的数据点作为新的聚类中心，选择的原则是：相似性较大的点，被选取作为聚类中心的概率较大；重复上面的步骤直到k个聚类中心被选出，利用这k个初始的聚类中心来运行标准的k-means算法。

本发明的有益效果：

本发明将自编码器和快速独立成分分析方法结合起来学习单细胞RNA-seq数据的表示，并对数据进行降维，使用零膨胀负二项分布重构数据降低了数据噪声对聚类结果的影响，同时，使用高斯混合聚类对低维的单细胞RNA-seq数据进行聚类，最终使用T分布随机近邻嵌入方法可视化聚类结果，不仅可以鉴定细胞亚型，还能帮助开展癌症的早期发现和相关诊疗。

附图说明

图1为本发明总体流程图。

具体实施方式

下面结合实施例对本发明作进一步详细说明。

如图1所示，显示了本发明基于深层降噪自编码器提升单细胞RNA-seq数据聚类效果的四大步骤，调整批次效应与数据标准化预处理、数据重构与降噪、数据降维、高斯混合聚类与数据可视化。

本发明提供了一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，本发明包括以下步骤：

步骤一、调整批次效应与数据标准化预处理。本发明选用从ArrayExpress和GEO数据库中下载的5个公开数据集来验证本发明的有效性，这5个公开数据集中的基因表达值取自于各种组织细胞，包括GSE60361、GSE65525、GSE72056、GSE76312和GSE103322。这些数据作为初始输入的原始数据，本发明使用层次贝叶斯模型来调整单细胞RNA-seq数据的批次效应，同时解决测量灵敏度带来的不确定性问题。本发明使用一个基于Python的单细胞基因表达分析包SCANPY对原始的单细胞RNA-seq数据进行有效筛选和过滤，去除测序质量差的数据，随后将数据做标准化处理，以方便后续的网络学习。

步骤二、数据重构与降噪。由于数据扩增和数据丢失等问题，会干扰单细胞RNA-seq数据的分析。本发明使用降噪自编码器技术将输入的单细胞RNA-seq数据映射到嵌入空间。在实验过程中，对使用的数据加入随机高斯噪声，并用全连接层构建自编码器。为了更好地捕获单细胞RNA-seq数据，本发明在解码器最后一个隐层的后面添加三个独立的全连接层，三个输出分别学习零膨胀负二项分布的脉冲函数调节因子、负二项分布的均值和稀疏程度。进一步将降噪自编码器的降噪部分的损失函数定义为零膨胀负二项分布公式的负对数。待分析单细胞RNA-seq数据使用X来表示，则自编码器中编码阶段表示为h(X)＝σ_h(WX+b)，W表示编码过程中的权重矩阵，b则表示偏置项，而自编码器的解码阶段则和编码阶段相对应，将编码后的数据进行重构。自编码器的输入维度和用于训练的单细胞RNA-seq数据的维度一致，编码器和解码器各有五层网络。近年来针对单细胞RNA-seq数据的科研进展表明，单细胞RNA-seq数据最接近负二项分布(Negative Binomial，NB)，公式化表示为

因为单细胞RNA-seq数据的离散通常是高度扭曲的，方差往往会大于均值，因此不适合采用泊松分布来近似，而单细胞RNA-seq数据的方差通常会随着均值的改变而改变。单细胞RNA-seq数据除了上面的特点外，还有一个特点就是零值特别多。由于基因表达数据中的零值既可能来自生物过程中不表达的基因(TrueZero)，还可能来自测序过程中由于技术原因导致的丢失(Dropout Zero)。为了更好地捕获单细胞RNA-seq数据，本发明改进了传统的降噪自编码器，在负二项分布(NB)模型的基础上加入一个零膨胀因子，也可以理解为在零点处增加一个脉冲函数，即用零膨胀负二项分布(Zero-Inflated Negative Binomial)来建模单细胞RNA-seq数据。公式化表示为ZINB(X|π,μ,θ)＝πδ₀(X)+(1-π)ΝΒ(X|μ,θ)，如果Y＝σ_o(W′h(X)+b′)表示解码器的最后一个隐层，本发明在其后添加三个独立的全连接层，也就是说整个自编码器会有三个输出，分别学习零膨胀负二项分布的零膨胀因子、均值和方差。本发明将降噪自编码器的降噪部分的损失函数表示为L_d＝-log(ZINB(X|π,μ,θ))。

步骤三、数据降维。在数据降维的过程中，本发明首先对高维的单细胞RNA-seq数据进行中心化以及白化预处理，在预处理后的数据基础上计算出分离矩阵，并对分离矩阵进行初始化；然后不断优化分离矩阵，并一直判断是否收敛，若已经收敛，则可求出最终的低维单细胞RNA-seq数据，若没有收敛，则持续优化分离矩阵。

假设重构后的单细胞RNA-seq数据服从模型X＝AS，其中的S是具有独立分量的未知源数据，A是未知混合矩阵，S中的每个独立成分和A中的每个混合系数均是未知的。独立成分分析方法只通过X中每个观测到的信号数据来预估混合系数与独立分量。该方法会先对原始的数据进行中心化以及白化预处理，经过上述预处理后，本发明采用快速独立成分分析的方法来处理样本数，首先是初始化向量w，定义W＝A^-1，w是W中的行向量。其次，令w⁺＝E{Xg(w^TX)}-E{g′(w^TX)}w，上式中的g是一个非线性标量函数，并令w＝w⁺/||w⁺||，如果上述过程没有收敛，则不断重复这一步骤。最后，使用快速独立成分分析方法对包含重要信息的几个独立分量进行估计，实现降低单细胞RNA-seq数据维度的目的。

步骤四、高斯混合聚类与数据可视化。使用高斯混合模型进行初始化的过程中，采用k-means++来解决质心初始化的问题。首先初始化高斯混合分布的模型参数，然后基于期望最大化算法反复迭代优化模型的参数；期望最大化算法中的E迭代步骤：基于第i个高斯混合分量计算出第i个样本数据的后验概率γ_ji：

；期望最大化算法中的M迭代步骤：反复迭代优化模型的其他参数μ_i，∑_i和α_i，基于以下公式计算：

当实验过程中达到最大的迭代次数时候停止迭代，如果不满足上述条件则继续迭代更新参数。最后，样本数据x_j的聚类标签λ_j使用λ_j＝argmaxγ_ji计算得到。最后使用T分布随机近邻嵌入方法对聚类结果进行可视化。

所述的使用高斯混合模型进行初始化过程中，采用k-means++来解决质心初始化的问题，其方法是从输入的数据点集合中随机选择一个点作为第一个聚类中心；对与数据集中的每一个对象，计算它与最近聚类中心的相似性；选择一个新的数据点作为新的聚类中心，选择的原则是：相似性较大的点，被选取作为聚类中心的概率较大；重复上面的步骤直到k个聚类中心被选出，利用这k个初始的聚类中心来运行标准的k-means算法。

对于高斯混合模型，本发明使用scikit-learn混合模块中的高斯混合函数，其中的参数保持默认值。本发明使用的聚类性能评价指标主要包括Normalized MutualInformation(标准化互信息)，Clustering Accuracy(聚类准确率)，Adjusted Rand Index(调整兰德指数)，这三个评价指标的数值越高，表明方法的聚类性能越好。

本发明具有以下特点：

1.降低单细胞RNA-seq数据批次效应对最终聚类效果的影响；

2.降低单细胞RNA-seq数据维度高、噪声大对聚类结果的影响；

3.对单细胞RNA-seq数据聚类过程中，能够有效学习数据的表示，具备强大的数据表征能力；

4.聚类结束后有良好的数据可视化能力。

自编码器模型是一种无监督的深度学习方法，该方法不仅可以对输入数据进行有效降维，还可以通过调整神经网络层数、优化网络训练过程等方式来学习所分析数据的隐含特征，并通过数据重建过程对数据进行恢复。降噪自编码器允许带有噪声信息的损坏数据成为网络的输入数据，使重构数据对输入数据中的噪声获得一定鲁棒性。

单细胞RNA-seq数据不仅维度较高，还含有较大噪声，这些噪声通常表现为单细胞RNA-seq数据比较稀疏，大量的零值一方面来源于某些确实没有表达的基因，另一方面则来源于因为测序等技术的缺陷导致没有测出有表达的基因数值。

在对单细胞RNA-seq数据进行深度聚类的过程中，本发明设计了一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，将自编码器、快速独立成分分析、高斯混合聚类、T分布随机近邻嵌入等方法结合起来解决单细胞RNA-seq数据聚类过程中遇到的特征学习、降维、聚类和数据可视化等难题，并通过引入零膨胀负二项分布重构数据来降低数据噪声对聚类效果的影响。

Claims (6)

1.一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，其特征在于，包括以下步骤；

1)调整批次效应与数据标准化预处理：

选用从ArrayExpress和GEO数据库中下载的5个公开的真实单细胞RNA-seq数据集来对单细胞进行聚类，这5个公开数据集中的基因表达值取自于各种组织细胞，包括GSE60361、GSE65525、GSE72056、GSE76312和GSE103322，读取原始的单细胞RNA-seq数据并对其进行批次效应调整和标准化预处理；

2)数据重构与降噪：

将经过对数标准化处理后的单细胞RNA-seq数据输入到深层降噪自编码器中，深层降噪自编码器使用零膨胀负二项分布来重构数据，重构出的数据能够较好地保存生物的原始特征；

3)数据降维：

深层降噪自编码器重构出的单细胞RNA-seq数据仍然是高维的，高维的单细胞RNA-seq数据给细胞亚型的鉴定带来了非常大的困难，使用快速独立成分分析的方法降低样本数据的维度，消除数据中的冗余部分，进一步避免因为数据中的冗余部分干扰癌症的早期发现和相关治疗；

4)高斯混合聚类与数据可视化：

得到低维低噪声的单细胞RNA-seq数据之后，使用高斯混合模型对细胞进行聚类并确定细胞类型，得到的细胞类型即为发现的潜在细胞亚型，采用T分布随机近邻嵌入方法可视化最终的聚类结果，并结合现有的细胞及癌症数据库上对聚类结果进行分析，帮助医生开展癌症的早期发现。

2.根据权利要求1所述的一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，其特征在于，所述步骤1)中对单细胞RNA-seq数据调整批次效应和标准化预处理的步骤包括：首先使用层次贝叶斯模型来调整单细胞RNA-seq数据的批次效应，同时解决测量灵敏度带来的不确定性问题；接着筛选出基因表达量正常的细胞；然后，采用对数标准化方法对数据的测序深度和基因长度进行标准化。

3.根据权利要求1所述的一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，其特征在于，所述步骤2)中所使用的深层降噪自编码器通过零膨胀负二项分布来重构单细胞RNA-seq数据，整个自编码器有三个输出，分别学习零膨胀负二项分布的零膨胀因子、均值和方差；

所述的待分析单细胞RNA-seq数据使用X来表示，自编码器中编码阶段表示为h(X)＝σ_h(WX+b)，W表示编码过程中的权重矩阵，b则表示偏置项，而自编码器的解码阶段则和编码阶段相对应，将编码后的数据进行重构，自编码器的输入维度和用于训练的单细胞RNA-seq数据的维度一致，编码器和解码器各有五层网络，在负二项分布(NB)模型的基础上加入一个零膨胀因子，也可以理解为在零点处增加一个脉冲函数，即用零膨胀负二项分布(Zero-Inflated Negative Binomial)来建模单细胞RNA-seq数据，公式化表示为ZINB(X|π,μ,θ)＝πδ₀(X)+(1-π)ΝΒ(X|μ,θ)，如果Y＝σ_o(W′h(X)+b′)表示解码器的最后一个隐层，在其后添加三个独立的全连接层，也就是说整个自编码器会有三个输出，分别学习零膨胀负二项分布的零膨胀因子、均值和方差，降噪自编码器的降噪部分的损失函数表示为L_d＝-log(ZINB(X|π,μ,θ))。

4.根据权利要求1所述的一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，其特征在于，所述步骤3)中使用快速独立成分分析的方法降低单细胞RNA-seq数据的维度，独立成分分析假设所有数据的各部分彼此独立，并认为所有成分同等重要，将原始的数据分解成统计意义相互独立的非高斯数据成分的线性组合；

假设重构后的单细胞RNA-seq数据服从模型X＝AS，其中的S是具有独立分量的未知源数据，A是未知混合矩阵，S中的每个独立成分和A中的每个混合系数均是未知的，独立成分分析方法只通过X中每个观测到的信号数据来预估混合系数与独立分量，该方法会先对原始的数据进行中心化以及白化预处理，经过上述预处理后，采用快速独立成分分析的方法来处理样本数据，首先是初始化向量w，定义W＝A^-1，w是W中的行向量。其次，令w⁺＝E{Xg(w^TX)}-E{g′(w^TX)}w，上式中的g是一个非线性标量函数，并令w＝w⁺/||w⁺||，如果上述过程没有收敛，则不断重复这一步骤，最后，使用快速独立成分分析方法对包含重要信息的几个独立分量进行估计，实现降低单细胞RNA-seq数据维度的目的。

5.根据权利要求1所述的一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，其特征在于，所述步骤4)中使用高斯混合模型对细胞进行聚类并确定细胞类型，具体步骤包括：

首先初始化高斯混合分布的模型参数，然后基于期望最大化算法反复迭代优化模型的参数；期望最大化算法中的E迭代步骤：基于第i个高斯混合分量计算出第i个样本数据的后验概率γ_ji：

期望最大化算法中的M迭代步骤：反复迭代优化模型的其他参数μ_i，∑_i和α_i，基于以下公式计算：

当实验过程中达到最大的迭代次数时候停止迭代，如果不满足上述条件则继续迭代更新参数，最后，样本数据x_j的聚类标签λ_j使用λ_j＝arg maxγ_ji计算得到，使用T分布随机近邻嵌入方法对最终的聚类结果进行可视化，使聚类结果显示到二维坐标上。

6.根据权利要求5所述的一种基于深层降噪自编码器的单细胞RNA-seq数据聚类方法，其特征在于，所述的使用高斯混合模型进行初始化过程中，采用k-means++来解决质心初始化的问题，其方法是从输入的数据点集合中随机选择一个点作为第一个聚类中心；对与数据集中的每一个对象，计算它与最近聚类中心的相似性；选择一个新的数据点作为新的聚类中心，选择的原则是：相似性较大的点，被选取作为聚类中心的概率较大；重复上面的步骤直到k个聚类中心被选出，利用这k个初始的聚类中心来运行标准的k-means算法。

Images