CN110797089B - 一种基于单细胞rna测序数据识别细胞类型的方法 - Google Patents

Abstract

本发明提供了一种基于单细胞RNA测序数据识别细胞类型的方法。该方法有效地将高维矩阵低秩表示(Low Rank Representation)模型和图正则化理论结合起来，同时考虑数据的全局结构和局部结构特征来构建优化模型，通过采用交替方向乘子法(ADMM)来求解模型得到可靠的细胞与细胞间相似性矩阵，然后采用谱聚类方法对相似性矩阵进行聚类，从而实现对单细胞进行聚类，识别细胞类型。该方法可以显著提高单细胞RNA测序数据的聚类效果。

Description

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

技术领域

本发明涉及数学与生物学交叉研究领域，具体涉及通过聚类算法进行细胞分类的方法。

背景技术

传统测序技术得到的是一群细胞(肿瘤细胞、免疫细胞、成纤维细胞和巨噬细胞的混合物)基因表达值的平均化，这样忽视了细胞之间基因表达的差异性，很难鉴别细胞之间表达的异质性。

近些年来，随着生物技术的不断发展，单细胞RNA测序技术能得到单个细胞内大量基因的表达信息，为辨别生物组织中各种细胞类型的转录组特征和全面揭示细胞之间基因表达的异质性提供了非常有力的工具。与传统的全基因组测序相比，单细胞测序不仅测量基因表达水平更加精确，而且还能检测到微量的基因表达子或罕见非编码RNA，其优势是全方位和多层次的。

通过对单细胞RNA测序数据进行聚类，可以有效地将不同类型细胞进行区分，有助于识别新的细胞类型。另外，准确地将单细胞进行聚类，对后面重构细胞伪轨迹、识别显著表达基因以及揭示未知细胞状态为进一步理解异质性有着重要的意义。

然而，由于技术原因，目前单细胞RNA测序得到的数据具有高维度、高噪声、数据缺失、稀疏的特点，这给对单细胞RNA测序数据的聚类算法设计带来了巨大的困难和挑战。传统基于批量平均化RNA测序数据的方法无法胜任单细胞RNA测序数据的分析。因此，迫切需要一种有效的方法来对单细胞RNA测序数据进行聚类。近几年来国内外学者提出了很多基于单细胞RNA测序数据进行聚类的方法，例如，Xu et al.,等学者提出了一种名为SNN-Cliq的方法对细胞进行有效分类，首先采用共享最邻近邻居(SNN)方法来计算细胞之间的相似性，然后使用图论中聚类理论模型来实现对单细胞RNA测序数据进行聚类。Wang et al.等学者认为不同的核函数可以刻画数据的不同信息表示，提出了一种融合多个核函数的方法SIMLR来构造细胞-细胞间相似性矩阵，然后基于构建的相似性矩阵采用传统谱聚类的方法来对细胞进行有效分类。Park et al.等学者基于SIMLR提出了一种改进的方法MPSSC来对单细胞RNA测序数据进行聚类，首先加入双随机相似性矩阵来得到相似性矩阵，然后通过对目标矩阵强加稀疏性约束来得到细胞分类。基于矩阵分解方面，Shao et al.等人采用非负矩阵分解(nonnegative matrix factorization(NMF))方法对单细胞RNA测序数据进行分解，得到基矩阵和系数矩阵，然后应用基矩阵来对单细胞进行分类，采用系数矩阵来识别特征基因。为了进一步提高聚类效果，Elyanow et al.,等人结合基因相互作用关系先验信息，提出了一种新的融合已知基因相互作用网络的非负矩阵分解方法NMFsc来将原始矩阵分解为两个低维矩阵：基因矩阵和细胞矩阵，认为已知网络中两基因如果相互连接那么在低维的基因矩阵中也有相似的表示。该方法不仅可以通过将降维的低阶细胞矩阵进行聚类，也可以用于原始数据矩阵中缺失元素的填充。然而，该方法的聚类效果依赖于高质量的已知基因网络。由于单细胞RNA测序数据存在严重的数据缺失(dropout)，这给传统降维、聚类的方法带来了巨大挑战，为了克服这个问题，一些学者提出了一些列数据填充的方法，如MAGIC，sclmpute,Mclmpute等来尝试通过各种方法先补全缺失数据，然后基于补全的数据再进行聚类。近来Li et al.等学者提出一个基于网络嵌入(network embedding)的表示学习方法SCRL，通过数据驱动的非线性映射和引进先验信息(比如基因相互作用信息，pathway information)来对细胞和基因学习一个更有意义的低维表示矩阵，然后基于这个低维表示矩阵来进行聚类。Zheng et al.等学者提出了一种新的基于非负低秩表示的方法SinNLRR来获取细胞-细胞间相似性矩阵，然后通过对相似性矩阵进行谱聚类的方法对单细胞进行有效分类，识别新的细胞类型。

