CN112967755B - 一种面向单细胞rna测序数据的细胞类型识别方法 - Google Patents

Abstract

本发明针对现有技术的局限性，提出了一种面向单细胞RNA测序数据的细胞类型识别方法，其采用概念分解提取了高维单细胞RNA测序数据的低维特征表示，获得了原始数据的主体信息，降低了数据噪声对模型性能的干扰；同时还采用自表示学习方法来从低维特征表示中学习细胞之间的相似性；因此，本发明的方案对噪声更为不敏感，且具有更强的适应能力；本发明还使用图正则化来整合利用多种子聚类算法的预测结果中蕴含的信息，使得提取的低维特征表示能契合子聚类结果中的关联结构，有助于更准确地从低维特征表示中学习出细胞之间的相似性。本发明改善了模型对原始数据的依赖，并且使得模型具有较好的抗噪性，能够得到更加稳健、准确的识别结果。

Description

一种面向单细胞RNA测序数据的细胞类型识别方法

技术领域

本发明涉及生物信息学技术领域，具体涉及生物数据挖掘；更具体的，涉及一种面向单细胞RNA测序数据的细胞类型识别方法。

背景技术

许多生物信息分析任务的展开需要知道细胞的类型，如细胞发育轨迹分析、肿瘤细胞亚型分析。但是细胞的真实类型往往是未知的，人工辨别细胞类型的方法一般依靠细胞表型特征如膜蛋白等元件的差异来进行区分，这些方法不仅速度慢、准确率低，还无法辨别细胞在遗传信息方面的差异。聚类分析算法凭借高速的运算能力以及先进的算法设计，成为了识别细胞类型的有效计算分析方法。

近年来，涌现了许多针对单细胞RNA测序数据的聚类算法，这些算法可以大致分为以下几类：基于矩阵分解的方法、基于图论的方法、基于概率模型的方法以及集成聚类方法。由于模型假设和优化策略的不同，这些算法有着各自的优点和缺陷。比如基于矩阵分解的方法把原始细胞样本投影到低维的特征空间中，以更具区分性的特征表示原始数据的样本点。不过矩阵分解模型只关注于缩小原始数据在高维空间中的重构误差，没有考虑数据在低维流形上的结构；而基于图论的方法通过特定的方法来构建细胞之间的关联图，如距离、核方法或近邻关系，然后通过图算法从关联图里提取聚类结构。这类方法根据图中的细胞关联度进行聚类分析，往往具有很好的抗噪性能。不过这类算法需要预定的图构建方法，但是对于不同类型的细胞数据集，很难确定哪种图构建方法能更准确的评估细胞之间的相似性；基于概率模型的方法按照先验知识假设单细胞基因表达数据的概率分布模型，再通过数据集学习模型中的参数，进而获得样本的类别概率。在拥有一定先验知识的前提下，基于概率模型的方法能取得比较好的性能。但在不熟悉数据分布的情况下很难确定样本的分布模型，限制了模型的准确性；集成聚类方法收集多种聚类算法的聚类结果并从中提炼更为准确的聚类结构，参与集成的聚类算法可以统称为子聚类算法，子聚类算法的聚类结果统称为子聚类结果。集成聚类算法可以整合利用各子聚类算法的优点，因此在不同类型的数据集上能取得比较稳定的聚类性能。但是集成聚类通常只考虑子聚类算法的结果，而忽视了原始数据所蕴含的信息，这就容易导致最终的结果过分依赖于子聚类算法的结果。

而公开日为2020.02.14，公开号为：CN110797089A的中国发明专利：一种基于单细胞RNA测序数据识别细胞类型的方法，试图通过将高维矩阵低秩表示模型和图正则化理论结合起来，同时考虑数据的全局结构和局部结构特征来构建优化模型，通过采用交替方向乘子法来求解模型得到可靠的细胞与细胞间相似性矩阵，然后采用谱聚类方法对相似性矩阵进行聚类，以此实现对单细胞进行聚类，识别细胞类型，提高单细胞RNA测序数据的聚类效果，但该方案仍无法解决上述技术问题。

