CN115223659A - 一种基于低秩表征和改进谱聚类的单细胞rna测序数据聚类方法 - Google Patents

Abstract

本申请涉及一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，它包括如下步骤：基于单细胞RNA测序数据矩阵，建立联合低秩表示和自适应图正则化的数学模型；采用交替方向乘子法，构建增广拉格朗日函数，对每个变量进行单独优化，得到细胞与细胞间相似性矩阵S；采用谱聚类方法对相似性矩阵S进行聚类，得到对单细胞的分类。本发明可有效的将低秩表征来刻画数据全局结构，结合谱聚类模型将数据非线性关联特征自适应的结合起来建立优化模型，并且实现对单细胞的细胞类型分类；实验数值实验结果表明，本发明可以对单细胞RNA测序数据进行有效聚类。

Description

一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类 方法

技术领域

本发明涉及生物信息学领域，具体涉及一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法。

背景技术

传统高通量测序技术也叫批量测序，得到的是组织内细胞群体(包含免疫细胞、成纤维细胞、肿瘤细胞和巨噬细胞等的混合物)基因表达的平均值，很难鉴别细胞个体之间表达的差异性。

近些年来，随着高通量生物学技术的快速发展，单细胞测序技术能得到单个细胞内大量基因的表达信息，为揭示生物组织中各种细胞类型的转录组特征和阐明细胞之间基因表达的异质性提供了非常有力的工具。与传统的批量化全基因组测序技术相比，单细胞测序技术不仅可以测量单个细胞内基因表达水平，而且还能检测到微量的基因表达子或罕见非编码RNA。

设计有效的计算方法对单细胞RNA测序数据进行聚类，是后续挖掘标志基因、构建细胞伪轨迹，阐明细胞异质性的内在分子机制的前提条件。对单细胞进行准确聚类，有助于鉴定肿瘤细胞亚型，对肿瘤的精确诊断和个性化医疗以及药物分子靶标设计具有非常重要的意义。

由于目前单细胞测序技术的局限，单细胞测序数据具有高噪声、稀疏、较高的dropout(很多低表达或中度表达的基因无法有效检测到，导致数据中大部分都是0)等特点，这给单细胞测序数据的分析带来了非常大的挑战。传统数据分析方法难以有效处理单细胞测序数据，近些年来针对单细胞测序数据的聚类分析方法研究已经成为热门的研究领域。目前已有大量的方法被提出，例如，Kiselev等学者提出了一种单细胞一致性聚类方法SC3，该方法通过将原数据空间进行PCA变换构造新的特征空间，然后基于kmeans算法得到在不同参数下的聚类结果，最后将不同聚类结果进行融合得到最终的细胞类型划分。Xu etal.等学者提出SNN-Cliq方法，该方法采用共享邻近邻居SNN方法来构建细胞间相似性矩阵，然后基于复杂网络模块识别方法实现对单细胞进行有效聚类。Wang et al.等学者提出一种名为SIMLR的单细胞聚类方法，该方法通过融合多个高斯核函数，构建低秩结构模型来实现得到细胞-细胞间相似性矩阵，然后采用谱聚类方法对单细胞进行聚类。由于非负矩阵分解(nonnegative matrix factorization(NMF))不仅具有较好的可解释性，同时具备非常好的扩展性，在数据降维以及聚类中有着非常重要的应用。Shao et al.等学者采用非负矩阵分解方法将单细胞RNA测序数据矩阵进行分解，得到基矩阵和系数矩阵，然后应用基矩阵来对单细胞进行分类。为了进一步提高聚类效果，Elyanow et al.等人提出了一种融合已有基因相互左右关系先验信息的矩阵分解方法NMFsc，该方法虽然能对单细胞进行较好的聚类，但是依赖于已知基因关联信息的可靠性。

随着无监督聚类方法在广泛兴起，Zheng et al.等学者提出了一种新的基于非负低秩表示(Low rank representation)的无监督聚类方法SinNLRR来对单细胞RNA测序数据进行聚类，首先通过建立基于非负低秩模型来获取细胞-细胞间相似性矩阵，然后基于相似性矩阵采用谱聚类的方法对单细胞进行有效分类，识别细胞类型。然而，该方法仅考虑数据的全局结构特征，忽视了数据的局部信息特征。Liang et al.等学者将细胞间关系进行建模，得到一个细胞-细胞间相似性的稀疏表征，然后融入三种经典的相似性度量(pearson相关系数，spearman相关系数，余弦相似性)计算细胞间相似性来增强相似性矩阵构建，从而达到对单细胞进行较准确的聚类。