虽然上述方法在对单细胞RNA测序数据聚类方面做了一些工作，但是目前针对单细胞RNA测序数据进行聚类的精度和效率仍有待提高。

发明内容

本发明旨在提供一种新的基于单细胞RNA测序数据识别细胞类型的方法，基于矩阵低秩表示和非负约束的方法来对带噪声高维稀疏单细胞RNA测序数据进行聚类，从而识别细胞类型。

本发明提供的一种基于单细胞RNA测序数据识别细胞类型的方法，包括以下步骤：

(1)基于单细胞RNA测序数据，定义一个基于低秩表示的同时刻画数据全局结构和局部关联特征结构的优化模型；

(2)建立增广的拉格朗日函数，采用交替方向乘子法(ADMM)来优化目标函数值，获得迭代优化后的细胞-细胞间相似性矩阵S；

(3)采用谱聚类方法对所述的相似性矩阵S进行聚类，从而得到对应的细胞类型划分。

优选的，步骤(1)所述的基于低秩表示的同时刻画数据全局结构和局部关联特征结构的优化模型具体为：

其中，

X表示含n个细胞m个基因的RNA测序数据；

E表示误差项，用于刻画数据噪声；

Z表示表征矩阵，用于度量细胞-细胞间相似性；

||Z||_*表示矩阵Z的核范数；

用于度量细胞i与细胞j间的相当距离关系，其中α为0-1之间参数。

优选的，步骤(2)所述的增广的拉格朗日函数具体为：

其中，

S和U表示引入的变量，初始值都设为Z；

C₁、C₂和C₃表示拉格朗日乘子；

λ₁、λ₂和μ表示正则化参数，防止过拟合；

<·>表示矩阵内积，例如<A,B>＝tr(A^TB)；

||·||_F表示Frobenius范数。

优选的，步骤(2)所述的采用交替方向乘子法(ADMM)来优化目标函数值包括以下具体步骤：

(2.1)通过如下公式更新Z

(2.2)通过求解如下问题来更新S

其中

d_i、r_i分别为矩阵D、矩阵R的第i行，s_i表示矩阵S的第i行；

(2.3)通过如下公式更新U

其中Θ表示奇异值阈值计算(singular valuethresholding operator)；

(2.4)通过如下公式更新E

其中υ表示收缩运算(shrinkage operator)；

(2.5)分别通过下列公式更新C1、C2、C3、u

C₁＝C₁+μ(X-XZ^k-E^k)

C₂＝C₂+μ(Z^k-S^k)

C₃＝C₃+μ(Z^k-U^k)

μ＝min(μρ,μ_max),其中ρ,μ_max皆为给定常数；

进行迭代，迭代误差计算如下：

