增强scRNA-seq数据基因表达相互作用的方法、设备和介质

摘要

本发明公开了一种增强scRNA‑seq数据基因表达相互作用的方法、设备和介质，属于数据处理技术领域。所述方法首先对细胞‑基因表达谱矩阵进行主成分分析，获得能够表明细胞区别的特征基因，基于特征基因计算细胞距离，再根据细胞距离获得细胞相似性，进一步基于细胞相似性得到马尔可夫转移概率，并对细胞‑基因表达谱矩阵进行多重插补。本发明的方法能够消除细胞‑基因表达谱矩阵中的表达量噪声并填补缺失的表达量，最终能够有效增强基因表达相互作用，进一步可以提高细胞的聚类效果，有效进行细胞分型。

权利要求

1.一种增强scRNA‑seq数据基因表达相互作用的方法，其特征在于，包括以下步骤：S1，获得单细胞转录组测序数据的细胞‑基因表达谱矩阵A；

S2，基于所述细胞‑基因表达谱矩阵A，利用主成分分析方法筛选N个主成分，其中N=20~

50，得到PCA矩阵；

S3，基于所述N个主成分的值计算任意两个细胞之间的距离 ，得到细胞距离矩阵D；

S4，基于所述细胞距离矩阵D利用核函数计算任意两个细胞之间的相似性 ，得到细胞相似性矩阵C；

S5，基于细胞相似性矩阵C，根据以下算法，计算任意两个细胞之间的转移概率，得到细胞转移概率矩阵 ：所述两个细胞包括第一细胞和第二细胞，所述第一细胞与所述第二细胞的转移概率为所述第一细胞和所述第二细胞的相似性与第一细胞与所有细胞的相似性的总和的比值；

S6，基于所述细胞转移概率矩阵 ，根据以下算法，对所述细胞‑基因表达谱矩阵A进行插补：将第 个细胞第 个基因的表达量转换其余各个细胞的第 个基因的表达量与所述第个细胞与相应细胞的转移概率的乘积的总和，其中， =1 且 =1 ， 代表所述细胞‑基因表达谱矩阵A中细胞数目，=1 ，~ ~ ~代表所述细胞‑基因表达谱矩阵A中基因数目，通过对转移概率矩阵 求幂，对所述细胞‑基因表达谱矩阵A进行多重插补：其中， 代表第 次插补后的细胞‑基因表达谱矩阵；代表迭代次数，从1开始并逐次加1进行多重插补；当进行第 次插补时，对于每个细胞的每个基因，转换所采用的表达量均为所述细胞‑基因表达谱矩阵A中各基因的表达量，利用以下方式确定迭代次数：

进行第 次插补后，先计算 中第 个细胞第 个基因的相对表达量：再计算第 个细胞插补后的平均基因表达量：

再基于第 次插补后和插补后的全部基因相对表达量的残差平方和与插补后的全部基因相对表达量的离差平方和，计算细胞 的决定系数 ，公式如下：若所有细胞的决定系数均小于预设阈值P1，则停止迭代，得到多重插补后的细胞‑基因表达谱矩阵 ，其中，=1 ，P1=3% 10%。

~ ~

2.根据权利要求1所述的一种增强scRNA‑seq数据基因表达相互作用的方法，其特征在于，步骤S3中，所述距离为欧氏距离，所述 的计算公式如下：其中， 代表第i个PC在第 个细胞中的值， 代表第i个PC在第 个细胞中的值；=1 N。

~

3. 根据权利要求1或2所述的一种增强scRNA‑seq数据基因表达相互作用的方法，其特征在于，步骤S4中，在计算相似性之前，进一步使用UMAP算法对PCA矩阵进行非线性降维，降维后的结果为二维坐标信息，细胞基于表型的相似性聚集到不同区域，对于处于不同区域的细胞，两两之间的相似性设置为0，对于同一区域的细胞，再利用核函数计算任意两个细胞之间的相似性，所述核函数为高斯核函数，所述 的计算公式如下其中， 是带宽，用于控制径向作用范围， =1 30。

~

4.根据权利要求3所述的一种增强scRNA‑seq数据基因表达相互作用的方法，其特征在于，根据以下方式确定 取值：先确定细胞 和细胞 所在区域的密度，当密度小于0.3时， 取值为20 30，密度大于~等于0.3且密度小于等于0.6时， 取值为5 20，密度大于0.6时， 取值为1 5。

~ ~