虽然上述方法在单细胞RNA测序数据聚类方面取得了一定的进展，但是基于单细胞RNA测序数据，设计快速、有效、可靠的聚类方法仍需要深入、系统的研究。

发明内容

本发明的目的在于，提供一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，结合矩阵低秩表征理论和对称非负矩阵分解对带噪声高维稀疏单细胞RNA测序数据进行聚类，实现对单细胞进行聚类。

本发明的采取的技术方案是：一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，包括如下步骤：

S1：基于单细胞RNA测序数据，定义联合低秩模型和谱聚类的优化模型，所述优化模型能有效刻画数据的全局特性和局部关联性质，具体描述为：

其中，X表示含n个细胞m个基因的RNA测序数据；E表示误差项，用于刻画数据噪声；Z表示表征矩阵；S表示细胞-细胞间相似性矩阵；||Z||*表示矩阵Z的核范数；λ₁、λ₂和λ₃表示正则化参数；|E||_2,1表示矩阵E的L₂₁范数；

表示F范数的平方；h_i，h_j分别表示表示H矩阵第i列和第j列；z_ij表示矩阵Z的(i,j)位置元素值；

表示L2范数的平方；H表示分解的正交矩阵；Z^T表示矩阵Z的转置；H^T表示矩阵H的转置；

S2：基于优化模型建立增广的拉格朗日函数，采用交替方向乘子法，即ADMM来单独优化每个变量，建立迭代格式，通过迭代循环得到优化后的细胞-细胞间相似性矩阵S；

S3：采用谱聚类方法对所述相似性矩阵S进行聚类，得到对单细胞的聚类结果。

进一步地，所述增广的拉格朗日函数具体描述为：

其中，G和U表示引入的变量，初始值都为Z；G_ij表示矩阵G的(i，j)位置元素值；C₁、C₂和C₃表示拉格朗日乘子；λ₁、λ₂、λ₃和μ表示正则化参数，防止过拟合；<·>表示矩阵内积，例如<A,B>＝tr(A^TB)；||·||_F表示Frobenius范数。

进一步地，所述步骤S2的具体步骤为：

S201：通过如下公式更新Z：

其中，I表示单位矩阵；k表示迭代步数；

S202：通过如下公式更新G：

其中，Dist矩阵中元素(i，j)的值为

S203：通过如下公式更新H：

其中，L_G表示矩阵G对应的拉普拉斯矩阵，

S204：通过如下公式更新U：

其中，Θ表示奇异值阈值计算；

S205：通过如下公式更新E：

其中，υ表示收缩运算；

S206：通过如下公式更新C1、C2、C3和μ：

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

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

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

μ＝min(μρ,μ_max)

其中，ρ,μ_max皆为给定常数；

S207：根据S201～S206更新的参数值进行迭代，计算迭代误差，具体计算公式如下：

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

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

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

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

进一步地，所述步骤S3的具体步骤为：

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

S302：计算规范化拉普拉斯矩阵L的特征值，得到c个最小特征值对应的特征向量V＝[v₁,v₂,…v_c]；

S303：对特征向量V按照下述公式进行规范化，得到规范化后的特征向量V1：

S304：采用Kmeans方法对规范化后的特征向量V1进行聚类，得到c个类，每个类对应一种细胞类型。

进一步地，所述最大迭代步数为100，误差值Error_value的设定值为1e-5。

本发明的有益效果在于：

(1)与现有的技术相比，本发明通过低秩表征来刻画数据的全局结构特征，通过自适应的正则化项结合对称非负矩阵分解来有效挖掘数据的局部嵌入特性；数值实验比较结果表明，本发明能有效提高单细胞聚类精度；

(2)本发明基于低维正交空间H下来度量细胞间距离

可以有效减少噪声的影响，相较于在原始数据空间下度量细胞间“距离”有明显优势；同时本发明通过引入对称分解的形式