发明内容

针对现有技术的局限，本发明提出一种面向单细胞RNA测序数据的细胞类型识别方法，本发明采用的技术方案是：

一种面向单细胞RNA测序数据的细胞类型识别方法，包括以下步骤：

S1，获取包括单细胞RNA测序数据的原始数据集，使用若干个子聚类算法分别对所述原始数据集进行聚类分析，获得对应的子聚类结果；

S2，运用概念分解的矩阵分解模型提取所述原始数据集的低维特征表示；

S3，构建用于描述细胞相似性的关联矩阵，通过自表示学习获取所述低维特征表示的结构信息，将所述结构信息保存到所述关联矩阵上；

S4，通过图正则化将所述子聚类结果整合到所述矩阵分解模型中，结合所述关联矩阵，构建目标函数；

S5，通过求解所述目标函数优化更新所述关联矩阵，通过对优化更新后的关联矩阵进行谱聚类获得细胞类型识别结果。

相较于现有技术，本发明采用概念分解方法，提取了高维单细胞RNA测序数据的低维特征表示，获得了原始数据的主体信息，降低了数据噪声对模型性能的干扰；同时还采用自表示学习方法来从低维特征表示中学习细胞之间的相似性：与直接从原始高维数据中学习到的细胞相似性相比，从低维特征表示中学习到的细胞相似性对噪声更为不敏感；相比直接利用特定的距离度量来计算细胞相似性，自表示学习能自适应地学习出细胞之间的相似性，使本发明的方案具有更强的适应能力；本发明还使用图正则化来整合利用多种子聚类算法的预测结果中蕴含的信息，使得提取到的低维特征表示能契合子聚类结果中的关联结构，有助于更准确地从低维特征表示中学习出细胞之间的相似性。本发明改善了模型对原始数据的依赖，并且使得模型具有较好的抗噪性，能够得到更加稳健、准确的识别结果。

作为一种优选方案，在所述步骤S2提取所述原始数据集的低维特征表示XW的过程中，所述原始数据集通过概念分解的矩阵分解模型进行拆解，目标式如下：

其中，X∈R^n×p为原始数据集对应的原始数据矩阵，n为所述单细胞RNA测序数据的样本数量，p为所述单细胞RNA测序数据的基因数量；W∈R^p×k为投影矩阵，V∈R^k×p为系数矩阵，k为特征数目。

进一步的，在结合所述关联矩阵后，所述目标式如下：

其中，λ₁为衡量低秩约束项重要性的超参数，diag(Z)为关联矩阵Z∈R^n×n的对角线元素。

进一步的，在所述步骤S4通过图正则化将所述子聚类结果整合到所述矩阵分解模型中的过程中的图正则化项公式如下：

/>

其中，L_S＝D-S为集成关联矩阵S∈R^n×n的拉普拉斯矩阵；D为度对角矩阵，D的对角线上元素的值为所述集成关联矩阵S各行元素之和；所述集成关联矩阵S通过以下公式获得：

S^r∈R^n×n为转化第r个子聚类算法结果获得的关联图，q为子聚类算法的总数量。

进一步的，在所述步骤S4通过图正则化将所述子聚类结果整合到所述矩阵分解模型中的过程中，通过以下公式将所述子聚类结果π^r＝{π^r(X₁),...,π^r(X_n)}转化为关联图S^r∈R^n×n：

其中，π^r(X_i)为第i个细胞在第r个子聚类结果中分配到的类标。

作为一种优选方案，所述目标函数如下：

其中，λ₁，λ₂以及λ₃分别是控制低秩约束项、图正则化项以及F范数约束项的超参数。

进一步的，在所述步骤S5各轮迭代更新的过程中，依次交替更新求解变量V、W、Z，直至所述目标函数收敛。

进一步的，在所述步骤S5各轮迭代更新的过程中，在求解变量V时，通过固定变量W、Z将所述目标函数简化为：

通过对L关于V求导，将导数置为0，获得闭式解如下：