5.根据权利要求3所述的一种增强scRNA‑seq数据基因表达相互作用的方法，其特征在于，得到细胞转移概率矩阵 后，进一步进行如下处理：（1）对于同一区域中任意细胞，确定与该两个细胞距离之和最小的15个细胞作为该两个细胞的邻居细胞，若不足15个，则全部其余细胞为邻居细胞；

（2）对于每个细胞，将其与非邻居细胞的细胞转移概率设置为0。

6.根据权利要求1所述的一种增强scRNA‑seq数据基因表达相互作用的方法，其特征在于，在进行第 次插补时，若 小于预设阈值P2，则先设置 ，再进行插补，其中，P2=0.00005 0.001。

~

7.一种计算机设备，其特征在于，包括：

存储器，用于存储计算机程序；

处理器，用于执行所述计算机程序时实现如权利要求1‑6任一所述的一种增强scRNA‑seq数据基因表达相互作用的方法的步骤。

8.一种计算机可读存储介质，其特征在于，

所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1‑6任一所述的一种增强scRNA‑seq数据基因表达相互作用的方法的步骤。

说明书

增强scRNA‑seq数据基因表达相互作用的方法、设备和介质

技术领域

[0001] 本发明属于数据处理技术领域，具体地，涉及一种增强scRNA‑seq数据基因表达相互作用的方法、设备和介质。

背景技术

[0002] 单细胞RNA测序（Single cell RNA sequencing，scRNA‑seq）又叫单细胞转录组测序，是一种在单细胞水平上利用RNA测序对特定细胞群体进行基因表达谱定量的高通量实验技术，是近年来的一项热门技术。对于多细胞生物来说，细胞与细胞之间是有差异性的（cell heterogeneity），即细胞异质性。这种细胞异质性可以体现为不同的遗传背景，不同的分化状态，不同的物理特征，不同的基因突变谱和转录组、蛋白质组表达谱等。

[0003] 对于单个细胞，由于mRNA分子采样不足，并不是所有的mRNA分子都能被捕获到，并且由于测序深度比较浅，一般每个细胞仅能检测到10%～50%的转录本，这导致细胞中许多基因计数为0，造成scRNA‑seq测序结果中细胞基因表达矩阵的稀疏性。这种细胞基因表达矩阵的稀疏性给后续分析工作增加计算难度，也可能会严重模糊重要的基因间的相互关系。

[0004] 为了克服单细胞RNA测序结果中细胞基因表达矩阵的稀疏性，目前大多数方法通过聚类的方式将成千上万的细胞聚类成少量的簇；或者通过其他方法合并基因（例如主成分分析[PCA]），创建“metagene”。虽然这些方法在一定程度上处理了稀疏性，但失去了单细胞或单基因的分辨率。

发明内容

[0005] 为了解决上述技术问题，发明人旨在提供一种新的处理单细胞测序数据中稀疏性的方法，提高细胞聚类效果。为此，本发明采用的技术方案如下：

[0006] 本发明第一方面提供一种增强scRNA‑seq数据基因表达相互作用的方法，包括以下步骤：

[0007] S1，获得单细胞转录组测序数据的细胞‑基因表达谱矩阵A；

[0008] S2，基于所述细胞‑基因表达谱矩阵A，利用主成分分析方法筛选N个主成分，其中N=20 50，得到PCA矩阵；~

[0009] S3，基于所述N个主成分的值计算任意两个细胞之间的距离 ，得到细胞距离矩阵D；

[0010] S4，基于所述细胞距离矩阵D利用核函数计算任意两个细胞之间的相似性，得到细胞相似性矩阵C；

[0011] S5，基于细胞相似性矩阵C，根据以下算法，计算任意两个细胞之间的转移概率，得到细胞转移概率矩阵 ：

[0012] 所述两个细胞包括第一细胞和第二细胞，所述第一细胞与所述第二细胞的转移概率为所述第一细胞和所述第二细胞的相似性与第一细胞与所有细胞的相似性的总和的比值；

[0013] S6，基于所述细胞转移概率矩阵 ，根据以下算法，对所述细胞‑基因表达谱矩阵A进行插补：

[0014] 将第 个细胞第 个基因的表达量转换其余各个细胞的第 个基因的表达量与所述第 个细胞与相应细胞的转移概率的乘积的总和，

[0015] 其中， =1 且 =1 ， 代表所述细胞‑基因表达谱矩阵A中细胞数目，=1~ ~ ~， 代表所述细胞‑基因表达谱矩阵A中基因数目。

[0016] 在本发明的一些实施方案中，在步骤S2进行主成分分析之前，进一步包括对所述对细胞‑基因表达谱矩阵A根据文库大小进行归一化的步骤，归一化的方法为：

[0017] ；