有助于刻画数据中非线性关联特性，有效提高精度；

(3)本发明的方法在单细胞聚类精确度上与已有方法相比具有明显优势，能有效的对实测的单细胞RNA测序数据进行聚类，将单细胞按照其表达特征进行准确聚类，有助于鉴定肿瘤细胞亚型，为辨别肿瘤类型，设计个性化医疗方案和药物提供理论依据。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1为本发明实施例的方法流程图。

图2为本发明实施例所述的方法(LRRS)与其它5种聚类方法(SinNLRR、SoptSC、SIMLR、MPSSC、Corr)在测试数据Ting下聚类结果可视化图的比较图；图中点表示细胞，聚为同一类的细胞用一种记号标注。

图3为本发明实施例所述的方法(LRRS)与其它5种聚类方法(SinNLRR、SoptSC、SIMLR、MPSSC、Corr)在测试数据Zheng下聚类结果可视化图的比较图；图中点表示细胞，聚为同一类的细胞用一种记号标注。

具体实施方式

为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明并不限于下面公开的具体实施例的限制。

除非另作定义，此处使用的技术术语或者科学术语应当为本申请所述领域内具有一般技能的人士所理解的通常意义。本专利申请说明书以及权利要求书中使用的“第一”、“第二”以及类似的词语并不表示任何顺序、数量或者重要性，而只是用来区分不同的组成部分。同样，“一个”或者“一”等类似词语也不表示数量限制，而是表示存在至少一个。“连接”或者“相连”等类似的词语并非限定于物理的或者机械的连接，而是可以包括电性的连接，不管是直接的还是间接的。“上”、“下”、“左”、“右”等仅用于表示相对位置关系，当被描述对象的绝对位置改变后，则该相对位置关系也相应地改变。

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

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

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

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

数据从GEO数据库中下载得到，编号为GSE38495，由平台Smarrt-Seq得到，处理后的单细胞RNA测序数据包含7种细胞类型，33个细胞，3575个基因。

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

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

Zheng数据来源于SRP数据库，编号为SRP073767,预处理后数据种包含500个细胞，4776个基因，细胞类型有3种。

CellLines数据由平台Illumina HiSeq 2000为描述大肠癌异质性的单细胞RNA测序数据，编号为GSE81861，原始数据包含630个细胞，预处理后的数据包含7个类型的561个细胞，涉及57241个基因。

本发明实施例所述的基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，具体步骤如下：

S1：将单细胞RNA测序数据整理成矩阵格式，得到矩阵X，其中X的行对应细胞，列对应基因；基于矩阵X，定义联合低秩模型和谱聚类的优化模型，所述优化模型能有效刻画数据的全局特性和局部关联性质；所述优化模型的具体描述为：

其中，X表示含n个细胞m个基因的RNA测序数据；E表示误差项，用于刻画数据噪声；Z表示表征矩阵；S表示细胞-细胞间相似性矩阵；||Z||*表示矩阵Z的核范数；λ₁、λ₂和λ₃表示正则化参数；|E||_2,1表示矩阵E的L₂₁范数；

表示F范数的平方；h_i，h_j分别表示表示H矩阵第i列和第j列；z_ij表示矩阵Z的(i，j)位置元素值；

表示L2范数的平方；H表示分解的正交矩阵；Z^T表示矩阵Z的转置；H^T表示矩阵H的转置。在本发明实施例中，λ₁、λ₂和λ₃的取值分别为10，10和1。

本发明实施例基于低维正交空间H下来度量细胞间距离

可以有效减少噪声的影响，相较于在原始数据空间下度量细胞间“距离”有明显优势；同时本发明实施例通过引入对称分解的形式

有助于刻画数据中非线性关联特性，有效提高精度。

S2：基于优化模型建立增广的拉格朗日函数，所述增广的拉格朗日函数具体描述为：

其中，G和U表示引入的变量，初始值都为Z；G_ij表示矩阵G的(i，j)位置元素值；C₁、C₂和C₃表示拉格朗日乘子；λ₁、λ₂、λ₃和μ表示正则化参数，防止过拟合；<·>表示矩阵内积，例如<A,B>＝tr(A^TB)；||·||_F表示Frobenius范数；

