CN112837754A - 一种基于特征基因的单细胞自动分类方法和装置 - Google Patents

Abstract

本发明提供了一种基于特征基因的单细胞自动分类方法和装置，所述方法包括：将单细胞基因表达矩阵进行数据过滤、标准化和选择高变异基因后，进行降维处理和聚类处理，获得一个或多个细胞类群；采用负二项分布，从细胞类群中筛选特征基因集；将筛选的特征基因集与细胞类型数据库进行比对，采用超几何分布对比对结果进行检验统计，完成单细胞类型注释和单细胞自动分类。本发明的基于特征基因的单细胞自动分类方法AngryCell涵盖数据预处理、降维聚类、基于负二项分布的特征基因筛选、基于超几何分布的细胞类型匹配以及下游富集分析以及结果可视化等部分，获得的结果可靠准确、可解释性强。

Description

一种基于特征基因的单细胞自动分类方法和装置

技术领域

本发明属于单细胞测序技术领域，涉及一种基于特征基因的单细胞自动分类方法和装置。

背景技术

目前，主流的单细胞转录组数据分析流程的核心步骤包括：获取单细胞基因表达定量矩阵并进行降维处理，将同一批实验测到的所有细胞分为数个细胞类群，这一步骤即聚类(Clustering)；接下来将各个细胞类群互相比较，找到每个细胞类群特征性高表达的基因，称为特征基因(Marker)；再根据特征基因，预测每个细胞类群所属的细胞类型，并对每个细胞类群进行标注(Cluster annotation)，标注的细胞数据用于科学问题的深度分析，这种对单细胞类群进行标注的过程称为单细胞分类。

单细胞分类在大部分单细胞数据分析过程中由人工完成。数据分析人员在获取特征基因后，通过查询数据库、文献或根据经验，对每个细胞类群的归属进行预测和判断。这一预测过程完全由人为判断，既耗时耗力，又受到操作者的经验限制，有时会出现错误。另外，由于人工行为不具有可重复性，往往导致结果的科学可靠性受到质疑。因此，近年来国际上逐步建立了单细胞数据的自动分类方法。已有的单细胞自动分类方法有两种思路：一种基于参考数据，一种基于参考基因集。

基于参考数据的自动分类方法包括SingleR(Aran D,Looney AP,Liu L,etal.Reference-based analysis of lung single-cell sequencing reveals atransitional profibrotic macrophage[J].Nature immunology,2019,20(2):163-172.)、SCMAP(Kiselev V Y,Yiu A,Hemberg M.scmap:projection of single-cell RNA-seq data across data sets[J].Nature methods,2018,15(5):359-362.)和SciBet(LiC,Liu B,Kang B,et al.SciBet as a portable and fast single cell typeidentifier[J].Nature communications,2020,11(1):1-8.)等。该方法使用人工标注的单细胞数据作为参考，进行模型拟合，再利用拟合的模型对新的单细胞数据进行预测，可以归类为使用机器学习的方法。基于参考数据的自动分类方法既可以对单个细胞进行预测，也可以对聚类产生的细胞类群进行预测。

基于参考基因集的自动分类方法目前仅有Garnett(Pliner H A,Shendure J,Trapnell C.Supervised classification enables rapid annotation of cell atlases[J].Nature methods,2019,16(10):983-986.)一种。它的原理是使用人工定义的各种细胞类型的特征基因集，先在待检测的单细胞数据中寻找符合特征基因定义的阳性细胞，再使用阳性细胞拟合模型，对剩下的单细胞数据进行预测，这一方法也归类于机器学习方法。

然而，基于参考数据的自动分类方法的最大问题是不能完全脱离人工监督。机器学习方法需要根据合适的已经进行人工鉴定的数据，才能够对新数据做出预测判断，如果没有已经标注好的单细胞数据，该方法将不能使用。对于单细胞测序方兴未艾的今天，很多单细胞数据来源于从未研究过的疾病和组织，研究者仅仅知道其中可能存在的细胞类型，没有可以用作拟合数据的单细胞数据集；其次，模型拟合方法给出的结果是否准确，完全取决于参考数据的准确性，如果标注的数据质量差或存在错误，则新数据的预测结果将存在错误；再次，该方法只能拟合参考数据里已经覆盖的细胞类型，对于参考数据里不存在的细胞类型，该方法无法学习其特征，鉴定出的细胞类型会较参考数据中的细胞类型更少。目前，单细胞测序研究的重要目的是发现新的细胞类群，而现有的基于参考数据的自动分类方法不能完全满足这一需要。

