基于拓扑特性的稀有类型细胞检测方法、系统及设备

1.本公开实施例涉及数据处理技术领域，尤其涉及一种基于拓扑特性的稀有类型细胞检测方法、系统及设备。


背景技术：

2.目前，与传统测序技术相比，单细胞测序技术通过提供高分辨率测序数据极大程度的提高生命医学领域的研究精度，为一些重要的医学问题提供新的解决方案，帮助研究人员实现从细胞水平对生物体的发育及疾病进展深入探索。发现稀有细胞类型，例如短暂存在的前体细胞、内皮祖细胞、抗原特异性t细胞、干细胞、循环肿瘤细胞或是癌细胞、不变型自然杀伤性t细胞等对于患者的临床诊断治疗有很大的帮助，对于深入理解正常和疾病状态下的组织生物学具有非常重要的意义。目前单细胞转录组数据的分析流程已经比较成熟，近年来已经涌现出大量针对单细胞转录组数据的聚类算法，他们能够较好的发现其中常见的细胞类型以进行后续标志物的筛选。然而研究一个组织或器官的发育变化及功能就需要清晰的知道其中所有的细胞类型，仅仅只了解其中的常见细胞类型是远远不够的。目前绝大多数异质性分析算法都是针对常见细胞类型进行分析的，因为常见细胞类型的细胞数量较多且表达量高，往往可以通过对其高可变基因，即细胞与细胞之间表达量差别最大的基因计算相似性进行聚类便可以达到较为理想的效果。对于一些稀有细胞类型来说，常见细胞类型的高可变基因特征选择并不适用，但是鉴定出在稀有类型细胞中特异表达并且可以进行类型划分的标记基因目前仍存在较大的挑战。而且稀有类型细胞相比于常见类型细胞会更容易受到样本质量或者批次效应等多种因素的干扰，不仅如此，大多数稀有细胞类型混合在常见细胞类型中导致稀有细胞类型和常见细胞类型之间很难区分。
3.目前绝大部分鉴定稀有类型细胞的算法都依赖于单细胞聚类算法，这样往往会暴露一些耗时久，假阳率高且准确率低等问题。目前主流的稀有类型细胞检测算法giniclust3、raceid3、cellsius、fire和gapclust及其优缺点。其中giniclust3、raceid3和cellsius算法核心都是基于聚类，所以只能输出单细胞表达谱的二元预测结果。giniclust3的设计思路是基于基尼指数选择基尼基因进行聚类。当单细胞表达谱中生物复杂度高时，giniclust3无法考虑稀有细胞类型簇内的关系，所以并不适用于稀有细胞亚型划分。cellsius是基于k-means先进行粗聚类再对类内进行细分，其预测性能较差且假阳率高。raceid3也是基于k-means算法设计，但是其耗时极高，当数据集细胞数量达到5k时，其计算消耗时间超过78h。fire和gapclust可以输出单细胞表达谱的连续稀有性分数，其中fire是基于随机投影得到的哈希码为每个细胞计算稀有性分数，但是其特征选择过程需要消耗大量的时间和内存。gapclust是基于近邻距离变化的二阶偏导为每个细胞计算稀有性分数，其参数近邻位置的选择可以直接大幅度影响预测结果精度。由此可见gapclust方法的参数敏感性较强且鲁棒性差，实验结果会较大程度受到数据集中噪声点的影响。
4.目前主流的稀有细胞类型检测方法存在以下问题：(1)如何设计快速高效的稀有细胞检测方法；(2)如何结合其他生物信息，提高稀有细胞类型划分的精度。(3)如何降低稀
有类型细胞检测的假阳性；(4)如何提高稀有细胞检测的召回率和准确率。
5.可见，亟需一种高效精准、适应性强的基于拓扑特性的稀有类型细胞检测方法。


技术实现要素：