采用交替方向乘子法(ADMM)，固定其它变量，每个变量进行单独优化，建立迭代格式，通过迭代更新的方法来得到模型的解Z，并计算得到优化后的细胞-细胞间相似性矩阵S，具体步骤如下：

S201：通过如下公式更新Z：

其中，I表示单位矩阵；k表示迭代步数；

S202：通过如下公式更新G：

其中，Dist矩阵中元素(i，j)的值为

S203：通过如下公式更新H：

其中，L_G表示矩阵G对应的拉普拉斯矩阵，

S204：通过如下公式更新U：

其中，Θ表示奇异值阈值计算；

S205：通过如下公式更新E：

其中，υ表示收缩运算；

S206：通过如下公式更新C1、C2、C3和μ：

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

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

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

μ＝min(μρ,μ_max)

其中，ρ,μ_max皆为给定常数；

S207：根据S201～S206更新的参数值进行迭代，计算迭代误差，具体计算公式如下：

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

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

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

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

S3：由步骤S2得到的相似性矩阵S和测试数据中提供的已知的类的数目，采用谱聚类方法对所述相似性矩阵S进行聚类，得到对单细胞的聚类结果，具体步骤如下：

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

S302：计算规范化拉普拉斯矩阵L的特征值，得到c个最小特征值对应的特征向量V＝[v1,v2,…v_c]；

S303：对特征向量V按照下述公式进行规范化，得到规范化后的特征向量V1：

S304：采用Kmeans方法对规范化后的特征向量V1进行聚类，得到c个类，每个类对应一种细胞类型。

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

通常，为了评估聚类方法的好坏，采用标准化互信息(NMI)和准确度(ACC)来度量，两个度量都是在0-1范围内，度量值越大意味着聚类方法的聚类效果越好，取值为1时意味着聚类结果与真实分类结果完全一致。

假设真实的聚类标签为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的边缘概率。

ACC的定义如下：

其中，t_i和y_i分别表示第i个细胞真实的类标签和预测的类标签，n为单细胞总数目，map(yi)表示最佳类标的映射函数，该函数一般通过匈牙利算法(Kuhn-Munkres orHungarian Algorithm)实现，从而在多项式时间内求解该标签最佳分配问题。