进一步的，在所述步骤S5各轮迭代更新的过程中，在求解变量W时，通过

固定变量V、Z将所述目标函数简化为：

通过对L关于W求导，获得关于W的梯度如下：

进一步的，在所述步骤S5各轮迭代更新的过程中，在求解变量Z时，通过固定变量V、W将所述目标函数简化为：

对上式运用奇异阈值运算符进行优化获得的更新式如下：

附图说明

图1为本发明实施例提供的面向单细胞RNA测序数据的细胞类型识别方法的步骤流程图；

图2为本发明实施例验证实验使用数据集的统计信息；

图3为本发明实施例验证实验中图正则化自表示集成聚类算法(GSREC)的超参数λ₁和λ₂在SimKumar8hard数据集上基于ACC评估指标的调参测试结果；

图4为本发明实施例验证实验中图正则化自表示集成聚类算法(GSREC)与用于集成的子聚类算法基于ACC评估指标在不同数据集上的对比结果；

图5为本发明实施例验证实验中图正则化自表示集成聚类算法(GSREC)与其它集成聚类算法基于ACC评估指标在不同数据集上的对比结果。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；

应当明确，所描述的实施例仅仅是本申请实施例一部分实施例，而不是全部的实施例。基于本申请实施例中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本申请实施例保护的范围。

在本申请实施例使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本申请实施例。在本申请实施例和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。还应当理解，本文中使用的术语“和/或”是指并包含一个或多个相关联的列出项目的任何或所有可能组合。

下面的描述涉及附图时，除非另有表示，不同附图中的相同数字表示相同或相似的要素。以下示例性实施例中所描述的实施方式并不代表与本申请相一致的所有实施方式。相反，它们仅是如所附权利要求书中所详述的、本申请的一些方面相一致的装置和方法的例子。在本申请的描述中，需要理解的是，术语“第一”、“第二”、“第三”等仅用于区别类似的对象，而不必用于描述特定的顺序或先后次序，也不能理解为指示或暗示相对重要性。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本申请中的具体含义。

此外，在本申请的描述中，除非另有说明，“多个”是指两个或两个以上。“和/或”，描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况。字符“/”一般表示前后关联对象是一种“或”的关系。以下结合附图和实施例对本发明做进一步的阐述。

为了解决现有技术的局限性，本实施例提供了一种技术方案，下面结合附图和实施例对本发明的技术方案做进一步的说明。

一种面向单细胞RNA测序数据的细胞类型识别方法，请参考图1，包括以下步骤：

S1，获取包括单细胞RNA测序数据的原始数据集，使用若干个子聚类算法分别对所述原始数据集进行聚类分析，获得对应的子聚类结果；

S2，运用概念分解的矩阵分解模型提取所述原始数据集的低维特征表示；

S3，构建用于描述细胞相似性的关联矩阵，通过自表示学习获取所述低维特征表示的结构信息，将所述结构信息保存到所述关联矩阵上；

S4，通过图正则化将所述子聚类结果整合到所述矩阵分解模型中，结合所述关联矩阵，构建目标函数；

S5，通过求解所述目标函数优化更新所述关联矩阵，通过对优化更新后的关联矩阵进行谱聚类获得细胞类型识别结果。

相较于现有技术，本发明采用概念分解方法，提取了高维单细胞RNA测序数据的低维特征表示，获得了原始数据的主体信息，降低了数据噪声对模型性能的干扰；同时还采用自表示学习方法来从低维特征表示中学习细胞之间的相似性：与直接从原始高维数据中学习到的细胞相似性相比，从低维特征表示中学习到的细胞相似性对噪声更为不敏感；相比直接利用特定的距离度量来计算细胞相似性，自表示学习能自适应地学习出细胞之间的相似性，使本发明的方案具有更强的适应能力；本发明还使用图正则化来整合利用多种子聚类算法的预测结果中蕴含的信息，使得提取到的低维特征表示能契合子聚类结果中的关联结构，有助于更准确地从低维特征表示中学习出细胞之间的相似性。本发明改善了模型对原始数据的依赖，并且使得模型具有较好的抗噪性，能够得到更加稳健、准确的识别结果。