[0018] 其中， 代表归一化后第 个细胞中第 个基因的表达量， 代表归一化前第 个细胞中第 个基因的表达量； 代表第 个细胞所有基因的表达量总和，=1 ； 代表所有细胞所有基因表达量总和的平均值。
~

[0019] 主成分分析（Principal Component Analysis，PCA）是一种使用最广泛的数据降维算法。降维就是一种对高维度特征数据预处理方法。降维是将高维度的数据保留下最重要的一些特征，去除噪声和不重要的特征，从而实现提升数据处理速度的目的。

[0020] 在本发明的一些实施方案中，选择能够反应原始细胞差异的前N个主成分（PC）作为后续分析的数据。当N=20 50时，能反应约80%以上的信息。在本发明的一些优选实施方案~中，选择N=30，即选择20个主成分进行分析，结果是足够稳健的。

[0021] 在本发明的一些实施方案中，步骤S3中，所述距离为欧氏距离，所述 的计算公式如下：

[0022] ；

[0023] 其中， 代表第i个PC在第 个细胞中的值， 代表第i个PC在第个细胞中的值；=1 N。
~

[0024] 在本发明的一些实施方案中，步骤S4中，在计算相似性之前，进一步使用UMAP（Uniform Manifold Approximation and Projection）算法对PCA矩阵进行非线性降维，降维后的结果为二维坐标信息，细胞基于表型的相似性聚集到不同区域。在本发明中，处于不同区域的细胞对稀疏矩阵的插补作用不明显，因此，对于处于不同区域的细胞，两两之间的相似性设置为0，对于同一区域的细胞，利用核函数计算任意两个细胞之间的相似性，[0025] 所述核函数为高斯核函数，所述 的计算公式如下：

[0026] ；

[0027] 其中， 是带宽，用于控制径向作用范围， =1 30。~

[0028] 不同的细胞具有不同的基因表达，不同的基因表达反映了细胞可能的表型，表型反映到细胞当中可以是细胞类型，如发育时期状态等。由于不同细胞的表型并不一致，因此针对不同细胞类型选择不同的 值， 太小（低于0.01时）会导致结果不稳定，准确性降低，即相同表型的细胞也难以分类为同类型细胞； 越大高斯核函数的局部影响范围就会越大，过大时（大于100）产生过拟合的情况，即导致不同表型、距离较远的细胞也会被平均到一起，失去数据的分辨率。

[0029] 在本发明的一些实施方案中，根据以下方式确定 取值：

[0030] 先确定细胞 和细胞 所在区域的密度，当密度小于0.3时， 取值为20 30，密度~大于等于0.3且密度小于等于0.6时， 取值为5 20，密度大于0.6时， 取值为1 5。
~ ~

[0031] 在本发明的一些实施方案中，得到细胞转移概率矩阵 后，进一步进行如下处理：

[0032] （1）对于同一区域中任意细胞，确定与该两个细胞距离之和最小的15个细胞作为该两个细胞的邻居细胞，若不足15个，则全部其余细胞为邻居细胞；

[0033] （2）对于每个细胞，将其与非邻居细胞的细胞转移概率设置为0。

[0034] 在本发明的一些实施方案中，步骤S6中，通过对转移概率矩阵 求幂，对所述细胞‑基因表达谱矩阵A进行多重插补：

[0035] ；

[0036] 其中， 代表第 次插补后的细胞‑基因表达谱矩阵；代表迭代次数，从1开始并逐次加1进行多重插补；当进行第 次插补时，对于每个细胞的每个基因，转换所采用的表达量均为所述细胞‑基因表达谱矩阵A中各基因的表达量。

[0037] 在本发明中，由于马尔可夫转移概率矩阵 的特征值均为[0,1]之间，特征值通过指数运算会逐渐减小，其范围也在[0,1]之间。随着马尔可夫转移概率矩阵 多次求幂，除1以外的所有特征值的大小不断减小，如此能够降低噪声的重要性，其解释能力接近于零。随着 的增加，细胞从它们的邻居细胞那里学习缺失的值，迅速获得生物学上非常相似的细胞之间的关系。

[0038] 在本发明的一些实施方案中，利用以下方式确定迭代次数：

[0039] 进行第 次插补后，先计算 中第 个细胞第 个基因的相对表达量：

[0040] ；

[0041] 再计算第 个细胞插补后的平均基因表达量：

[0042] ；

[0043] 再基于第 次插补后和插补后的全部基因相对表达量的残差平方和与插补后的全部基因相对表达量的离差平方和，计算细胞 的决定系数 ，公式如下：

[0044] ；

[0045] ；

[0046] ；