为了评估本发明方法的有效性，将本发明方法与已有文献中代表性方法(Kmeans(Alsabti K,Ranka S,Singh V.An efficient k-means clustering algorithm[J].1997.)、Spectral(Cristianini N,Shawe-Taylor J,Kandola J S.Spectral kernelmethods for clustering[C]//Advances in neural information processingsystems.2002:649-655.)、tsne(Van der Maaten L,Hinton G.Visualizing data usingt-SNE[J].Journal of Machine Learning Research,2008,9(2579-2605):85.)、SIMLR(Wang B,Zhu J,Pierson E,et al.Visualization and analysis of single-cell RNA-seq data by kernel-based similarity learning[J].Nature methods,2017,14(4):414.)、Corr(Jiang H,Sohn L L,Huang H,et al.Single cell clustering based oncell-pair differentiability correlation and variance analysis[J].Bioinformatics,2018,34(21):3684-3694.)、MPSSC(PSSC)(Park S,Zhao H.Spectralclustering based on learning similarity matrix[J].Bioinformatics,2018,34(12):2069-2076.)、SinNLRR(Zheng R,Li M,Liang Z,et al.SinNLRR:a robust subspaceclustering method for cell type detection by non-negative and low-rankrepresentation[J].Bioinformatics,2019.)、Corr(H.Jiang,L.L.Sohn,et al.Singlecell clustering based on cell-pair differentiability correlation and varianceanalysis,Bioinformatics,2018,34(21):3684-3694.)、SoptSC(Wang S,MacLean A L,NieQ.SoptSC:similarity matrix optimization for clustering,lineage,and signalinginference[J].bioRxiv,2018:168922.))在常用的单细胞RNA测序数据下进行比较，采用NMI和ACC两种度量方法对聚类效果进行评估，比较各种方法的优劣。

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

表1结果显示，本发明实施例所述的方法与其它9种聚类方法在7组测试的单细胞RNA测序数据下聚类结果NMI指标的比较，本发明方法在所有数据下的聚类结果均优于其它几种聚类方法，其中，在

Li_islet，Zheng数据下聚类结果与真实细胞类型完全一致。

表1本发明实施例与其它聚类方法在测试单细胞RNA测序数据下聚类结果NMI的比较

B.基于ACC指标的比较

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

表2本发明实施例与其它聚类方法在测试单细胞RNA测序数据下聚类结果ACC的比较

与现有的技术相比，本发明实施例通过低秩表征来刻画数据的全局结构特征，通过自适应的正则化项结合对称非负矩阵分解来有效挖掘数据的局部嵌入特性；数值实验比较结果表明，本发明实施例能有效提高单细胞聚类精度，与已有方法相比具有明显优势，能有效的对实测的单细胞RNA测序数据进行聚类，将单细胞按照其表达特征进行准确聚类，有助于鉴定肿瘤细胞亚型，为辨别肿瘤类型，设计个性化医疗方案和药物提供理论依据。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，其特征在于，包括如下步骤：

S1：基于单细胞RNA测序数据，定义联合低秩模型和谱聚类的优化模型，所述优化模型能有效刻画数据的全局特性和局部关联性质，具体描述为：

其中，X表示含n个细胞m个基因的RNA测序数据；E表示误差项，用于刻画数据噪声；Z表示表征矩阵；S表示细胞-细胞间相似性矩阵；||Z||*表示矩阵Z的核范数；λ₁、λ₂和λ₃表示正则化参数；|E||_2,1表示矩阵E的L₂₁范数；

表示F范数的平方；h_i，h_j分别表示表示H矩阵第i列和第j列；z_ij表示矩阵Z的(i,j)位置元素值；

表示L2范数的平方；H表示分解的正交矩阵；Z^T表示矩阵Z的转置；H^T表示矩阵H的转置；

S2：基于优化模型建立增广的拉格朗日函数，采用交替方向乘子法，即ADMM来单独优化每个变量，建立迭代格式，通过迭代循环得到优化后的细胞-细胞间相似性矩阵S；

S3：采用谱聚类方法对所述相似性矩阵S进行聚类，得到对单细胞的聚类结果。

2.根据权利要求1所述的一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，其特征在于，所述增广的拉格朗日函数具体描述为：

其中，G和U表示引入的变量，初始值都为Z；G_ij表示矩阵G的(i，j)位置元素值；C₁、C₂和C₃表示拉格朗日乘子；λ₁、λ₂、λ₃和μ表示正则化参数，防止过拟合；<·>表示矩阵内积，例如<A,B>＝tr(A^TB)；||·||_F表示Frobenius范数。

3.根据权利要求2所述的一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，其特征在于，所述步骤S2的具体步骤为：

S201：通过如下公式更新Z：

其中，I表示单位矩阵；k表示迭代步数；

S202：通过如下公式更新G：

其中，Dist矩阵中元素(i，j)的值为

S203：通过如下公式更新H：

其中，L_G表示矩阵G对应的拉普拉斯矩阵，

S204：通过如下公式更新U：

其中，Θ表示奇异值阈值计算；

S205：通过如下公式更新E：

其中，υ表示收缩运算；

S206：通过如下公式更新C1、C2、C3和μ：

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

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

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

μ＝min(μρ,μ_max)

其中，ρ,μ_max皆为给定常数；

S207：根据S201～S206更新的参数值进行迭代，计算迭代误差，具体计算公式如下：

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

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

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

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

4.根据权利要求3所述的一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，其特征在于，所述步骤S3的具体步骤为：

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

S302：计算规范化拉普拉斯矩阵L的特征值，得到c个最小特征值对应的特征向量V＝[v₁,v₂,…v_c]；

S303：对特征向量V按照下述公式进行规范化，得到规范化后的特征向量V1：

S304：采用Kmeans方法对规范化后的特征向量V1进行聚类，得到c个类，每个类对应一种细胞类型。

5.根据权利要求3所述的一种基于低秩表征和改进谱聚类的单细胞RNA测序数据聚类方法，其特征在于，所述最大迭代步数为100，误差值Error_value的设定值为1e-5。

Images