基于参考基因集的自动分类方法，Garnett的问题是使用了复杂的网络模型设置细胞身份，计算缓慢，效果也不甚理想。同时，该方法需要使用者手工指定细胞类型的参考基因列表，尽管有一些数据库可供查询，但是很难得到一个高质量的参考基因列表。在两次分析中，如果不同的分析人员输入了不同的参考基因列表，计算结果可能截然不同。因此，这一方法受到操作者人为因素的影响较大。

以上方法都使用了机器学习的核心思想，模型是拟合得出的，必然无法脱离参考数据，且存在一定程度的难以解释的成分。对每个自动判断的细胞类群，既不输出其他潜在的候选类型，也不输出判断依据，结果的可解释性较差。总而言之，现有的自动化分析方法无法完全复现分析人员在进行细胞标注时的思考过程，结果的可利用性较低。

发明内容

针对现有技术的不足和实际需求，本发明提供了一种基于特征基因的单细胞自动分类方法和装置，所述方法基于特征基因进行单细胞自动分类，避免了参考数据的限制性，模拟分析人员的决策过程，给出判断依据，结果的可解释性强，分类标准统一，结果的准确性和可靠性高。

为达此目的，本发明采用以下技术方案：

第一方面，本发明提供了一种基于特征基因的单细胞自动分类方法，所述方法包括：

将单细胞基因表达矩阵进行数据过滤、标准化和选择高变异基因后，进行降维处理和聚类处理，获得一个或多个细胞类群；

采用负二项分布，从细胞类群中筛选特征基因集；

将筛选的特征基因集与细胞类型数据库进行比对，采用超几何分布对比对结果进行检验统计，完成单细胞类型注释和单细胞自动分类。

本发明中，基于特征基因的单细胞自动分类方法AngryCell将单细胞基因表达矩阵进行预处理、降维聚类处理后，采用基于负二项分布的特征基因筛选方法和基于超几何分布的细胞类型匹配方法，不依赖于参考数据和机器学习思想，实现了基于特征基因的单细胞自动分类，各个步骤各司其职、承上启下、相辅相成，保证了分类结果的准确性和可解释性。

优选地，所述基于特征基因的单细胞自动分类方法AngryCell首先将单细胞基因表达矩阵和细胞样本信息作为输入数据，过滤去除质量不合格的细胞和基因，并进行数据标准化，移除测序深度差异的影响；随后基于方差值的排序，选择具体参数(例如，前2000个基因)作为高变异基因；

采用传统的线性降维方法，将基于高变基因X细胞的零均值化表达矩阵进行降维处理，例如采用奇异值分解方法或主成分分析(PCA)获取低维(<100d)空间信息，引用Seurat基于社区划分的Louvain算法进行聚类处理，获得一个或多个细胞类群；

从细胞类群中筛选特征基因集，并与细胞类型数据库进行比对，完成单细胞类型注释和单细胞自动分类。

优选地，所述特征基因集的筛选方法根据单细胞基因表达量平均数与基因零表达比例呈负二项分布的特性，确定每个细胞类群的零表达比例阈值，根据零表达比例阈值筛选特征基因集，包括：

根据细胞类群的基因表达数据获取基因表达量平均数μ和基因表达量方差Var，将所述基因表达量平均数μ和基因表达量方差Var进行数学拟合，获得细胞类群的负二项分布参数

根据细胞类群的基因表达量平均数μ和负二项分布参数

计算负二项分布拟合的零表达比例

并根据初始零表达比例α确定零表达比例阈值cutoff；

根据零表达比例阈值cutoff筛选细胞类群的特征基因集。

优选地，所述基因表达量平均数和基因表达量方差进行数学拟合的公式为：

其中，μ为基因表达量平均数，

为负二项分布参数，Var为基因表达量方差；

所述负二项分布拟合的零表达比例的计算公式为：

其中，μ为基因表达量平均数，

为负二项分布参数，

为负二项分布拟合的零表达比例；

所述零表达比例阈值为：

其中，α为初始零表达比例，

为负二项分布参数，cutoff为零表达比例阈值。

在一个具体实施例中，AngryCell根据管家基因负二项分布拟合的零表达比例

和初始零表达比例α，筛选特征基因集。

本发明中，细胞分群特征基因集的筛选是AngryCell的关键部分，根据单细胞基因表达量平均数与基因零表达比例之间特殊的负二项分布关系，AngryCell用基因零表达比例直接推断每个细胞类群的特征基因，避免了进行细胞类群间的差异基因分析，显著节省了程序运行时间。