[0047] 若所有细胞的决定系数均小于预设阈值P1，则停止迭代，得到多重插补后的细胞‑基因表达谱矩阵 ，

[0048] 其中，=1 ，P1=3% 10%。~ ~

[0049] 在本发明的一些实施方案中，在进行第 次插补时，若 小于预设阈值P2，则先设置 ，再进行插补，其中，P2=0.00005 0.001。~

[0050] 本发明第二方面提供一种计算机设备，包括：

[0051] 存储器，用于存储计算机程序；

[0052] 处理器，用于执行所述计算机程序时实现如本发明第一方面任一所述的一种增强scRNA‑seq数据基因表达相互作用的方法的步骤。

[0053] 本发明第三方面提供一种计算机可读存储介质，

[0054] 所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如本发明第一方面任一所述的一种增强scRNA‑seq数据基因表达相互作用的方法的步骤。

[0055] 本发明的有益效果

[0056] 相对于现有技术，本发明的技术效果如下：

[0057] 本发明的方法基于马尔可夫随机模型对单细胞转录组测序数据进行多重插补，从而能够消除细胞‑基因表达谱矩阵中的表达量噪声并填补缺失的转录本计数。

[0058] 利用本发明的方法对单细胞转录组测序数据得到的细胞‑基因表达谱矩阵进行多重插补后，能够有效增强基因表达相互作用，进一步可以提高细胞的聚类效果，有效进行细胞分型。

附图说明

[0059] 图1示出了单细胞转录组测序得到的细胞‑基因表达谱矩阵（部分）。

[0060] 图2示出了本发明实施例1中使用UMAP算法对PCA矩阵信息再次进行非线性降维后的聚类图。

[0061] 图3示出了本发明实施例1中人PBMC单细胞转录组测序数据进行加权多重插补前后CD3D基因表达分析结果。

[0062] 图4示出了本发明实施例1中人PBMC单细胞转录组测序数据进行加权多重插补前后MS4A1基因表达分析结果。

[0063] 图5示出了本发明实施例1中人PBMC单细胞转录组测序数据进行加权多重插补前后CD14基因与ITGAM基因相互关系分析结果。

[0064] 图6示出了本发明实施例2中人PBMC单细胞转录组测序数据进行加权多重插补前后CD3D基因表达分析结果。

[0065] 图7示出了本发明实施例2中人PBMC单细胞转录组测序数据进行加权多重插补前后MS4A1基因表达分析结果。

[0066] 图8示出了本发明实施例2中人PBMC单细胞转录组测序数据进行加权多重插补前后CD14基因与ITGAM基因相互关系分析结果。

[0067] 图9示出了本发明实施例3中人PBMC单细胞转录组测序数据进行加权多重插补前后CD3D基因表达分析结果。

[0068] 图10示出了本发明实施例3中人PBMC单细胞转录组测序数据进行加权多重插补前后MS4A1基因表达分析结果。

[0069] 图11示出了本发明实施例3中人PBMC单细胞转录组测序数据进行加权多重插补前后CD14基因与ITGAM基因相互关系分析结果。

具体实施方式

[0070] 除非另有说明、从上下文暗示或属于现有技术的惯例，否则本申请中所有的份数和百分比都基于重量，且所用的测试和表征方法都是与本申请的提交日期同步的。在适用的情况下，本申请中涉及的任何专利、专利申请或公开的内容全部结合于此作为参考，且其等价的同族专利也引入作为参考，特别这些文献所披露的关于本领域中的相关术语的定义。如果现有技术中披露的具体术语的定义与本申请中提供的任何定义不一致，则以本申请中提供的术语定义为准。

[0071] 为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。

[0072] 实施例

[0073] 以下例子在此用于示范本发明的优选实施方案。本领域内的技术人员会明白，下述例子中披露的技术代表发明人发现的可以用于实施本发明的技术，因此可以视为实施本发明的优选方案。但是本领域内的技术人员根据本说明书应该明白，这里所公开的特定实施例可以做很多修改，仍然能得到相同的或者类似的结果，而非背离本发明的精神或范围。

[0074] 除非另有定义，所有在此使用的技术和科学的术语，和本发明所属领域内的技术人员所通常理解的意思相同，在此公开引用及他们引用的材料都将以引用的方式被并入。

[0075] 那些本领域内的技术人员将意识到或者通过常规试验就能了解许多这里所描述的发明的特定实施方案的许多等同技术。这些等同将被包含在权利要求书中。

[0076] 下述实施例中的实验方法，如无特殊说明，均为常规方法。下述实施例中所用的仪器设备，如无特殊说明，均为实验室常规仪器设备；下述实施例中所用的试验材料，如无特殊说明，均为自常规生化试剂商店购买得到的。