Error_value＝max(leq1,max(abs(L3(:)))；

leq1＝max(max(abs(L1(:))),max(abs(L2(:))))；

L1＝X-XZ^k-E^k；L2＝Z^k-S^k；L3＝Z^k-U^k

当迭代次数满足设定的最大迭代步数(一般设定为100)或者误差值Error_value小于设定1e-5时终止迭代，得到优化后的Z；然后按照公式S＝(|Z|+|Z^T|)/2计算得到迭代优化后的细胞-细胞间相似性矩阵S。

优选的，步骤(3)所述的采用谱聚类方法对所述的相似性矩阵S进行聚类包括以下具体步骤：

(3.1)构建规范化拉普拉斯矩阵L＝D^-1/2SD^-1/2，其中D为对角矩阵

(3.2)计算矩阵L的特征值，得到k个最小特征值对应的特征向量V＝[v₁,v₂,…v_k]；

(3.3)采用L2范数对V进行规范化

(3.4)采用Kmeans方法对规范化的矩阵V进行聚类，得到k个类，每个类对应一种细胞类型。

本发明的技术效果是：与现有的技术相比，本发明的方法在考虑数据的全局结构特征基础上，有效的挖掘细胞与细胞间内在关联特性，该方法能有效提高聚类效果；本发明的方法在聚类效果上与已有方法相比具有明显优势，能有效的对实测的单细胞RNA测序数据进行聚类，得到较准确的细胞类型的划分，有助于识别肿瘤单细胞RNA测序数据中的新的细胞类型，为阐明肿瘤细胞异质性提供可靠依据。

附图说明

图1为本发明方法(SCCLRR)与其它四种聚类方法(SinNLRR、tsne、SIMLR、MPSS)在测试数据Kold下聚类结果可视化图的比较，图中点表示细胞，聚为同一类的细胞用一种记号标注。

图2为本发明方法(SCCLRR)与其它四种聚类方法(SinNLRR、tsne、SIMLR、MPSS)在测试数据Ting下聚类结果可视化图的比较，图中点表示细胞，聚为同一类的细胞用一种记号标注。

具体实施方式

下面将结合实例详细说明本发明所具有的有益效果，旨在帮助阅读者更好地理解本发明的实质，但不能对本发明的实施和保护范围构成任何限定。

本发明公提供的一种基于单细胞RNA测序数据识别细胞类型的方法，基于矩阵低秩表示模型和图正则化约束的方法来对带噪声高维稀疏单细胞RNA测序数据进行聚类，有效挖掘单细胞RNA测序数据的全局结构特征和局部关联特性，提出了新的预测关键蛋白质的计算方法。该方法的主要步骤包括：

(1)基于单细胞RNA测序数据X，将细胞与细胞间相似性矩阵构建转化为优化问题，得到该优化问题的数学模型；

(2)基于所构建的优化问题模型，采用交替方向乘子法(ADMM)对问题进行优化求解，选定一个变量，令其它变量固定，单独优化选定的这个变量，然后通过迭代更新的方法来进行优化求解；基于迭代优化后的Z，按照公式S＝(|Z|+|Z^T|)/2计算得到最终的细胞-细胞间相似性矩阵S；

(3)由上述步骤得到的相似性矩阵S及已知的类的个数，采用谱聚类方法基于相似性矩阵S进行聚类，得到对单细胞的分类。

为了测试本发明方法的有效性，采用文献中提供的单细胞RNA测序数据，具体七组单细胞RNA测序数据说明及数据来源如下：

Kold数据从ArrayExpress database数据库中序列号为E-MTAB-2600下载得到，为老鼠胚胎干细胞多能状态在三种条件下的生物学实验数据，过滤掉全为0和缺失的数据后，得到一个包含3个类型涵盖704个细胞10685个基因的单细胞RNA测序数据。

Darmanis数据从GEO(Gene Expression Omnibus)数据库中编号为GSE67835下载得到，去掉混合细胞后得到一个包含8个细胞类型中420个细胞22085个基因的单细胞RNA测序数据。

Treutlein数据从GEO数据库中编号为GSE52583下载得到，过滤掉全部为0的数据后，单细胞RNA测序包含80个肺上皮细胞和5种细胞类型，涉及9352个基因。

Ting数据从GEO数据库中编号为GSE51372下载得到，单细胞RNA测序包括114个细胞14405个基因，共5种细胞类型。

Deng数据从GEO数据库中下载得到，编号为GSE45719，小鼠胚胎单细胞在不同植入期的转录数据，处理后的单细胞RNA测序数据包含7种细胞类型，135个细胞，12548个基因。

Yan_human数据从GEO数据库中下载得到，编号为GSE36552，对人类植入前胚胎和胚胎干细胞的单细胞测序数据，处理后的单细胞RNA测序数据包含8种细胞类型，124个细胞，3840个基因。

Li_islet数据从GEO数据库中下载得到，编号为GSE73727，去掉12个未定义的单细胞及基因的少于20的细胞后，得到包含6种细胞类型的60个单细胞，4494个基因的单细胞RNA测序数据。

首先将这些收集的单细胞RNA测序数据，整理成矩阵格式X，X的行对应细胞，列对应基因。对单细胞RNA测序数据X进行聚类，具体步骤如下：

步骤1.建立联合低秩表示模型和图正则化的数学模型如下：

其中，X表示含n个细胞m个基因的RNA测序数据；E表示误差项，用于刻画数据噪声；Z为表征矩阵，λ1和λ2分别表示正则化参数；||Z||_*表示矩阵Z的核范数；F(x_i,x_j)度量细胞i与细胞j间的相当距离关系，定义如下：

其中第一项为两细胞间规范化的欧拉矩阵，第二项为两细胞间余弦相似性，α为0-1之间参数。

步骤2.对上述所建立的优化模型采用交替方向乘子法(ADMM)进行优化求解。首先基于所建立模型引入变量S,U得到的增广的拉格朗日函数如下：

其中，C1、C2、C3为拉格朗日乘子；λ₁、λ₂和μ为正则化参数；<·>表示矩阵内积，<A,B>＝tr(A^TB)；||·||_F表示Frobenius范数。

基于上述增广拉格朗日函数L，分别选定一个变量，然后固定其它变量情况下对所选变量进行单独优化，得到如下迭代的方法优化步骤(1)中目标函数：

(1).通过如下公式更新Z，第k+1步

(2).通过求解如下问题来更新S

其中

d_i、r_i分别为矩阵D及R的第i行，s_i表示矩阵S的第i行；

(3).通过如下公式更新U

其中Θ表示奇异值阈值计算(singular valuethresholding operator)；

(4).通过如下公式更新E

其中υ表示收缩运算(shrinkage operator)；

(5).更新C₁、C₂、C₃、μ

C₁＝C₁+μ(X-XZ^k-E^k)

C₂＝C₂+μ(Z^k-S^k)

C₃＝C₃+μ(Z^k-U^k)

μ＝min(μρ,μ_max)，其中ρ,μ_max皆为给定常数。

进行迭代，迭代误差计算如下：

Error_value＝max(leq1,max(abs(L3(:)))；

leq1＝max(max(abs(L1(:))),max(abs(L2(:))))；

L1＝X-XZ^k-E^k；L2＝Z^k-S^k；L3＝Z^k-U^k；

当迭代次数满足最大迭代步数100或者误差值Error_value小于1e-5时终止迭代，得到优化后的Z。基于迭代优化后的Z，按照公式S＝(|Z|+|Z^T|)/2计算得到最终的细胞-细胞间相似性矩阵S。

步骤3.基于上述步骤得到的相似性矩阵S，采用谱聚类方法对单细胞数据进行有效聚类，具体操作步骤如下：

(1)构建规范化拉普拉斯矩阵L＝D^-1/2SD^-1/2，其中D为对角矩阵

(2)计算矩阵L的特征值，得到k个最小特征值对应的特征向量V＝[v₁,v₂,…v_k]；

(3)采用L2范数对V进行规范化

(4)采用Kmeans方法对规范化的矩阵V进行聚类，得到k个类，每个类对应一种细胞类型。

本发明的方法可以归结为三个步骤，基于单细胞RNA测序数据矩阵，建立联合低秩表示和自适应图正则化的数学模型，然后采用交替方向乘子法，基于已建立优化模型，构建增广拉格朗日函数L分别选定一个变量然后固定其它变量情况下对所选变量进行单独优化，得到求解优化问题的迭代更新算法，求解出表示矩阵Z，得到细胞与细胞间相似性矩阵S，最后采用谱聚类方法对相似性矩阵S进行聚类，得到对单细胞的分类。

通常，为了评估聚类方法的好坏，采用标准化互信息(NMI)和兰德指数(ARI)来度量，两个度量越大意味着该方法聚类效果越好。

假设真实的聚类标签T，和预测的聚类标签Y，NMI的定义如下：

其中，MI(T,Y)是聚类标签T和Y之间的互信息，H(Y)、H(T)分别是聚类标签Y和T对应的熵，P(t,y)表示t和y的联合概率分布，p(t)和p(y)分别表示t和y的边缘概率。

ARI的定义如下：

其中a_ty表示在真实类和预测的类中都属于同一类的点对数，a_t表示在真实类种属于同一类但是在预测聚类种不属于同一类的点对数，a_y表示在预测的类中属于同一类，但是在真实类中不属于同一类的点对数，a表示在真实类和预测的类中都不属于同一类的点对数。

为了评估本发明方法的有效性，将本发明方法与文献中最新的代表性方法Kmeans、Spectral、tsne、SIMLR、Corr、MPSSC、SinNLRR{Alsabti K.An efficient k-meansclustering algorithm[J].1997；Cristianini N.Spectral kernel methods forclustering[C]//Advances in neural information processing systems.2002:649-655；Van der Maaten L.Visualizing data using t-SNE[J].Journal of MachineLearning Research,2008,9(2579-2605):85；Wang B.Visualization and analysis ofsingle-cell RNA-seq data by kernel-based similarity learning[J].Naturemethods,2017,14(4):414；Jiang H.Single cell clustering based on cell-pairdifferentiability correlation and variance analysis[J].Bioinformatics,2018,34(21):3684-3694；Park S.Spectral clustering based on learning similarity matrix[J].Bioinformatics,2018,34(12):2069-2076；Zheng R.SinNLRR:a robust subspaceclustering method for cell type detection by non-negative and low-rankrepresentation[J].Bioinformatics,2019}在已有的单细胞RNA测序数据下进行比较，采用NMI和ARI两种度量方法对聚类效果进行评估，比较各种方法的优劣。

A.基于聚类结果NMI指标的比较

表1结果显示，本发明方法与其它六种聚类方法7种测试的单细胞RNA测序数据下聚类结果NMI指标的比较，本发明方法在所有数据下的聚类结果均优于其它几种聚类方法，其中，在Kold,Ting，Li_islet数据下聚类结果与真实细胞类型完全一致。

表1.新发明方法与其它几种聚类方法在测试单细胞RNA测序数据下聚类结果NMI的比较。

B.基于ARI指标的比较

表2为本发明方法在七组测试的单细胞RNA测序数据下聚类结果的ARI指标与其它几种方法聚类结果的比较，ARI越大意味着该方法聚类效果越好，从表2中可以看出本发明方法聚类效果明显优于其它方法。

表2.新发明方法与其它几种聚类方法在测试单细胞RNA测序数据下聚类结果ARI的比较。

由此可知，本发明方法具有很好的普适性，可以适用于不同物种下单细胞RNA测序数据的聚类，适用范围广、精确度高。

以上所述的仅仅是对本发明的实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (4)

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

(1)基于单细胞RNA测序数据，定义一个基于低秩表示的同时刻画数据全局结构和局部关联特征结构的优化模型；

其中，

X表示含n个细胞m个基因的RNA测序数据；

E表示误差项，用于刻画数据噪声；

Z表示表征矩阵，用于度量细胞-细胞间相似性；

||Z||_*表示矩阵Z的核范数；

用于度量细胞i与细胞j间的相当距离关系，其中α为0-1之间参数；

(2)基于优化模型建立增广的拉格朗日函数，采用交替方向乘子法ADMM来优化目标函数值，获得迭代优化后的细胞-细胞间相似性矩阵S；

(3)采用谱聚类方法对所述的相似性矩阵S进行聚类，从而得到对应的细胞类型划分。

2.根据权利要求1所述的方法，其特征在于：步骤(2)所述的增广的拉格朗日函数具体为：

其中，

S和U表示引入的变量，初始值都为Z；

C₁、C₂和C₃表示拉格朗日乘子；

λ₁、λ₂和μ表示正则化参数，防止过拟合；

<·>表示矩阵内积；

||·||_F表示Frobenius范数。

3.根据权利要求2所述的方法，其特征在于：步骤(2)所述的采用交替方向乘子法ADMM来优化目标函数值包括以下具体步骤：

(2.1)通过如下公式更新Z

(2.2)通过求解如下问题来更新S

其中

d_i、r_i分别为矩阵D、矩阵R的第i行，s_i表示矩阵S的第i行；

(2.3)通过如下公式更新U

其中Θ表示奇异值阈值计算singular value thresholdingoperator；

(2.4)通过如下公式更新E

其中υ表示收缩运算shrinkage operator；

(2.5)分别通过下列公式更新C1、C2、C3、μ

C₁＝C₁+μ(X-XZ^k-E^k)

C₂＝C₂+μ(Z^k-S^k)

C₃＝C₃+μ(Z^k-U^k)

μ＝min(μρ,μ_max)，其中ρ,μ_max皆为给定常数；

进行迭代，迭代误差计算如下：

Error_value＝max(leq1,max(abs(L3(:)))；

leq1＝max(max(abs(L1(:))),max(abs(L2(:))))；

L1＝X-XZ^k-E^k；L2＝Z^k-S^k；L3＝Z^k-U^k；

当迭代次数满足设定的最大迭代步数设定为100或者误差值Error_value小于设定1e-5时终止迭代,得到优化后的Z；然后按照公式S＝(|Z|+|Z^T|)/2计算得到迭代优化后的细胞-细胞间相似性矩阵S。

4.根据权利要求1所述的方法，其特征在于：步骤(3)所述的采用谱聚类方法对所述的相似性矩阵S进行聚类包括以下具体步骤：

(3.1)构建规范化拉普拉斯矩阵L＝D^-1/2SD^-1/2，其中D为对角矩阵

(3.2)计算矩阵L的特征值，得到k个最小特征值对应的特征向量V＝[v₁,v₂,…v_k]；

(3.3)采用L2范数对V进行规范化

(3.4)采用Kmeans方法对规范化的矩阵V进行聚类，得到k个类，每个类对应一种细胞类型。

Images