因此，可以将本发明实施例所提供的方法，视为一种新的图正则化自表示集成聚类算法(GSREC)，其用于从高维度、高缺失的单细胞RNA测序数据中挖掘出细胞聚类结构。

具体的，利用矩阵分解提取原始单细胞RNA测序数据的低维特征表示，同时通过图正则化约束低维特征表示契合多个子聚类结果整合得到的关联结构，在于使得提取的低维特征表示能保留原始数据和子聚类结果两种信息共同刻画的细胞关联结构，从而实现两种信息的有效整合。

自表示学习是一种基于研究对象自身的相似性学习方法，即从一组样本构成的子空间中重构出该组样本并学习出样本之间的相似性。

作为一种优选实施例，在所述步骤S2提取所述原始数据集的低维特征表示XW的过程中，所述原始数据集通过概念分解的矩阵分解模型进行拆解，目标式如下：

其中，X∈R^n×p为原始数据集对应的原始数据矩阵，n为所述单细胞RNA测序数据的样本数量，p为所述单细胞RNA测序数据的基因数量；W∈R^p×k为投影矩阵，V∈R^k×p为系数矩阵，k为特征数目。

具体的，在本实施例中，大写字母表示矩阵，小写字母表示标量。那么原始数据矩阵X的(i,j)项元素写成X_ij；X_i表示原始数据矩阵X的第i个行向量；X^T表示X的转置矩阵；tr(X)表示X的迹；||X||_F、||X||_*分别表示X的Frobenius范数和核范数；I_n表示尺寸为n行n列的单位矩阵。

在所述步骤S2中，首先将原始数据集的样本数据通过投影矩阵W∈R^p×k将原始数据投影到低维特征子空间，得到低维特征表示XW；再利用系数矩阵V∈R^k×p重组拟合原始数据。因为投影矩阵W的维度远小于基因的数量，由此实现了数据压缩的目的。获得的低维特征保留了原始数据中的主体信息。其中，特征数目k小于n和p，默认值为300，这样既能保证储存尽可能多的压缩信息，还能减少运算量。

进一步的，在结合所述关联矩阵后，所述目标式如下：

其中，λ₁为衡量低秩约束项重要性的超参数，diag(Z)为关联矩阵Z∈R^n×n的对角线元素。

具体的，在所述关联矩阵Z∈R^n×n中，刻画了第i个细胞与第j个细胞之间的相似性；由于属于相同类别的细胞的低维特征表示往往是比较相似的，在所述步骤S3中，为所述关联矩阵Z添加低秩约束，使得关联矩阵具有低秩结构。将diag(Z)约束为0可以令每个细胞的低维特征表示都通过其它细胞的低维特征表示来重构，避免无效解的产生。在低维特征子空间进行的自表示学习可以减轻原始数据的噪声干扰。

进一步的，在所述步骤S4通过图正则化将所述子聚类结果整合到所述矩阵分解模型中的过程中的图正则化项公式如下：

其中，L_S＝D-S为集成关联矩阵S∈R^n×n的拉普拉斯矩阵；D为度对角矩阵，D的对角线上元素的值为所述集成关联矩阵S各行元素之和；所述集成关联矩阵S通过以下公式获得：

S^r∈R^n×n为转化第r个子聚类算法结果获得的关联图，q为子聚类算法的总数量。

进一步的，在所述步骤S4通过图正则化将所述子聚类结果整合到所述矩阵分解模型中的过程中，通过以下公式将所述子聚类结果π^r＝{π^r(X₁),...,π^r(X_n)}转化为关联图S^r∈R^n×n：

其中，π^r(X_i)为第i个细胞在第r个子聚类结果中分配到的类标。δ(y₁,y₂)中的y₁,y₂仅表示该函数的输入值，即细胞在子聚类结果中分配到的类标，并非指引入了其它变量，特此说明。

