一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法
阅读说明：本技术 一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法 (Single cell transcriptome sequencing data clustering recommendation method based on two-dimensional distribution structure judgment ) 是由 李敏 田宇 郑瑞清 于 2021-01-18 设计创作，主要内容包括：本发明公开了一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法,包括：获取多个细胞的单细胞转录组测序数据得到的基因表达矩阵,在过滤和标准化处理后,构建二维特征矩阵并进行线性归一化；根据归一化后的二维特征矩阵计算细胞间的欧式距离,从而建立细胞最小生成树；通过自适应阈值对细胞最小生成树进行切割,并以切割后构成的簇的平衡性来确定数据的二维分布结构；对具有模糊簇间边界和连续二维分布结构的数据,推荐并应用层次聚类算法,而对具有明显簇间边界和分块二维分布结构的数据,推荐并应用谱聚类算法。本发明能为单细胞转录组测序数据在层次聚类和谱聚类中推荐更适合其二维分布结构的方法作为下游聚类分析方法,提高聚类准确性。(The invention discloses a single cell transcriptome sequencing data clustering recommendation method based on two-dimensional distribution structure judgment, which comprises the following steps: acquiring gene expression matrixes obtained from single-cell transcriptome sequencing data of a plurality of cells, and after filtering and standardization processing, constructing a two-dimensional characteristic matrix and performing linear normalization; calculating the Euclidean distance between cells according to the normalized two-dimensional feature matrix, thereby establishing a minimum spanning tree of the cells; cutting the minimum spanning tree of the cells through an adaptive threshold, and determining a two-dimensional distribution structure of data according to the balance of clusters formed after cutting; a hierarchical clustering algorithm is recommended and applied to data with fuzzy inter-cluster boundaries and a continuous two-dimensional distribution structure, and a spectral clustering algorithm is recommended and applied to data with obvious inter-cluster boundaries and a blocking two-dimensional distribution structure. The method can recommend a method more suitable for the two-dimensional distribution structure of the sequencing data of the single-cell transcriptome in hierarchical clustering and spectral clustering as a downstream clustering analysis method, and improves the clustering accuracy.)
技术领域
本发明涉及生物信息学领域，涉及一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法。
背景技术
在细胞生物学领域，单细胞分析是在单细胞水平上对基因组学、转录组学、蛋白质组学和代谢组学的研究。它提供了一种超灵敏的工具来阐明特定的分子机制和途径，并揭示了细胞异质性的本质。随着技术的发展和成本的下降，应用于单细胞全基因组的转录组测序(scRNA-seq)技术正在迅速成为生物学和生物医学研究等许多领域的选择。在单细胞分辨率下研究全基因组的基因表达克服了传统RNA测序的内在局限性，单细胞转录组测序使研究人员能够更严格地解决有关组织的细胞组成，转录物异质性和细胞类型的问题，使单细胞转录组测序已用于研究癌症，宏基因组学以及调控和进化网络。近些年来，单细胞异质性分析研究成为了生物信息领域的一大研究热点。为了剖析细胞异质性，就有必要在单细胞水平进行基因表达分析。单细胞转录组测序技术能获得单个细胞内近万个基因的表达信息，为辨别生物组织中各种细胞类型的转录特征和全面揭示细胞之间的基因表达差异提供了有力的工具。在单细胞转录组测序数据分析中，鉴定种群结构是至关重要的一步。对单细胞转录组测序数据进行的下游分析涉及多个步骤，包括质量控制，定量，归一化，聚类，寻找轨迹和鉴定差异表达的基因等。目前已经提出了很多单细胞聚类方法来对细胞进行划分，其中层次聚类和谱聚类是下游聚类分析中应用最多的两个聚类方法：
层次聚类：将单个数据迭代的合并为较大的簇(凝聚聚类)或将一个大类迭代的划分为较小的簇(分裂聚类)，并通过得到的树状图中层次树的分支结构来表示最终的聚类结果。层次聚类可以通过切割树状图来获得不同的聚类结果，但方法缺少对簇间边界的精准定义，这将导致其在具有模糊边界结构的数据上的聚类效果不够理想。
谱聚类：将所有数据看作空间中的点，点之间用边连接起来。距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，通过对所有数据点组成的图进行切图，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的。因此，谱聚类可以处理更为复杂的数据分布结构，例如模糊边界问题，但方法缺点在于十分依赖相似性矩阵的准确性。
两种聚类方法由于基于不同的理论和策略，因而可能在具有不同分布结构的数据上的聚类结果存在差异。在通过分析每组数据的二维分布结构以选择更合适的聚类方法方面，还有进一步提升的空间。
因此，有必要提供一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法。
发明内容
本发明所要解决的技术问题是，针对现有技术不足，提供一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法，能够为具有不同二维分布结构的数据在层次聚类和谱聚类中推荐更合适的聚类方法，提高聚类准确性。
为解决上述技术问题，本发明所采用的技术方案是：
一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法，包括以下步骤：
步骤1，获取N个细胞的单细胞转录组测序数据，通过统计基因在每个细胞中的表达量得到基因表达矩阵，并进行过滤和标准化处理；然后通过均匀流形近似和投影构建二维特征矩阵，并将二维特征矩阵进行线性归一化；
步骤2，根据归一化后的二维特征矩阵，计算每两个细胞间的欧式距离，并建立基于距离矩阵的细胞最小生成树；
步骤3，通过自适应阈值对构建的细胞最小生成树进行切割，并以切割后构成的簇的平衡性来确定单细胞转录组测序数据的二维分布结构；所述自适应阈值根据细胞数N定义得到；
步骤4，根据数据二维分布结构的差异并结合层次聚类和谱聚类的方法特性，对计算确定具有模糊簇间边界和连续二维分布结构的数据，推荐并使用层次聚类算法进行下游聚类分析，而对计算确定具有明显簇间边界和分块二维分布结构的数据，推荐并使用谱聚类算法进行下游聚类分析。
在更优的技术方案，所述步骤1对二维特征矩阵进行线性归一化的公式为：
其中，Y′表示通过均匀流形近似和投影构建的二维特征矩阵，Ymax和Ymin分别表示初始基因表达矩阵中表达量的最大值和最小值。
在更优的技术方案，所述步骤2对归一化后的二维特征矩阵计算细胞间距离的公式为：
其中，dij表示细胞i到细胞j之间的欧式距离，yim和yjm分别为细胞i和细胞j的第m个基因的表达量。
在更优的技术方案，所述步骤2中，基于距离矩阵，通过快速EMST双树Boruvka算法，建立基于距离矩阵的细胞最小生成树。
在更优的技术方案，所述步骤3中根据细胞数N定义自适应阈值的公式为：
Threshold＝λ×pCR
其中，Threshold为定义的自适应阈值，λ＝5为默认常数值，pCR为单个细胞分辨率，且
在更优的技术方案，所述步骤3中，若细胞最小生成树中所有边长均小于阈值，则确定数据为具有模糊簇间边界和连续的二维分布结构；反之若存在大于阈值的树边，则将满足条件的边进行切割并计算重新构成的簇间的平衡性，当最大簇与次大簇中细胞的数量比超过预设比例时，确定数据为具有模糊簇间边界和连续的二维分布结构，否则将数据确定为具有明显簇间边界和分块的二维分布结构。
在更优的技术方案，所述预设比例为8：2。
有益效果
本发明能为单细胞转录组测序数据集在层次聚类和谱聚类中推荐更适合其二维分布结构的方法作为下游聚类分析方法，并提高聚类准确性。
附图说明
图1为本发明的流程图；
图2、图3为本发明针对实施例中6套数据集确定的数据二维分布结构与使用tSNE和UMAP两种可视化方法在具有真实标签的分化数据集和亚型数据集上进行二维可视化的结果比较图；
图4、图5分别为本发明针对实施例中另外6套数据集确定的数据二维分布结构与使用tSNE和UMAP两种可视化方法在具有真实标签的分化数据集和亚型数据集上进行二维可视化的结果比较图。
图6和图7为用选择的10个经典相关方法的预处理策略，分别结合层次聚类和谱聚类作为下游聚类方法所得到的聚类结果的NMI比较图。
具体实施方式
下面对本发明的实施例作详细说明，本实施例以本发明的技术方案为依据开展，给出了详细的实施方式和具体的操作过程，对本发明的技术方案作进一步解释说明。
本实施例提供一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法，包括以下步骤：
步骤1，获取N个细胞的单细胞转录组测序数据，通过统计测序数据中基因在每个细胞中的表达量得到基因表达矩阵X＝[x1,x2,…,xN]，xi＝[xi1,xi2,…,xim],i＝1,2,…,N，m表示细胞中的基因数，xi1,xi2,…,xim表示细胞i分别在m个基因中的表达量；将基因表达矩阵X中表达量为0的基因删除以完成过滤后，再对过滤后的基因表达矩阵进行标准化处理；再通过对标准化处理后的基因表达矩阵进行均匀流形近似和投影[1]，以构建二维特征矩阵Y，并将二维特征矩阵Y按以下公式进行线性归一化处理：
其中，Y′为归一化后的二维特征矩阵，Ymax和Ymin分别表示初始基因表达矩阵中表达量矩阵表达值的最大值和最小值。
步骤2，根据归一化后的二维特征矩阵Y′，计算每两个细胞间的欧式距离；然后基于欧式距离构成的距离矩阵，通过快速EMST双树Boruvka(fast EMST Dual-Tree Boruvka)算法[2]，建立基于距离矩阵的细胞最小生成树；
每两个细胞间的欧式距离计算公式为：
其中，dij表示细胞i到细胞j之间的欧式距离，yim和yjm分别为细胞i和细胞j的第m个基因的表达量。
基于欧式距离矩阵的快速EMST双树Boruvka算法的步骤为：
其中，Cq表示包含q的分支，d(Cq)表示Cq中当前最近邻的距离(初始值为无穷大)，e(Cq)表示Cq与其最近邻构成的边集合，d(Q,R)表示Q与R中结点的最小距离。
步骤3，通过自适应阈值对构建的细胞最小生成树进行切割，并以切割后构成的簇的平衡性来确定单细胞转录组测序数据的二维分布结构；
其中的自适应阈值，根据细胞数N按以下公式定义得到：
Threshold＝λ×pCR
其中，Threshold为定义的自适应阈值，λ＝5为默认常数值，pCR为单个细胞分辨率(定义为随机平铺在1x1平面的每个细胞间距离)，本实施例中取
另外，本实施例中单细胞转录组测序数据的二维分布结构的确定方法具体为：若细胞最小生成树中所有欧式距离计算所得的边长值均小于阈值，则确定数据为具有模糊簇间边界和连续的二维分布结构；反之若存在大于阈值的树边，则将满足条件的边进行切割并计算重新构成的簇间的平衡性，当最大簇与次大簇中细胞的数量比超过预设比例8：2时，确定数据为具有模糊簇间边界和连续的二维分布结构，否则将数据确定为具有明显簇间边界和分块的二维分布结构。
步骤4，根据数据二维分布结构的差异并结合层次聚类和谱聚类的方法特性，对计算确定具有模糊簇间边界和连续二维分布结构的数据，推荐并使用层次聚类算法进行下游聚类分析，而对计算确定具有明显簇间边界和分块二维分布结构的数据，推荐并使用谱聚类算法进行下游聚类分析。
[1]McInnes L,Healy J and Melville J.UMAP:uniform manifoldapproximation and projection.arXiv 2018；1802.03426.
[2]March WB,Ram P and Gray AG.Fast euclidian minimum spanning tree:algorithm analysis,and applications.16th ACM SIGKDD International Conferenceon Knowledge Discovery and Data mining 2010.
实验验证：
为了验证本发明方法的有效性，从已发表文献和公共数据库(GEO、ArrayExpress)收集了12套来自不同物种或不同组织或不同规模的单细胞数据集，包括细胞的基因表达矩阵和真实类别标签，并在这12套真实数据集上测试本发明提出的方法。12套数据集的详细信息如表1所示。为客观评价本发明方法的有效性，选用常用指标NMI作为聚类评价标准，同时选择10个经典的分别以相似性计算、填充和降维方法作为预处理策略的单细胞相关方法来验证推荐的准确性，其中包括SPEARMAN、UMAP、SC3、RAFSIL、SIMLR、MPSSC、CIDR、SINCERA、SEURAT和ZIFA。对两组不同的标签F1和F2，NMI的计算公式为：
其中，I(F1,F2)表示F1和F2两者的互信息，H(F1)、H(F2)分别表示两者的信息熵。
为了评价本发明方法分类的准确性及功能可解释性，进行以下三项分析：
1、分析确定的数据二维分布结构及推荐的聚类方法
该部分针对12套数据集中的样本分别进行分析。使用本发明确定每套数据集的二维分布结构和推荐的聚类方法。实验结果如表2所示。
2、分析确定的结构与可视化的一致性
该部分针对12套数据集中的样本分别进行分析。对每套数据集，分别通过tSNE和UMAP结合数据真实标签进行二维可视化，结果如图2至图5所示。从图2至图5可以看出，两种可视化的结果均与本发明确定的数据二维分布结构显示出一致性。
3、分析推荐的聚类方法的准确性
该部分针对12套数据集中的样本分别进行分析。使用选择的10个经典相关方法的预处理策略，分别结合层次聚类和谱聚类作为下游聚类方法。采用NMI指标评价聚类结果。NMI的实验结果如图6和图7所示，从图6和图7可以看出，本发明能够对12套数据集分别从层次聚类和谱聚类中推荐更合适的聚类方法，并以此进一步提高聚类的准确性。
实验结果表明，本发明方法具有很好的鲁棒性和准确性。
表1数据集信息
数据集
细胞数
基因数
类别数
计量单位
Ting
114
14405
5
RPM
Buettner
182
8989
3
FPKM
Pollen
249
14805
11
TPM
Ginhoux
251
11834
3
RPKM
LaManno
337
14703
13
UMI
Darmanis
420
22085
8
CPM
Leng
460
19084
4
TPM
Camp
465
18999
6
FPKM
Gokce
1208
16379
10
TPM
Nestorowa
1645
3991
3
UMI
Close
1733
23045
4
TPM
Zeisel
3005
4412
9
UMI
表2本方法确定的数据二维分布结构和推荐方法
数据集
细胞数
阈值
切割边数
簇平衡性
数据分布结构
推荐方法
Ting
114
0.4683
1
平衡
分块结构
层次聚类
Buettner
182
0.3706
0
-
连续结构
谱聚类
Pollen
249
0.3169
2
平衡
分块结构
层次聚类
Ginhoux
251
0.3156
0
-
连续结构
谱聚类
LaManno
337
0.2724
0
-
连续结构
谱聚类
Darmanis
420
0.2440
2
平衡
分块结构
层次聚类
Leng
460
0.2331
0
-
连续结构
谱聚类
Camp
465
0.2319
2
平衡
分块结构
层次聚类
Gokce
1208
0.1439
4
平衡
分块结构
层次聚类
Nestorowa
1645
0.1233
1
不平衡
连续结构
谱聚类
Close
1733
0.1201
1
不平衡
连续结构
谱聚类
Zeisel
3005
0.0912
3
平衡
分块结构
层次聚类
以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保护的范围之内。