[0077] 实施例 1 单细胞测序基因矩阵稀疏性处理

[0078] 1. 获得细胞‑基因表达谱矩阵

[0079] 采集人周血单核细胞（Peripheral Blood Mononuclear Cell，PBMC），进行单细胞测序，获得细胞‑基因表达谱矩阵A，部分如图1和表1所示，图1中，列为不同的细胞（由Barcode表示，Barcode是由16个随机核苷酸组成的），行表示基因在不同细胞中的表达量（转录本count数），“·”代表表达量为0，即该基因转录本count数为0。可见细胞‑基因表达谱矩阵A中很多基因的表达量0。

[0080] 在本实施例中，细胞‑基因表达谱中共包括2887个细胞36477个基因的表达量。

[0081] 表1 细胞‑基因表达谱矩阵

[0082]

[0083] 2. 细胞‑基因表达谱矩阵A归一化及主成分分析

[0084] 为了确保细胞之间的距离反映的是真实生物学差异而不受人为因素影响，发明对原始表达矩阵进行以下两步处理：

[0085] （1）对细胞‑基因表达谱矩阵A根据文库大小进行归一化

[0086] 归一化方式为：

[0087] ；

[0088] 其中， 代表归一化后第 个细胞中第 个基因的表达量（转录本count数）， =1 ， 代表细胞‑基因表达谱矩阵A中细胞数目， 代表归一化前第~
个细胞中第 个基因的表达量； 代表第 个细胞所有基因的表达量总和（转录本count数总和），=1 且 =1 ，n代表细胞‑基因表达谱矩阵A中基因数目；LS代表细胞~ ~
文库大小集合，即每个细胞所有基因表达量总和（转录本count数总和）构成的集合，代表所有细胞转录本count数总和的平均值。

[0089] （2）主成分分析

[0090] 主成分分析（Principal Component Analysis，PCA）是一种使用最广泛的数据降维算法。降维就是一种对高维度特征数据预处理方法。降维是将高维度的数据保留下最重要的一些特征，去除噪声和不重要的特征，从而实现提升数据处理速度的目的。

[0091] 通过对归一化后的数据进行主成分分析，将n个基因即n维数据降为50个主成分（PC），进一步选择能够反应原始细胞差异的前30个主成分（PC）作为后续分析的数据。发明人研究证实，30个主成分能反应约80%以上的信息，因此选择30个主成分进行分析，结果是足够稳健的。

[0092] 对于细胞1和细胞2，30个PC值如表2所示。

[0093] 表2 细胞1和细胞2的PCA矩阵

[0094]

[0095] 部分细胞的前10个PC值如表3所示。

[0096] 表3 部分细胞部分PCA矩阵

[0097]PC_1 PC_2 PC_3 PC_4 PC_5 PC_6 PC_7 PC_8 PC_9 PC_10
细胞1 3.19 ‑1.73 1.68 ‑0.21 3.41 ‑11.20 ‑5.54 6.87 3.36 ‑2.59细胞2 6.90 ‑1.82 ‑4.35 1.11 ‑4.18 ‑3.19 ‑0.25 0.38 ‑0.16 0.30细胞3 ‑0.54 0.15 2.12 ‑0.83 4.14 ‑3.08 ‑3.92 ‑0.45 0.17 ‑0.46细胞4 1.89 ‑1.14 4.43 0.94 1.48 2.20 ‑1.95 ‑0.63 ‑0.89 ‑0.82细胞5 9.00 ‑2.56 ‑5.70 ‑0.60 ‑3.67 2.59 0.72 ‑1.05 ‑0.16 ‑0.65细胞6 4.09 ‑3.68 7.84 3.73 ‑7.94 ‑2.75 ‑5.11 ‑0.06 0.28 0.65细胞7 3.54 ‑0.03 2.20 0.06 1.97 ‑6.59 ‑2.46 1.88 0.80 1.34
细胞8 9.11 ‑2.86 ‑6.11 ‑0.64 ‑5.07 ‑2.04 1.27 1.80 0.74 2.07细胞9 8.04 ‑1.49 ‑4.29 0.50 ‑3.27 2.24 ‑0.68 0.28 ‑0.65 ‑0.63细胞10 3.25 ‑1.82 1.90 ‑0.31 3.74 3.62 ‑0.61 0.93 ‑0.78 ‑0.75细胞11 0.18 ‑0.32 1.95 ‑1.58 4.03 3.06 ‑0.14 0.64 0.14 ‑0.05
细胞12 2.27 ‑3.40 ‑4.47 ‑0.44 ‑3.39 ‑1.42 3.97 1.94 ‑0.61 ‑0.79细胞13 3.96 ‑1.77 ‑1.14 ‑2.80 0.60 ‑3.15 0.36 4.20 ‑0.67 4.84细胞14 4.02 ‑1.11 1.23 ‑0.95 3.31 2.80 ‑0.69 0.23 ‑0.20 ‑1.88细胞15 7.63 ‑1.99 ‑4.59 0.14 ‑2.59 0.61 ‑0.36 ‑1.29 ‑0.07 ‑0.05细胞16 7.31 ‑2.70 ‑2.74 0.23 ‑3.38 ‑3.23 0.72 0.90 0.25 0.34
细胞17 4.76 ‑1.50 ‑2.77 ‑0.68 0.67 ‑3.34 0.91 2.16 ‑0.63 0.80细胞18 3.92 ‑1.70 1.43 ‑1.00 2.93 3.56 ‑1.52 0.73 ‑0.63 ‑0.97细胞19 2.25 ‑1.71 2.00 ‑1.52 3.83 4.17 ‑2.67 ‑0.11 ‑0.68 ‑0.46细胞20 4.82 ‑1.70 ‑1.61 ‑0.93 7.99 ‑12.74 ‑3.87 1.79 2.87 ‑1.61[0098] 3. 计算细胞‑细胞距离矩阵D