具体的，本发明实施例的矩阵分解模型基于提取的特征来学习细胞之间的相似性；但是低维特征的提取比较依赖原始数据的质量，原始数据中的高噪声以及缺失问题都会干扰到细胞相似性的学习。

对于本发明实施例中用到的q个子聚类算法，第r个子聚类算法产生的聚类结果是π^r＝{π^r(X₁),...,π^r(X_n)}，由于S^r∈R^n×n刻画的是第r个聚类算法产生的聚类结果，并不能代表真实的结果；不过如果一对细胞在多个算法的聚类结果中都被划为同一类，则表明这对细胞有着较强的相似性。基于这种思想，本发明实施例通过图正则化引入所述子聚类结果中较稳定的细胞结构信息，整合多个子聚类结果所刻画的关联图S^r∈R^n×n，最终构建出聚类结果集的集成关联矩阵S∈R^n×n；所述集成关联矩阵中的元素值越大，意味着对应的细胞对属于同一类的可能性越高。

获得集成关联矩阵S后，本发明实施例利用这个矩阵提供的信息来指导低维特征的提取，使得提取的低维特征具有更清晰的聚类结构。通过使所述图正则化项公式最小化，则S_ij的数值越大，其对应的低维特征表示越相似。

作为一种优选的实施例，所述目标函数如下：

其中，λ₁，λ₂以及λ₃分别是控制低秩约束项、图正则化项以及F范数约束项的超参数。

具体的，F范数惩罚主要为了避免所述投影矩阵W和系数矩阵V的数值在迭代更新过程中变得过小。

由于所述目标函数中包括多个未知变量，因此，作为一种优选的实施例，在所述步骤S5各轮迭代更新的过程中，依次交替更新求解变量V、W、Z，直至所述目标函数收敛。

进一步的，在所述步骤S5各轮迭代更新的过程中，在求解变量V时，通过固定变量W、Z将所述目标函数简化为：

通过对L关于V求导，将导数置为0，获得闭式解如下：

进一步的，在所述步骤S5各轮迭代更新的过程中，在求解变量W时，通过固定变量V、Z将所述目标函数简化为：

通过对L关于W求导，获得关于W的梯度如下：

具体的，获得W相对于L的梯度，可以采用梯度下降方法降低目标函数。而梯度下降的难点在于更新步长的确定，作为一种优选实施例，可以参考AdaDelta方法，通过累积梯度和累积更新值来确定每次变量更新的步长。

进一步的，在所述步骤S5各轮迭代更新的过程中，在求解变量Z时，通过固定变量V、W将所述目标函数简化为：

对上式运用奇异阈值运算符进行优化获得的更新式如下：

具体的，对上式运用奇异阈值运算符进行优化，首先基于低维特征XW作奇异值分解XW＝UΛV^T；其中，Λ∈R^n×k是由XW的奇异值按从大至小顺序排列组成的对角矩阵，σ_i则是XW的第i个奇异值；再将奇异值按照阈值

把下标分为两个集合：

同时将对应的矩阵按列划分为两部分U＝[U₁,U₂],V＝[V₁,V₂],Λ＝diag[Λ₁,Λ₂]；最终得到上述关于变量Z的更新式。

进一步的，对于第t次迭代时所述目标函数的数值l^t，迭代停止的收敛条件为

以下将结合具体的验证实验对本实施的面向单细胞RNA测序数据的细胞类型识别方法进行说明；为便于实验结果的展示，后续说明中，将以缩写GSREC指代本实施的面向单细胞RNA测序数据的细胞类型识别方法：

实验使用数据集的统计信息请参阅图2，本实施例选取了5个来自小鼠、人类的单细胞RNA测序数据集作为基准数据集进行验证实验，保证了验证实验的全面性。由于原始数据含有较多的异常数值，分析前需要进行预处理，我们采用Seurat的策略来过滤数据。所述子聚类算法有6种，分别是CIDR、ascend、FlowSom、PCAKmeans、TSCAN以及SC3。