本发明中，细胞类群的零表达比例阈值可以设置为统一的数值，也可以根据负二项分布的标准化模型拟合计算，拟合计算的零表达比例阈值反应了细胞类群的异质性，异质性越大，拟合的零表达比例阈值越低，获取的特征基因数量也越少，标准化模型拟合的零表达比例阈值可以使异质性大的细胞类群偏向于大类的注释，提高细胞注释的准确度。

优选地，AngryCell的细胞注释基于细胞类型特征基因比对的原理进行，细胞类型特征基因收集于名为AngryCellDB的细胞类型数据库中，共分为三层注释：第一层为最大类，第二层为细胞类型，第三层为细分的亚类，用户根据数据的来源和质量情况选择不同层次的注释；AngryCell根据筛选的特征基因集和选择的数据库进行比对，并采用超几何分布对比对结果进行检验统计，获取的细胞类型匹配值P进行Benjamini-Hochberg校正，细胞类型匹配值P越小，代表预测的细胞类型越可信。

优选地，所述细胞类型匹配值P的计算公式为：

其中，N为所有基因数目，M为细胞类型数据库基因数目，n为细胞类群的特征基因数目，k为匹配的基因数目，P为细胞类型匹配值。

本发明中，AngryCell的另一关键部分是对细胞类群的注释匹配，结合特征基因集和细胞类型数据库AngryCellDB，对每个细胞类群进行快速且准确的细胞类型注释，输出数据类型丰富，包括表格、基因表达图或网络富集图，可视化和可读性高。

第二方面，本发明提供了一种基于特征基因的单细胞自动分类装置，所述装置包括：

数据预处理模块，用于将单细胞基因表达矩阵进行数据过滤和标准化；

高变异基因选择模块，用于根据数据的方差值去除噪音数据；

降维处理模块，用于对过滤和标准化的单细胞基因表达矩阵进行奇异值分解和/或主成分分析；

聚类处理模块，用于基于Louvain算法处理降维的单细胞基因表达数据；

细胞注释模块，包括特征基因集筛选单元和数据比对单元，用于基于负二项分布筛选细胞类群的特征基因集，并将筛选的特征基因集与细胞类型数据库进行比对。

优选地，所述特征基因集筛选单元用于：

根据细胞类群的基因表达数据获取基因表达量平均数和基因表达量方差，将所述基因表达量平均数和基因表达量方差进行数学拟合，获得细胞类群的负二项分布参数；

根据细胞类群的基因表达量平均数和负二项分布参数，计算负二项分布拟合的零表达比例，并根据初始零表达比例确定零表达比例阈值；

根据零表达比例阈值筛选细胞类群的特征基因集；

所述基因表达量平均数和基因表达量方差进行数学拟合的公式为：

其中，μ为基因表达量平均数，

为负二项分布参数，Var为基因表达量方差；

所述负二项分布拟合的零表达比例的计算公式为：

其中，μ为基因表达量平均数，

为负二项分布参数，

为负二项分布拟合的零表达比例；

所述零表达比例阈值为：

其中，α为初始零表达比例，

为负二项分布参数，cutoff为零表达比例阈值。

优选地，所述数据比对单元包括：

采用超几何分布计算比对结果的细胞类型匹配值，并对所述细胞类型匹配值进行Benjamini-Hochberg校正；

所述细胞类型匹配值的计算公式为：

其中，N为所有基因数目，M为细胞类型数据库基因数目，n为细胞类群的特征基因数目，k为匹配的基因数目，P为细胞类型匹配值。

与现有技术相比，本发明具有如下有益效果：

(1)本发明的基于特征基因的单细胞自动分类方法AngryCell从单细胞基因表达矩阵开始，到细胞类型注释和结果可视化，过程涵盖数据预处理、降维聚类、特征基因筛选、细胞类型匹配以及下游富集分析以及结果可视化等部分，包含多项自主研发的算法和数据架构，如基于负二项分布的特征基因选择算法和基于超几何分布的细胞类型匹配算法，各个算法各司其职、承上启下、相辅相成，保证结果的可靠性；

(2)AngryCell特征基因筛选不依赖于参考数据和机器学习思想，根据单细胞基因表达量平均数与基因零表达比例呈负二项分布关系，用基因零表达比例直接推断每个细胞类群的特征基因，避免了进行细胞类群间的差异基因分析，显著节省了程序运行时间；