6.有鉴于此，本公开实施例提供一种基于拓扑特性的稀有类型细胞检测方法、系统及设备，至少部分解决现有技术中存在检测效率、精准度和适应性较差的问题。
7.第一方面，本公开实施例提供了一种基于拓扑特性的稀有类型细胞检测方法，包括：
8.步骤1，输入一个单细胞表达谱、n个细胞和m个基因，通过数据预处理保留表达散度排名靠前的m1个高可变基因构成特征矩阵，并对得到的特征矩阵计算knn最近邻距离矩阵x；
9.步骤2，对得到n个细胞与k个近邻细胞距离构成的近邻距离矩阵x，根据每个细胞与其他k个近邻细胞的距离向量，计算其熵值h(x)、平均值中位数m
0.5
和偏度g作为该细胞的四个领域拓扑特性，并以此构建细胞领域拓扑特征矩阵r；
10.步骤3，基于拓扑特征矩阵r，结合局部异常因子检测算法计算每个细胞的稀有性分数，并根据自适应参数优化方法设定阈值，输出稀有类型细胞的二元预测结果；
11.步骤4，将稀有类型细胞的二元预测结果结合社区发现算法进一步对稀有细胞进行亚群划分。
12.根据本公开实施例的一种具体实现方式，所述步骤2具体包括：
13.计算熵值拓扑特征h(x)的公式表示为：
[0014][0015]
其中p(x
ij
)表示第i个细胞的第j位近邻距离的概率，其公式表示为：
[0016][0017]
领域拓扑特征代表平均值拓扑特征，其公式表示为：
[0018][0019]
采用x
0.5
代表中位数拓扑特征，其计算方式表示为:
[0020]
当近邻数k为奇数时，
[0021][0022]
当近邻数k为偶数时，
[0023][0024]
采用skew(x)代表偏度拓扑特征，其计算方式表示为：
[0025]
skew(xi)＝g([x
i1
,x
i2
,x
i3
,
…
,x
ik
]),i∈(1,2,
…
,n)
[0026]
其中k表示近邻数，n表示细胞数量，x
ij
表示第i个细胞的第j个近邻距离，h(x)表示
熵值，表示平均值，m
0.5
表示中位数和g表示偏度；
[0027]
根据以上的四个拓扑特征构建拓扑特征矩阵r：
[0028][0029]
根据本公开实施例的一种具体实现方式，所述局部异常因子检测算法的表达式为表示细胞i领域点nk(i)的局部可达密度与细胞i的局部可达密度之比的平均数；
[0030]
其中，lrd(i)表示细胞i的局部可达密度，公式表示为：
[0031][0032]
可达距离distancek(i,p)表示为：
[0033]
distancek(i,p)＝max{distancek(p),d(i,p)}
[0034]
其中，distancek(i,p)表示细胞p到细胞i的第k可达距离，至少为p的k近邻距离或者是i和p的距离，d(i,p)表示为i和p之间的距离。
[0035]
根据本公开实施例的一种具体实现方式，所述步骤3具体包括：
[0036]
利用数据驱动的方式输入预设的稀有细胞比例；
[0037]
将输出所有细胞的连续稀有性分数并将其从大到小进行排序，绘制散点图；
[0038]
根据散点图计算细胞i和其相邻细胞之间的稀有性分数差值δwi；
[0039]
根据稀有性分数差值δwi计算稀有性分数变化率即二阶偏导数δδwi，并将δδwi的极值点位置的横坐标即稀有度分数排名靠前的细胞个数。
[0040]
第二方面，本公开实施例提供了一种基于拓扑特性的稀有类型细胞检测系统，包括：
[0041]
预处理模块，用于输入一个单细胞表达谱、n个细胞和m个基因，通过数据预处理保留表达散度排名靠前的m1个高可变基因构成特征矩阵，并对得到的特征矩阵计算knn最近邻距离矩阵x；
[0042]
计算模块，用于对得到n个细胞与k个近邻细胞距离构成的近邻距离矩阵x，根据每个细胞与其他k个近邻细胞的距离向量，计算其熵值h(x)、平均值中位数m
0.5
和偏度g作为该细胞的四个领域拓扑特性，并以此构建细胞领域拓扑特征矩阵r；
[0043]
预测模块，用于基于拓扑特征矩阵r，结合局部异常因子检测算法计算每个细胞的稀有性分数，并根据自适应参数优化方法设定阈值，输出稀有类型细胞的二元预测结果；
[0044]
划分模块，用于将稀有类型细胞的二元预测结果结合社区发现算法进一步对稀有细胞进行亚群划分。
[0045]
第三方面，本公开实施例还提供了一种电子设备，该电子设备包括：
[0046]
至少一个处理器；以及，
[0047]
与该至少一个处理器通信连接的存储器；其中，
[0048]
该存储器存储有可被该至少一个处理器执行的指令，该指令被该至少一个处理器执行，以使该至少一个处理器能够执行前述第一方面或第一方面的任一实现方式中的基于拓扑特性的稀有类型细胞检测方法。
[0049]
本公开实施例中的基于拓扑特性的稀有类型细胞检测方案，包括：步骤1，输入一个单细胞表达谱、n个细胞和m个基因，通过数据预处理保留表达散度排名靠前的m1个高可变基因构成特征矩阵，并对得到的特征矩阵计算knn最近邻距离矩阵x；步骤2，对得到n个细胞与k个近邻细胞距离构成的近邻距离矩阵x，根据每个细胞与其他k个近邻细胞的距离向量，计算其熵值h(x)、平均值中位数m
0.5
和偏度g作为该细胞的四个领域拓扑特性，并以此构建细胞领域拓扑特征矩阵r；步骤3，基于拓扑特征矩阵r，结合局部异常因子检测算法计算每个细胞的稀有性分数，并根据自适应参数优化方法设定阈值，输出稀有类型细胞的二元预测结果；步骤4，将稀有类型细胞的二元预测结果结合社区发现算法进一步对稀有细胞进行亚群划分。
[0050]
本公开实施例的有益效果为：通过本公开的方案，从稀有类型细胞的近邻差异出发构建拓扑特征矩阵，其优点是有效降低了数据维度，提升了预测效率、精准度和适应性。
附图说明
[0051]
为了更清楚地说明本公开实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。
[0052]
图1为本公开实施例提供的一种基于拓扑特性的稀有类型细胞检测方法的流程示意图；
[0053]
图2为本公开实施例提供的另一种基于拓扑特性的稀有类型细胞检测方法的流程示意图；
[0054]
图3为本公开实施例提供的一种利用本方法在真实数据集上的稀有细胞类型的异质性分析性能示意图；
[0055]
图4为本公开实施例提供的一种基于拓扑特性的稀有类型细胞检测系统的结构示意图；
[0056]
图5为本公开实施例提供的电子设备示意图。
具体实施方式
[0057]
下面结合附图对本公开实施例进行详细描述。
[0058]
以下通过特定的具体实例说明本公开的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本公开的其他优点与功效。显然，所描述的实施例仅仅是本公开一部分实施例，而不是全部的实施例。本公开还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本公开的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。基于本公开中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本公开保护的范围。
[0059]
需要说明的是，下文描述在所附权利要求书的范围内的实施例的各种方面。应显
而易见，本文中所描述的方面可体现于广泛多种形式中，且本文中所描述的任何特定结构及/或功能仅为说明性的。基于本公开，所属领域的技术人员应了解，本文中所描述的一个方面可与任何其它方面独立地实施，且可以各种方式组合这些方面中的两者或两者以上。举例来说，可使用本文中所阐述的任何数目个方面来实施设备及/或实践方法。另外，可使用除了本文中所阐述的方面中的一或多者之外的其它结构及/或功能性实施此设备及/或实践此方法。
[0060]
还需要说明的是，以下实施例中所提供的图示仅以示意方式说明本公开的基本构想，图式中仅显示与本公开中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。
[0061]
另外，在以下描述中，提供具体细节是为了便于透彻理解实例。然而，所属领域的技术人员将理解，可在没有这些特定细节的情况下实践所述方面。
[0062]
本公开实施例提供一种基于拓扑特性的稀有类型细胞检测方法，所述方法可以应用于生命医学领域的单细胞检测过程中。
[0063]
参见图1，为本公开实施例提供的一种基于拓扑特性的稀有类型细胞检测方法的流程示意图。如图1和图2所示，所述方法主要包括以下步骤：
[0064]
步骤1，输入一个单细胞表达谱、n个细胞和m个基因，通过数据预处理保留表达散度排名靠前的m1个高可变基因构成特征矩阵，并对得到的特征矩阵计算knn最近邻距离矩阵x；
[0065]
具体实施时，可以输入一个单细胞表达谱，n个细胞，m个基因，通过数据预处理保留表达散度排名靠前的m1个高可变基因构成特征矩阵，基于上一步得到的特征矩阵计算knn最近邻(k-nearest neighbor，简称knn)。朴素的knn计算方法是计算样本点和数据集中其他所有样本点的距离，然后排序选择距离最小的样本。假设本发明计算每个细胞的k个近邻细胞(不包含自己)，当细胞数目较多时朴素的计算方法消耗的时间将呈线性增长。因此topo采用轻量快速的hnsw(hierarchical navigable small world)近邻搜索算法构建knn图。hnsw近邻搜索算法是传统的构建knn图的基础上，通过在构建过程中保留一些非近邻的边，从而实现相距较远点之间的跳跃，并且通过“跳表结构”实现多层次的近邻搜索。经过这一步骤处理可以得到近邻距离矩阵x：
[0066][0067]
矩阵x表示单细胞表达谱中n个细胞与k个近邻细胞距离构成的近邻距离矩阵，其中行向量[x
i1
,x
i2
,x
i3
,
…
,x
ik
],i∈(1,2,
…
,n)表示第i个细胞与其他k个近邻细胞的距离向量其knn近邻计算方法采用轻量快速的hnsw方法。
[0068]
步骤2，对得到n个细胞与k个近邻细胞距离构成的近邻距离矩阵x，根据每个细胞与其他k个近邻细胞的距离向量，计算其熵值h(x)、平均值中位数m
0.5
和偏度g作为该细胞的四个领域拓扑特性，并以此构建细胞领域拓扑特征矩阵r；
[0069]
进一步的，所述步骤2具体包括：
[0070]
计算熵值拓扑特征h(x)的公式表示为：
[0071][0072]
其中p(x
ij
)表示第i个细胞的第j位近邻距离的概率，其公式表示为：
[0073][0074]
领域拓扑特征代表平均值拓扑特征，其公式表示为：
[0075][0076]
采用x
0.5
代表中位数拓扑特征，其计算方式表示为:
[0077]
当近邻数k为奇数时，
[0078][0079]
当近邻数k为偶数时，
[0080][0081]
采用skew(x)代表偏度拓扑特征，其计算方式表示为：
[0082]
skew(xi)＝g([x
i1
,x
i2
,x
i3
,
…
,x
ik
]),i∈(1,2,
…
,n)
[0083]
其中k表示近邻数，n表示细胞数量，x
ij
表示第i个细胞的第j个近邻距离，h(x)表示熵值，表示平均值，m
0.5
表示中位数和g表示偏度；
[0084]
根据以上的四个拓扑特征构建拓扑特征矩阵r：
[0085][0086]
具体实施时，从上一步可以得到n个细胞与k个近邻细胞距离构成的近邻距离矩阵x，取矩阵每一行[x
i1
,x
i2
,x
i3
,
…
,x
ik
],i∈(1,2,
…
,n)，即每个细胞与其他k个近邻细胞的距离向量，对这个向量计算熵值h(x)，平均值中位数m
0.5
和偏度g构成四列拓扑矩阵特征，其中熵值拓扑特征h(x)的公式表示为：
[0087][0088]
其中p(x
ij
)表示第i个细胞的第j位近邻距离的概率，其公式表示为：
[0089][0090]
领域拓扑特征代表平均值拓扑特征，其公式表示为：
[0091][0092]
采用x
0.5
代表中位数拓扑特征，其计算方式表示为:
[0093]
a.如果近邻数k为奇数，
[0094][0095]
b.当近邻数k为偶数时，
[0096][0097]
采用skew(x)代表偏度拓扑特征，其计算方式表示为：
[0098]
skew(xi)＝g([x
i1
,x
i2
,x
i3
,
…
,x
ik
]),i∈(1,2,
…
,n)
[0099]
其中k表示近邻数，n表示细胞数量，x
ij
表示第i个细胞的第j个近邻距离。h(x)表示熵值，表示平均值，m
0.5
表示中位数和g表示偏度。
[0100]
根据以上基于熵值，平均数，中位数和偏度的四种计算公式，最终得到了行为n,列为4的n
×
4的拓扑特征矩阵r:
[0101][0102]
步骤3，基于拓扑特征矩阵r，结合局部异常因子检测算法计算每个细胞的稀有性分数，并根据自适应参数优化方法设定阈值，输出稀有类型细胞的二元预测结果；
[0103]
进一步的，所述局部异常因子检测算法的表达式为进一步的，所述局部异常因子检测算法的表达式为表示细胞i领域点nk(i)的局部可达密度与细胞i的局部可达密度之比的平均数；
[0104]
其中，lrd(i)表示细胞i的局部可达密度，公式表示为：
[0105][0106]
可达距离distancek(i,p)表示为：
[0107]
distancek(i,p)＝max{distancek(p),d(i,p)}
[0108]
其中，distancek(i,p)表示细胞p到细胞i的第k可达距离，至少为p的k近邻距离或者是i和p的距离，d(i,p)表示为i和p之间的距离。
[0109]
在上述实施例的基础上，所述步骤3具体包括：
[0110]
利用数据驱动的方式输入预设的稀有细胞比例；
[0111]
将输出所有细胞的连续稀有性分数并将其从大到小进行排序，绘制散点图；
[0112]
根据散点图计算细胞i和其相邻细胞之间的稀有性分数差值δwi；
[0113]
根据稀有性分数差值δwi计算稀有性分数变化率即二阶偏导数δδwi，并将δδwi的极值点位置的横坐标即稀有度分数排名靠前的细胞个数。
[0114]
具体实施时，基于拓扑特征矩阵r，已经可以从特征上可以区分出稀有细胞和常见细胞，所以这一步仅需要采用一种简单的异常值检测算法，将细胞进行二分类。本发明采用
的是异常值检测算法lof(local outlier factor)，这是一种基于密度的无监督异常点检测算法，相比于其他异常值检测算法，lof并不需要要求数据满足任何分布，并且可以量化每个数据点的异常程度(outlierness)。lof算法的主要思想是给定数据点局部邻域数据点比较密集，那么它为正常点的可能性较高，即常见类型细胞，若该点邻域中数据点较为稀疏，密度较小，那么它是异常点的概率较高，即稀有细胞。其中lof公式表示为：
[0115][0116]
上式表示细胞i领域点nk(i)的局部可达密度(local reachability density)与细胞i的局部可达密度之比的平均数。如果这个比值越接近1，说明i与领域点密度相近，i可能和领域同属一簇；如果这个比值越小于1，说明i的密度高于其领域点密度，i为密集点；如果这个比值越大于1，说i的密度小于其领域点密度，i越可能是异常点。其中，lrd(i)表示细胞i的局部可达密度，公式表示为：
[0117][0118]
lrd(i)表示密度值，如果细胞i和领域细胞属于同一簇，那么可达距离较小，密度值较高；如果可达距离较大，那么密度值较小，细胞i是离群点的可能性越大。其中可达距离distancek(i,p)表示为：
[0119]
distancek(i,p)＝max{distancek(p),d(i,p)}
[0120]
其中distancek(i,p)表示细胞p到细胞i的第k可达距离，至少为p的k近邻距离，或者是i和p的距离。d(i,p)表示为i和p之间的距离。
[0121]
因为离群因子越大，数据点是异常点的概率越高，所以本发明将离群因子定义为稀有度分数，稀有度分数越高，细胞属于稀有类型的概率越大。
[0122]
本发明实施例可以通过两种参数选择方式输出lof稀有性分数。一种方式是输入默认稀有细胞比例0.02，另一种方式是本发明设计的自适应参数优化方式。它是通过以数据驱动的方式输入默认稀有细胞比例0.02，然后将输出所有细胞的连续稀有性分数并将其从大到小进行排序，绘制散点图。假设wi表示细胞i的稀有性分数，其中i∈(1,2,
…
,n)。那么
[0123]
δwi＝w
i-w
i 1
[0124]
上式中δwi表示细胞i和其相邻细胞之间的稀有性分数差值。稀有性分数变化率即二阶偏导数表示为：
[0125][0126]
其中δδwi的极值点位置的横坐标即稀有度分数排名靠前的细胞个数，本发明将其设定为lof方法中的异常值比例，使模型达到最佳实验效果。
[0127]
步骤4，将稀有类型细胞的二元预测结果结合社区发现算法进一步对稀有细胞进行亚群划分。
[0128]
具体实施时，在得到稀有类型细胞的二元预测结果后，可以将稀有类型细胞的二
元预测结果结合社区发现算法进一步对稀有细胞进行亚群划分。
[0129]
本实施例提供的基于拓扑特性的稀有类型细胞检测方法，通过从稀有类型细胞的近邻差异出发构建拓扑特征矩阵，其优点是有效降低了数据维度，提升了预测效率、精准度和适应性。
[0130]
下面将结合一个实施例对本方案进行说明，为了验证本发明方法的有效性，我们在r包splatter产生的12个模拟数据集上和5套标准单细胞测序数据集baron_human、shekhar、hcf-spleen、macosko和pbmc～68k进行了测试，数据集的详细信息被列在表1中。并和giniclust3、fire、cellsius、gapclust和raceid3五种稀有细胞检测方法进行比较分析。本发明引用fire方法的模拟数据集，即采用r包splatter混合jurkat和293t数据集。按照类别不变，即只分为常见细胞类型和稀有细胞类型两类，控制常见类型细胞数量为1500，调整稀有类型细胞占比为0.5％、1％、1.5％、2％、2.5％、5％、10％、15％，模拟生成了8组数据集。这8组数据集按照细胞个数分别命名为simulate_1508、simulate_1516、simulate_1523、simulate_1531、simulate_1539、simulate_1579、simulate_1667和simulate_1765，其中稀有类型细胞数量分别为8、16、23、31、39、79、167和265。本发明在这8组模拟数据集上比较了topo方法和giniclust3、fire、cellsius、gapclust和raceid3方法的recall、precision和f1值。因常见细胞类型数目需要发生变化，采用单一的真实数据集无法完成本实验的数据模拟。所以不同于之前通过真实数据集进行混合生成模拟数据集的方式，本实验采用splatter随机生成细胞数目为5000，基因数目为4000的模拟数据集。在模拟数据集中，固定稀有细胞比例为1％，调整其中的常见细胞类型个数为1，4，7和10，即总细胞类型数量为2，5，8和11共计4套模拟数据集。本发明将这4套模拟数据集按照细胞类型数量命名为simulate_2clusters、simulate_5clusters、simulate_8clusters和simulate_11clusters。因raceid3方法当数据集细胞数目达到3000时，其计算消耗时间超过48小时，所以本发明在这8组模拟数据集上比较了topo方法和giniclust3、fire、cellsius和gapclust方法的recall、precision和f1值。为了分析topo框架对于稀有类型细胞亚型划分的实验效果，本发明设计了基于baron_human、shekhar、hcf-spleen、macosko和pbmc～68k共5套数据集的15个稀有细胞类簇进行对比分析实验。
[0131]
数据集细胞数量基因数量稀有细胞数量稀有比例组织来源baron_human8569201251181.3％the pancreasshekhar26830131661390.5％the retinahcf-spleen13180141331891.4％the spleenmacosko44808232882690.6％the retinapbmc～68k68579327383740.5％a healthy donor
[0132]
表1
[0133]
为客观评价本方法的有效性，选用两个常用指标recall、precision和f1值三个指标作为评价标准。常规的二分类预测模型通常将结果分为正类和负类，并根据真实的正类和负类及预测的正类和负类构建混淆矩阵(表2)，其中tp(true positive)是指将真实的正样本正确预测为正，tn(true negative)是指将真实的负样本正确预测为负，fp(false positive)是指将真实的负样本错误预测为正，fn(false negative)是指将真实的正样本错误的预测为负。
[0134][0135]
表2
[0136]
在二分类模型中，召回率recall定义为：
[0137][0138]
准确率precision定义为：
[0139][0140]
综合评价指标f1值由recall和precision计算得到，其定义为：
[0141][0142]
实验中topo方法默认参数设置见表3。
[0143]
参数默认取值特征数2000距离计算方式cosine近邻数量k1000稀有细胞比例0.02
[0144]
表3
[0145]
实验结果表明，本发明的方法在不同样本下，基于recall、precision及f1值均优于其他的方法。由此可见我们提出的方法topo框架可以精准预测稀有细胞并实现稀有细胞类型的有效异质性分析。
[0146][0147]
表4
[0148][0149]
表5
[0150][0151]
表6
[0152]
表4、表5和表6为稀有细胞类型占数据集中不同比例时，topo方法与其他稀有细胞检测方法的recall、precision和f1值性能表现，其中，表4表示本发明的方法和5个比较方法在八套不同稀有细胞比例的模拟数据集上的recall值，表5表示本发明的方法和5个比较方法在八套不同稀有细胞比例的模拟数据集上的precision值，表6表示本发明的方法和5个比较方法在八套不同稀有细胞比例的模拟数据集上的f1值。例如simulate_1508表示模拟数据集共包含1508个细胞，其中8个为稀有细胞，其稀有度为0.5％。从结果可以看出，在fire工具模拟的数据集中，topo方法的性能不仅不会受到稀有细胞占比变化的影响，并且其recall、precision和f1值都可以稳定保持在100％。简单来说，数据集真实标签中的所有稀有细胞在预测结果中全部被检测出，并且没有一个预测错误，目前没有任何稀有细胞检测方法可以达到这样的性能。同时实验结果也验证了topo方法的性能表现并不会受到数据集中稀有细胞比例变化的影响。
[0153][0154]
表7
[0155][0156]
表8
[0157][0158][0159]
表9
[0160]
表7、表8和表9表示当稀有细胞比例不变，数据集中常见细胞类型种类增加时topo方法的性能表现。实验结果表明topo方法的稀有细胞检测性能并不受常见细胞类型数量的影响，即无论常见细胞类型数量为多少，都不影响topo方法的实验性能，并且与上个模拟实验结果相同，在这组模拟实验中topo方法的recall、precision和f1值仍然稳定保持在100％，并且无假阳率。
[0161][0162]
表10
[0163][0164]
表11
[0165][0166][0167]
表12
[0168]
表10、表11和表12是topo方法和giniclust3、fire、cellsius、gapclust共4种方法在五套真实数据集上的recall、precision和f1值性能表现。实验结果表明，topo方法相比于其他稀有细胞检测方法的预测结果有显著提升，topo方法的预测性能远优于其他方法。
[0169]
图3是本发明结合了seurat工具的louvain方法进一步对稀有细胞类型进行分类，并比较了单独采用seurat工具和toomanycells工具在真实数据集baron_human、shekhar、hcf-spleen、macosko和pbmc～68k聚类时稀有细胞亚型划分的结果。其中baron_human数据集的稀有细胞类型有macrophage、mast、epsilon、schwann和t_cell，shekahr数据集的稀有细胞类型有rods和cones，hcf-spleen数据集的稀有细胞类型有megakaryocytes和afp_albpositivecells，macosko数据集的稀有细胞类型有fibroblasts、microglia、pericytes
和astrocytes，pbmc～68k数据集的稀有细胞类型有cd34 和cd4 t，以上所有数据集的稀有细胞类型共计15个。通过对于这15个稀有细胞类型的聚类结果进行分析，在所有稀有细胞类型上，采用topo框架进行亚型划分的f1值都远高于直接采用经典聚类工具seurat和toomanycells的结果，甚至对于大多数稀有细胞类型，topo框架预测结果的f1值都可以达到100％，这表明topo框架可以可靠准确的实现稀有类型细胞的检测正确率和亚型划分正确率。
[0170]
与上面的方法实施例相对应，参见图4，本公开实施例还提供了一种基于拓扑特性的稀有类型细胞检测系统40，包括：
[0171]
预处理模块400，用于输入一个单细胞表达谱、n个细胞和m个基因，通过数据预处理保留表达散度排名靠前的m1个高可变基因构成特征矩阵，并对得到的特征矩阵计算knn最近邻距离矩阵x；
[0172]
计算模块402，用于对得到n个细胞与k个近邻细胞距离构成的近邻距离矩阵x，根据每个细胞与其他k个近邻细胞的距离向量，计算其熵值h(x)、平均值中位数m
0.5
和偏度g作为该细胞的四个领域拓扑特性，并以此构建细胞领域拓扑特征矩阵r；
[0173]
预测模块403，用于基于拓扑特征矩阵r，结合局部异常因子检测算法计算每个细胞的稀有性分数，并根据自适应参数优化方法设定阈值，输出稀有类型细胞的二元预测结果；
[0174]
划分模块404，用于将稀有类型细胞的二元预测结果结合社区发现算法进一步对稀有细胞进行亚群划分。
[0175]
图4所示系统可以对应的执行上述方法实施例中的内容，本实施例未详细描述的部分，参照上述方法实施例中记载的内容，在此不再赘述。
[0176]
参见图5，本公开实施例还提供了一种电子设备50，该电子设备包括：至少一个处理器以及与该至少一个处理器通信连接的存储器。其中，该存储器存储有可被该至少一个处理器执行的指令，该指令被该至少一个处理器执行，以使该至少一个处理器能够执行前述方法实施例中的基于拓扑特性的稀有类型细胞检测方法。
[0177]
本公开实施例还提供了一种非暂态计算机可读存储介质，该非暂态计算机可读存储介质存储计算机指令，该计算机指令用于使该计算机执行前述方法实施例中的基于拓扑特性的稀有类型细胞检测方法。
[0178]
本公开实施例还提供了一种计算机程序产品，该计算机程序产品包括存储在非暂态计算机可读存储介质上的计算程序，该计算机程序包括程序指令，当该程序指令被计算机执行时，使该计算机执行前述方法实施例中的基于拓扑特性的稀有类型细胞检测方法。
[0179]
下面参考图5，其示出了适于用来实现本公开实施例的电子设备50的结构示意图。本公开实施例中的电子设备可以包括但不限于诸如移动电话、笔记本电脑、数字广播接收器、pda(个人数字助理)、pad(平板电脑)、pmp(便携式多媒体播放器)、车载终端(例如车载导航终端)等等的移动终端以及诸如数字tv、台式计算机等等的固定终端。图5示出的电子设备仅仅是一个示例，不应对本公开实施例的功能和使用范围带来任何限制。
[0180]
如图5所示，电子设备50可以包括处理装置(例如中央处理器、图形处理器等)501，其可以根据存储在只读存储器(rom)502中的程序或者从存储装置508加载到随机访问存储器(ram)503中的程序而执行各种适当的动作和处理。在ram 503中，还存储有电子设备50操
作所需的各种程序和数据。处理装置501、rom 502以及ram 503通过总线504彼此相连。输入/输出(i/o)接口505也连接至总线504。
[0181]
通常，以下装置可以连接至i/o接口505：包括例如触摸屏、触摸板、键盘、鼠标、图像传感器、麦克风、加速度计、陀螺仪等的输入装置506；包括例如液晶显示器(lcd)、扬声器、振动器等的输出装置507；包括例如磁带、硬盘等的存储装置508；以及通信装置509。通信装置509可以允许电子设备50与其他设备进行无线或有线通信以交换数据。虽然图中示出了具有各种装置的电子设备50，但是应理解的是，并不要求实施或具备所有示出的装置。可以替代地实施或具备更多或更少的装置。
[0182]
特别地，根据本公开的实施例，上文参考流程图描述的过程可以被实现为计算机软件程序。例如，本公开的实施例包括一种计算机程序产品，其包括承载在计算机可读介质上的计算机程序，该计算机程序包含用于执行流程图所示的方法的程序代码。在这样的实施例中，该计算机程序可以通过通信装置509从网络上被下载和安装，或者从存储装置508被安装，或者从rom 502被安装。在该计算机程序被处理装置501执行时，执行本公开实施例的方法中限定的上述功能。
[0183]
需要说明的是，本公开上述的计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质或者是上述两者的任意组合。计算机可读存储介质例如可以是——但不限于——电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。计算机可读存储介质的更具体的例子可以包括但不限于：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机访问存储器(ram)、只读存储器(rom)、可擦式可编程只读存储器(eprom或闪存)、光纤、便携式紧凑磁盘只读存储器(cd-rom)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本公开中，计算机可读存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。而在本公开中，计算机可读信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读信号介质还可以是计算机可读存储介质以外的任何计算机可读介质，该计算机可读信号介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括但不限于：电线、光缆、rf(射频)等等，或者上述的任意合适的组合。
[0184]
上述计算机可读介质可以是上述电子设备中所包含的；也可以是单独存在，而未装配入该电子设备中。
[0185]
上述计算机可读介质承载有一个或者多个程序，当上述一个或者多个程序被该电子设备执行时，使得该电子设备可以执行上述方法实施例的相关步骤。
[0186]
或者，上述计算机可读介质承载有一个或者多个程序，当上述一个或者多个程序被该电子设备执行时，使得该电子设备可以执行上述方法实施例的相关步骤。
[0187]
可以以一种或多种程序设计语言或其组合来编写用于执行本公开的操作的计算机程序代码，上述程序设计语言包括面向对象的程序设计语言—诸如java、smalltalk、c ，还包括常规的过程式程序设计语言—诸如“c”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、
部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(lan)或广域网(wan)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
[0188]
附图中的流程图和框图，图示了按照本公开各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分，该模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连地表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
[0189]
描述于本公开实施例中所涉及到的单元可以通过软件的方式实现，也可以通过硬件的方式来实现。
[0190]
应当理解，本公开的各部分可以用硬件、软件、固件或它们的组合来实现。
[0191]
以上所述，仅为本公开的具体实施方式，但本公开的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本公开揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本公开的保护范围之内。因此，本公开的保护范围应以权利要求的保护范围为准。