[0099] 采用基于细胞的各PC值计算细胞与细胞之间的欧氏距离，计算公式如下：

[0100] ；

[0101] 其中， 代表第 个细胞和第 个细胞之间的欧式距离， =1 且 =1~ ~
， 代表细胞数目， 代表第i个PC在第 个细胞中的值， 代表第i个PC
在第 个细胞中的值；=1 N，N代表用于计算细胞欧式距离的主成分数目，N=30。当p=q时，~
即对于同一细胞，欧式距离值为0。计算可知，上述细胞1和细胞2的欧氏距离为22.95，进而可以计算得到2887个细胞任意两个细胞之间的欧氏距离。

[0102] 部分细胞之间的欧式距离矩阵D如表4所示：

[0103] 表4 细胞欧式距离矩阵

[0104]

[0105] 4. 采用自适应高斯核函数将距离矩阵转换为相似性矩阵

[0106] 使用UMAP算法对PCA矩阵信息再次进行非线性降维，降维后的结果为二维坐标信息，该二维坐标信息在低维坐标上可以呈现出细胞具体的位置，细胞基于表型的相似性聚集到不同区域，如图2所示。对于处于不同区域的细胞，两两之间的相似性设置为0。

[0107] 对于聚集到同一细胞区域的细胞之间的相关性，以往一般采用spearman/person的算法，通过比较细胞中归一化的基因表达量计算两两细胞之间的相关性系数。但是，由于单细胞转录组测序数据具有强稀疏性的特点，即很多基因检测不到转录本，计算得到的相关性均低于0.3，难以区分细胞与细胞之间的相似性强弱。

[0108] 因此，在本实施例中，发明人引入自适应高斯核函数表示聚集到同一区域的两两细胞之间的相关性。

[0109] 高斯核函数也称径向基（RBF）函数，是常用的一种核函数。它可以将有限维数据映射到高维空间，以达到区分两个向量的目的。

[0110] 高斯核函数公式如下：

[0111] ；

[0112] 其中， 代表细胞 和细胞 之间的相似性； 是带宽，控制径向作用范围。

[0113] 由此可见，高斯核函数是两个细胞之间欧式距离的单调函数，两个细胞之间的相关性随两个细胞的欧式距离增加呈现单调递减关系。

[0114] 不同的细胞具有不同的基因表达，不同的基因表达反映了细胞可能的表型，表型反映到细胞当中可以是细胞类型，如发育时期状态等。由于不同细胞的表型并不一致，因此针对不同细胞类型选择不同的 值， 太小（低于0.01时）会导致结果不稳定，准确性降低，即相同表型的细胞也难以分类为同类型细胞； 越大高斯核函数的局部影响范围就会越大，过大时（大于100）产生过拟合的情况，即导致不同表型、距离较远的细胞也会被平均到一起，失去数据的分辨率。

[0115] 由图2可以看出，不同区域细胞（不同细胞类型）聚集的密度不同，如成熟细胞的数目远大于干细胞的数目，因此其对应密度也不一致。密度的范围为[0,1]，密度较大区域的细胞其邻居会比密度小的区域的邻居。为了保证邻居数目尽可能一致，密度大的区域选择更小的 值，密度小的区域选择更大的 值。当密度小于0.3时， 取值为25，密度大于等于0.3且密度小于等于0.6时， 取值为10，密度大于0.6时， 取值为3，通过该方式可以为低维坐标上不同局部区域的细胞确定较为接近的邻居细胞。