(3)AngryCell细胞类型注释基于筛选的特征基因和细胞类型数据库，输出结果包含匹配细胞类型的特征基因以及文献来源等，与人工注释结果匹配度高，结果准确，可视化效果丰富多样，可解释性强。

附图说明

图1为单细胞自动分类方法AngryCell的主流程图；

图2为单细胞自动分类方法AngryCell的框图；

图3为不同类型免疫细胞的特征基因集；

图4A为PBMC的细胞类型注释结果的网络富集图，图4B为PBMC的细胞类型注释结果的基因表达图。

具体实施方式

为进一步阐述本发明所采取的技术手段及其效果，以下结合实施例和附图对本发明作进一步地说明。可以理解的是，此处所描述的具体实施方式仅仅用于解释本发明，而非对本发明的限定。

实施例中未注明具体技术或条件者，按照本领域内的文献所描述的技术或条件，或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可通过正规渠道商购获得的常规产品。

实施例1

基于特征基因的单细胞自动分类方法AngryCell的主流程图如图1所示，框图如图2所示，步骤如下：

(1)数据预处理：将单细胞基因表达矩阵和细胞样本信息作为输入数据，过滤去除质量不合格的细胞和基因，并进行数据标准化，移除测序深度差异的影响；

(2)高变异基因选择：基于方差值的排序，选择前2000个基因作为高变异基因；

(3)降维：采用传统的线性降维方法，例如奇异值分解方法或Seurat的主成分分析(PCA)方法，将基于高变基因X细胞的零均值化表达矩阵进行降维处理，获取低维(<100d)空间信息；

(4)聚类：引用Seurat基于社区划分的Louvain算法进行聚类处理，获得多个细胞类群；

(5)细胞类型注释：基于负二项分布和零表达比例的标准化模型，从细胞类群中筛选特征基因集，并与收集有公共数据库(Cellmarker、NCBI)和文献报道的细胞类型特征基因的细胞类型数据库AngryCellDB进行比对，基于超几何分布进行检验统计，完成单细胞类型注释和单细胞自动分类；

(6)输出表格、基因表达图或网络富集图等可视化数据。

其中，基于负二项分布和零表达比例的标准化模型：

基于超几何分布进行检验统计：

实施例2

本实施例以公共10×PBMC数据为测试数据，进行过滤、标准化、选择高变异基因、降维聚类处理后，获得多个细胞类群，对所述细胞类群进行细胞类型注释，步骤如下：

(1)特征基因集筛选

根据细胞类群的管家基因表达数据获取管家基因表达量平均数μ和管家基因表达量方差Var，将所述管家基因表达量平均数μ和管家基因表达量方差Var进行数学拟合，获得细胞类群的负二项分布参数

根据细胞类群的管家基因表达量平均数μ和负二项分布参数

计算负二项分布拟合的零表达比例

并根据初始零表达比例α确定零表达比例阈值cutoff；

根据零表达比例阈值cutoff筛选细胞类群的特征基因集；

如图3所示，不同类型免疫细胞的特征基因集位于负二项分布曲线的右下方；

(2)特征基因集与细胞类型数据库的比对

采用超几何分布计算比对结果的细胞类型匹配值，并对所述细胞类型匹配值进行Benjamini-Hochberg校正；

所述细胞类型匹配值的计算公式为：

其中，N为所有基因数目，M为细胞类型数据库基因数目，n为细胞类群的特征基因数目，k为匹配的基因数目，P为细胞类型匹配值，P越小，代表预测的细胞类型越可信。

表1

如表1、图4A和图4B所示为输出结果，从不同角度反映了PBMC的细胞类型注释结果，表1中，Celltype_predicted是AngryCell的自动注释结果，可以匹配人工注释的结果Orig_Idents，还展示了特征基因(Specific marker)、参考文献(Reference)等重要信息。

综上所述，本发明的基于特征基因的单细胞自动分类方法根据单细胞基因表达量平均数与基因零表达比例呈负二项分布关系，用基因零表达比例直接推断每个细胞类群的特征基因，避免了参考数据的限制性，模拟分析人员的决策过程，给出判断依据，结果的可解释性强。

申请人声明，本发明通过上述实施例来说明本发明的详细方法，但本发明并不局限于上述详细方法，即不意味着本发明必须依赖上述详细方法才能实施。所属技术领域的技术人员应该明了，对本发明的任何改进，对本发明产品各原料的等效替换及辅助成分的添加、具体方式的选择等，均落在本发明的保护范围和公开范围之内。