选取评估指标：在本实施例的验证实验中，采用评估预测准确度的准确率(ACC)评估指标来评价聚类结果的好坏。

参数设定及效果评估：低维特征W的维度k对结果的影响不大，不过要稍大于数据集的聚类数目c以保证特征矩阵能保存细胞的聚类信息。为了降低模型的运算量，降维维度k取得越低越好，本实施例的验证实验将其默认值设为300。此外，还要考察超参数对聚类性能的影响，两个超参数λ₁、λ₂分别在范围中进行调参测试。而超参数λ₃对聚类性能的影响不大，亦可设为默认值10^-2。具体的，本实施例的验证实验选取SimKumar8hard数据集进行调参测试，实验结果如图3所示。可见模型在超参数的调节范围内有一个比较明显的峰值，预测性能经过一定的上升后趋于平稳，这个规律说明模型对于数据集有一个相对合适的最优参数组合。且在最优参数范围里，模型的性能比较稳定没有出现波动，这说明模型性能对于参数不太敏感。对于不同的数据集，模型的最佳参数略有差异，建议将其分别设为1和10³。其他用于对比的集成聚类算法也会选用调参后的最优结果进行对比。

本实施例的验证实验将用于集成的子聚类算法结果与本发明实施例的提供的GSREC进行对比以说明GSREC的稳定性。结果如图4所示，可以看出，在各式各样的数据集中，最佳的子聚类算法会随着数据集的变化而改变，这是因为子聚类算法基于特定的策略分析数据，在不同类型的数据集中性能会出现波动。本发明实施例的提供的GSREC在所有数据集中都是高于或接近最佳的结果，这表明了GSREC不仅能降低数据噪声的干扰，还能稳定地从单细胞数据中挖掘结构信息，具有更准确的聚类性能。

为了进一步验证本发明的有效性，本实施例的验证实验选择了两个集成聚类算法进行比较，它们分别为LWEC和LWGP。实验结果如5所示，从图中可以看出，本发明实施例的提供的GSREC除了在KohTCC数据集上性能略低于LWGP算法，在其他数据集上都达到了最佳。在总体性能上，本发明实施例的提供的GSREC在平均得分上也是优于所有的对比算法。这进一步验证了GSREC相对于其他算法在模型设计上的优越性。这也证明了GSREC能挖掘高维数据背后隐藏的低维结构信息。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (3)

1.一种面向单细胞RNA测序数据的细胞类型识别方法，其特征在于，包括以下步骤：

S1，获取包括单细胞RNA测序数据的原始数据集，使用若干个子聚类算法分别对所述原始数据集进行聚类分析，获得对应的子聚类结果；

S2，运用概念分解的矩阵分解模型提取所述原始数据集的低维特征表示；

S3，构建用于描述细胞相似性的关联矩阵，通过自表示学习获取所述低维特征表示的结构信息，将所述结构信息保存到所述关联矩阵上；

S4，通过图正则化将所述子聚类结果整合到所述矩阵分解模型中，结合所述关联矩阵，构建目标函数；

S5，通过求解所述目标函数优化更新所述关联矩阵，通过对优化更新后的关联矩阵进行谱聚类获得细胞类型识别结果。

2.根据权利要求1所述的面向单细胞RNA测序数据的细胞类型识别方法，其特征在于，在所述步骤S2提取所述原始数据集的低维特征表示XW的过程中，所述原始数据集通过概念分解的矩阵分解模型进行拆解，目标式如下：

其中，X∈R^n×p为原始数据集对应的原始数据矩阵，n为所述单细胞RNA测序数据的样本数量，p为所述单细胞RNA测序数据的基因数量；W∈R^p×k为投影矩阵，V∈R^k×p为系数矩阵，k为特征数目。

3.根据权利要求2所述的面向单细胞RNA测序数据的细胞类型识别方法，其特征在于，在结合所述关联矩阵后，所述目标式如下：

其中，λ₁为衡量低秩约束项重要性的超参数，diag(Z)为关联矩阵Z∈R^n×n的对角线元素。

Images