[0116] 上述细胞1和细胞2共同聚集于同一区域，该区域的密度为0.4，设置 取值为15，由于获取细胞1和细胞2的高斯核函数为： 0.10。

[0117] 由此，可以得到任意两个细胞之间的相似性，获得相似性矩阵 ，相同细胞相似性为1。部分细胞之间的相似性矩阵 如表5所示：

[0118] 表5 细胞之间的相似性矩阵

[0119]

[0120] 5. 计算马尔可夫转移概率矩阵

[0121] 单细胞转录组测序过程中，由于mRNA的捕获是随机的，因此不同的细胞具有随机mRNA缺失的特点，该过程可以表征为马尔可夫过程。马尔可夫过程具有以下特性：在已知目前状态（现在）的条件下，它未来的演变（将来）不依赖于它以往的演变（过去），每一个状态的转移只依赖于其之前的那一个状态。对于单细胞转录组测序数据来讲，不同的细胞转换为其他细胞的概率是不同的。

[0122] 基于相似性矩阵 可获得马尔可夫转移概率矩阵 。细胞 转换为细胞 的概率（马尔可夫转移概率）是细胞 与细胞 的相似性占细胞 与其他每个细胞相似性总和的比值：

[0123] ；

[0124] 其中， 代表细胞 与细胞 的马尔可夫转移概率， =1 且 =1 ，~ ~
代表细胞 与细胞 的马尔可夫转移概率， 代表第 个细胞（
）， 代表细胞个数； 代表细胞 与其他每个细胞相似
性总和。

[0125] 由此可见，对于马尔可夫转移概率矩阵的每一列或者每一行，其转移概率总和为1。部分细胞之间的马尔可夫转移概率矩阵M如表6所示：

[0126] 表6 细胞之间的马尔可夫转移概率矩阵

[0127]

[0128] 6. 利用马尔可夫转移概率矩阵对基因‑细胞表达谱矩阵进行加权多重插补[0129] 对于单细胞转录组测序数据而言，用于区分细胞表型的主要依据为高变基因（也称为特征向量），而数据的真实结构特征主要由顶部特征向量所体现，其余特征向量可能为噪声。马尔可夫转移概率矩阵中的特征值也反映了该特征信息，因此针对上述步骤首先基于马尔可夫转移概率矩阵 的特征值大小保留顶部特征向量，顶部的特征向量主要来源于低频的特征值。

[0130] 由于数据中噪声一般在数据中呈现较高的频率，因此对于马尔可夫转移概率矩阵进行求幂，可以对马尔可夫转移概率矩阵 中的特征值（转移概率）进行低通滤波，即让低频率通过而滤掉或衰减高频，由此可以过滤掉包含在高频中的噪声。

[0131] 进一步，基于求幂后的马尔可夫转移概率矩阵对细胞‑基因表达谱矩阵A进行加权多重插补：

[0132] ；

[0133] 其中， 代表第 次插补后的细胞‑基因表达谱矩阵；代表迭代次数，也为幂次，从1开始并逐次加1进行多重插补；当进行第 次插补时，对于每个细胞的每个基因，转换所采用的表达量均为所述细胞‑基因表达谱矩阵A中各基因的表达量。

[0134] 进行第 次插补后，细胞‑基因表达谱矩阵 中第 个细胞第 个基因的表达量 会变为细胞‑基因表达谱矩阵 中第 个细胞的各个相邻细胞的该基因的表达量乘第 个细胞转换为该各个相邻细胞的概率（幂次为）之和：

[0135] ；

[0136] 由于马尔可夫转移概率矩阵 的特征值均为[0,1]之间，特征值通过指数运算会逐渐减小，其范围也在[0,1]之间。随着马尔可夫转移概率矩阵 多次求幂，除1以外的所有特征值的大小不断减小，如此能够降低噪声的重要性，其解释能力接近于零。随着 的增加，细胞从它们的邻居细胞那里学习缺失的值，迅速获得生物学上非常相似的细胞之间的关系。

[0137] 当噪声去除转变为真实生物信息信号去除，则 达到最佳迭代次数。由于噪声通常与真实信号本身具有不同的频率（即分别为高频和低频），因此随着高频信息的去除，数据将发生快速变化，然后变为较慢地变化，达到收敛。当数据变化率低于首次低于阈值（本实施例设置为5%）时，可以认为达到收敛，达到最佳迭代次数。最终获得细胞‑基因加权表达插补矩阵 ，以降低数据噪声并有效计算缺失的转录本表达量，而不会过度拟合数据。