Claims (10)

1.一种基于特征基因的单细胞自动分类方法，其特征在于，所述方法包括：

将单细胞基因表达矩阵进行数据过滤、标准化和选择高变异基因后，进行降维处理和聚类处理，获得一个或多个细胞类群；

采用负二项分布，从细胞类群中筛选特征基因集；

将筛选的特征基因集与细胞类型数据库进行比对，采用超几何分布对比对结果进行检验统计，完成单细胞类型注释和单细胞自动分类。

2.根据权利要求1所述的方法，其特征在于，所述特征基因集的筛选方法包括：

根据细胞类群的基因表达数据获取基因表达量平均数和基因表达量方差，将所述基因表达量平均数和基因表达量方差进行数学拟合，获得细胞类群的负二项分布参数；

根据细胞类群的基因表达量平均数和负二项分布参数，计算负二项分布拟合的零表达比例，并根据初始零表达比例确定零表达比例阈值；

根据零表达比例阈值筛选细胞类群的特征基因集。

3.根据权利要求2所述的方法，其特征在于，所述基因表达量平均数和基因表达量方差进行数学拟合的公式为：

其中，μ为基因表达量平均数，

为负二项分布参数，Var为基因表达量方差；

所述负二项分布拟合的零表达比例的计算公式为：

其中，μ为基因表达量平均数，

为负二项分布参数，

为负二项分布拟合的零表达比例；

所述零表达比例阈值为：

其中，α为初始零表达比例，

为负二项分布参数，cutoff为零表达比例阈值。

4.根据权利要求1所述的方法，其特征在于，所述采用超几何分布对比对结果进行检验统计包括：

采用超几何分布计算比对结果的细胞类型匹配值，并对所述细胞类型匹配值进行Benjamini-Hochberg校正。

5.根据权利要求4所述的方法，其特征在于，所述细胞类型匹配值的计算公式为：

其中，N为所有基因数目，M为细胞类型数据库基因数目，n为细胞类群的特征基因数目，k为匹配的基因数目，P为细胞类型匹配值。

6.根据权利要求1所述的方法，其特征在于，所述降维处理包括对过滤和标准化的单细胞基因表达矩阵进行奇异值分解和/或主成分分析。

7.根据权利要求1所述的方法，其特征在于，所述聚类处理包括采用Louvain算法处理降维的单细胞基因表达数据。

8.一种基于特征基因的单细胞自动分类装置，其特征在于，所述装置包括：

数据预处理模块，用于将单细胞基因表达矩阵进行数据过滤和标准化；

高变异基因选择模块，用于根据数据的方差值去除噪音数据；

降维处理模块，用于对过滤和标准化的单细胞基因表达矩阵进行奇异值分解和/或主成分分析；

聚类处理模块，用于基于Louvain算法处理降维的单细胞基因表达数据；

细胞注释模块，包括特征基因集筛选单元和数据比对单元，用于基于负二项分布筛选细胞类群的特征基因集，并将筛选的特征基因集与细胞类型数据库进行比对。

9.根据权利要求8所述的装置，其特征在于，所述特征基因集筛选单元用于：

根据细胞类群的基因表达数据获取基因表达量平均数和基因表达量方差，将所述基因表达量平均数和基因表达量方差进行数学拟合，获得细胞类群的负二项分布参数；

根据细胞类群的基因表达量平均数和负二项分布参数，计算负二项分布拟合的零表达比例，并根据初始零表达比例确定零表达比例阈值；

根据零表达比例阈值筛选细胞类群的特征基因集；

所述基因表达量平均数和基因表达量方差进行数学拟合的公式为：

其中，μ为基因表达量平均数，

为负二项分布参数，Var为基因表达量方差；

所述负二项分布拟合的零表达比例的计算公式为：

其中，μ为基因表达量平均数，

为负二项分布参数，

为负二项分布拟合的零表达比例；

所述零表达比例阈值为：

其中，α为初始零表达比例，

为负二项分布参数，cutoff为零表达比例阈值。

10.根据权利要求8所述的装置，其特征在于，所述数据比对单元包括：

采用超几何分布计算比对结果的细胞类型匹配值，并对所述细胞类型匹配值进行Benjamini-Hochberg校正；

所述细胞类型匹配值的计算公式为：

其中，N为所有基因数目，M为细胞类型数据库基因数目，n为细胞类群的特征基因数目，k为匹配的基因数目，P为细胞类型匹配值。

Images