[0138] 细胞 的数据变化率通过决定系数 进行量化：

[0139] 对于迭代次数 ，先计算插补后基因 的相对表达量（基因表达量除以全部基因表达量的总和），

[0140] ；

[0141] 再计算细胞 插补后的平均基因表达量：

[0142] ；

[0143] 再基于第 次插补后和插补前的全部基因相对表达量的残差平方和（ ）与插补后的全部基因相对表达量的离差平方和（ ）计算细胞 的决定系数 ，公式如下：

[0144] ；

[0145] ；

[0146] ；

[0147] 在本实施中，当 时，所有细胞的决定系数均 。由此获得了加权插补后的细胞‑基因表达谱矩阵 ，有效降低了数据噪声并有效地填补了缺失的转录本表达量，并且不会过度拟合数据。

[0148] 利用多重加权插补前后的细胞‑基因表达谱矩阵 分别进行分析，发现加权插补后相同基因的表达量在相邻细胞中更加一致，更加真实。例如对于T细胞的marker基因CD3D与B细胞的marker基因MS4A1，基因表达情况对比分别如图3和图4所示。可以看出，进行多重加权插补后的细胞‑基因表达谱矩阵 中，基因表达量在相邻的细胞间表达更加清晰，可以直观判断细胞类型的特征并确定具体的细胞类型。相反，在多重加权插补前，基因表达量即使是同一cluster中细胞中一致性也较差，难以基于细胞类型特异marker基因的表达量区分不同的细胞群。

[0149] 另外，由于细胞中的各基因表达量随着真实信号的增强，各个细胞中基因的稀疏性降低，表达量为0的细胞数减少，基因之间的相关作用随之增强。如CD14基因和ITGAM基因均为髓系细胞的marker基因，在正常的髓系细胞中该两种基因表达相关性较强。但是利用多重加权插补前的细胞‑基因表达谱矩阵 难以发现这种相关性，利用多重加权插补后的细胞‑基因表达谱矩阵 进行分析，可以看出两者的相关性明显地展现出来，如图5所示。因此，利用本发明的方法，对于细胞类型的鉴定也起到非常积极的作用。

[0150] 实施例 2 单细胞测序基因矩阵稀疏性处理的优化

[0151] 本实施例对实施例1的方案进行进一步优化，具体地：

[0152] 在计算马尔可夫转移概率矩阵时，进行如下处理：

[0153] 对于同一区域（簇、类型）中任意两个细胞，确定与该两个细胞距离之和最小的15个细胞作为该两个细胞的邻居细胞，若不足15个，则全部其余细胞为邻居细胞。

[0154] 对于每个细胞，将其与非邻居细胞的马尔可夫转移概率设置为0。

[0155] 利用本实施例优化后的马尔可夫转移概率矩阵对基因‑细胞表达谱矩阵进行加权多重插补，进一步进行分析。

[0156] 同样，对于T细胞的marker基因CD3D与B细胞的marker基因MS4A1，基因表达情况对比分别如图6和图7所示。可见优化后的方法，基因表达量在相邻的细胞间表达进一步清晰，可以直观判断细胞类型的特征并确定具体的细胞类型。对于细胞之间的相互作用，同样如CD14基因和ITGAM基因，两者的相关性进一步得到加强，如图8所示。

[0157] 实施例 3 单细胞测序基因矩阵稀疏性处理的进一步优化

[0158] 本实施例对实施例2的方案进行进一步优化，具体地：

[0159] 在利用马尔可夫转移概率矩阵对基因‑细胞表达谱矩阵进行加权多重插补步骤中，对于马尔可夫转移概率矩阵 进行求幂时，当幂次为时，若某个特征值低于0.001，则将该特征值赋值0，再进行插补。

[0160] 利用本实施例优化后的方法对基因‑细胞表达谱矩阵进行加权多重插补后进行分析。结果发现，优化后的方法能够去除数据中的大部分噪声。

[0161] 同样，对于T细胞的marker基因CD3D与B细胞的marker基因MS4A1，基因表达情况对比分别如图9和图10所示。可见优化后的方法，基因表达量在相邻的细胞间表达进一步清晰，可以直观判断细胞类型的特征并确定具体的细胞类型。对于细胞之间的相互作用，同样如CD14基因和ITGAM基因，两者的相关性进一步得到加强，如图11所示。

[0162] 在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。