CN120319326B

CN120319326B - 基于元细胞的空间转录组细胞类型注释方法、设备和介质

Info

Publication number: CN120319326B
Application number: CN202510814919.2A
Authority: CN
Inventors: 陈哲名; 车天栋; 乔亚娟; 张子豪
Original assignee: Hangzhou Lianchuan Biotechnology Co ltd
Current assignee: Hangzhou Lianchuan Biotechnology Co ltd
Priority date: 2025-06-18
Filing date: 2025-06-18
Publication date: 2025-10-10
Anticipated expiration: 2045-06-18
Also published as: CN121171373A; CN120319326A

Abstract

本申请公开了一种基于元细胞的空间转录组细胞类型注释方法、设备和介质，属于生物信息技术领域。所述方法包括以下步骤：将具有细胞类型注释的单细胞组学数据按细胞类型拆分，将同一细胞类型的细胞聚类并生成元细胞，获得元细胞特征计数矩阵；计算各细胞类型相对于其他细胞类型特异性的标志特征；计算各细胞类型标志特征的基矩阵，并基于待注释的空间组学数据的特征计数矩阵，计算每个空间位置对应的细胞类型权重值，取权重值最高的细胞类型作为该空间位置的细胞类型注释。利用本申请的方法、设备或介质进行空间细胞类型注释，大幅度减少了分析过程的计算量，提高了分析效率，增加了注释的准确性和稳定性。

Description

基于元细胞的空间转录组细胞类型注释方法、设备和介质

技术领域

本申请涉及生物信息技术领域，尤其涉及一种基于元细胞的空间转录组细胞类型注释方法、设备和介质。

背景技术

近年来，空间组学和单细胞组学测序技术快速发展和普及。单细胞组学技术能够以单个细胞的分辨率揭示基因表达或染色质开放程度、蛋白质峰度等的异质性，同时保留较丰富的基因表达或其他信息，但需通过组织解离等手段分离出单个细胞，导致细胞的空间位置信息完全丢失。而空间组学技术则可以保留细胞原位信息，将基因表达或蛋白等与空间坐标结合，包含了基于显微解剖、空间条形码或成像的各种方法（如10x Visium、MERFISH、Stereo-seq、10x Visium HD等产品）。然而，空间组学存在分辨率低或者检测深度不足等缺陷，导致其分析过程中细胞鉴定环节往往需要使用单细胞组学数据辅助鉴定。以转录组为例，当前使用单细胞转录组数据辅助注释空间转录组细胞类型的方法有多种，但是通常使用完整的单细胞转录组数据计算，而单细胞数据通常有数万、数十万乃至上千万个细胞，进而导致分析计算量巨大、耗时漫长；或者对单细胞数据进行随机抽样，使用少量细胞进行分析，但常规的单细胞转录组表达谱为稀疏矩阵，且可能存在较大噪声，这会导致分析结果容易受到随机因素（如噪声或未捕获低表达的标志性基因）影响，特别是在低质量数据中容易产生较大误差。

为应对大规模单细胞测序数据带来的计算开销，研究人员提出了多种高效的单细胞分析工具，主要用于数据插补、整合、聚类和细胞类型注释等任务。然而，这些工具通常专门针对特定任务设计，难以直接集成到现有的单细胞数据分析框架中。为实现更通用且高效的单细胞数据处理，一种解决方案是对原始数据进行压缩，进而降低数据冗余，赋能传统分析工具以更高效地处理大规模测序数据。针对单细胞数据压缩，一种代表性方法是元细胞（Metacell）推断，其通过聚合生物学上相似的细胞群体，将若干个单细胞压缩为单个代表性的元细胞，从而在最大程度上保留生物信息的情况下有效减少了细胞数量。

元细胞推断方法在大规模数据处理中具有显著优势。一方面，元细胞带来的数据压缩减少了测序数据分析的计算开销；其次，通过聚合特征相似的细胞，元细胞缓解了数据的稀疏性，在一定程度上提升了下游分析（如细胞类型注释、发育轨迹推断等）的鲁棒性。然而，尽管元细胞推断方法在一些应用场景下取得了令人满意的效果，其在大规模数据集上的准确性和效率仍存在不足。例如，当前最优的SEACell算法通过构建全局的邻接矩阵进行单细胞聚类，并根据聚类结果来推断元细胞。该算法在处理较小规模数据时取得了良好的效果，但当处理10万个单细胞时则需要大于一天的时间开销，且由于其指数级的内存开销，难以处理更大规模的单细胞数据。换而言之，现有的元细胞推断方法本质上将计算瓶颈从下游分析转移到了元细胞推理阶段，而未真正解决计算复杂度问题。

发明内容

为了解决上述技术问题中的至少一个，本申请采用的技术方案如下。

本申请第一方面提供一种基于元细胞的空间组学数据细胞类型注释方法，包括以下步骤：

S101，将具有细胞类型注释的单细胞组学数据按细胞类型拆分，将同一细胞类型的细胞聚类并生成元细胞，获得元细胞特征计数矩阵；

S102，基于所述元细胞特征计数矩阵，计算各细胞类型相对于其他细胞类型特异性的标志特征，将所有细胞类型的标志特征合并生成标志特征集；

S103，基于所述元细胞特征计数矩阵，计算各细胞类型标志特征的基矩阵；

S104，基于待注释的空间组学数据的特征计数矩阵，利用所述基矩阵，计算每个空间位置对应的细胞类型权重值，取权重值最高的细胞类型作为该空间位置的细胞类型注释。

元细胞是一种通过聚类单细胞数据生成“超细胞”单元的计算策略，旨在通过聚合相似细胞，降低单细胞数据的高噪声和稀疏性，同时保留生物学异质性。

在本申请的一些实施方案中，步骤S101通过以下步骤生成元细胞特征计数矩阵：

S1011，对于同一细胞类型的单细胞组学数据，对细胞-特征计数矩阵进行降维处理，得到细胞降维矩阵；

S1012，基于所述细胞降维矩阵，计算矩阵中任意两个细胞间的距离，并基于细胞间的距离，将细胞聚类，进一步计算每个细胞距离最小即最相似的k个细胞，其中k=10~30；

S1013，随机选取N个细胞，所述N个细胞互不重叠并且在所有聚类中均有分布，其中N=100~1000，对于所述N个细胞中的任意一个细胞的任意一个特征，计算该特征在该细胞及与该细胞最相似的k个细胞中的平均计数，由此得到N个特征平均计数矩阵。

对于每个细胞类型，均重复步骤S1011~S1013，由此得到所述元细胞特征计数矩阵。

在本申请的一些实施方案中，利用主成分分析（PCA）对所述细胞-特征计数矩阵进行降维处理。

在本申请的一些具体实施方案中，将细胞聚成10类。

在本申请中，N的取值根据细胞数量多少进行调整，通常取该细胞类型的细胞数量的1%；但当该细胞类型细胞数量少于10000时，N设为100；当该细胞类型细胞数量大于100000时，N设为1000。

在本申请的一些实施方案中，对于所述N个细胞中的任意一个细胞，利用以下公式计算该细胞中各特征的平均计数：

其中，表示该细胞中第i个特征的平均计数；表示该细胞中第i个特征的计数，表示该细胞中第j个最相似细胞中第i个特征的计数，j=1~k，由此得到N个特征平均表达计数矩阵。

在本申请的一些可实施方案中，步骤S102中，计算所述元细胞特征计数矩阵中每种细胞类型每个特征相对于其他细胞类型的平均差异和p值，在每个细胞类型中，将特征按p值从小到大进行排序，对于p值相同的特征，则按平均差异从大到小排序，每个细胞类型仅保留前10个特征作为标志特征。

在本申请的一些实施方案中，平均差异表示该特征在该细胞类型中的计数相对于在其他细胞类型中的计数的差异大小的平均值，利用差异大小log₂FC表征，即差异倍数（FC）经Log₂处理后的值，由此平均差异利用avg_log₂FC表示。

在本申请的一些可实施方案中，步骤S103具体如下：

将所述元细胞特征计数矩阵中所有非标志特征集中的特征剔除，得到元细胞标志特征矩阵，并分解得到所述基矩阵和第一系数矩阵。

其中，，m表示标志特征数量，表示元细胞数量，

将分解成一个维的基矩阵W和维的第一系数矩阵，即，其中，i表示细胞类型的种类。

在本申请的一些具体实施方案中，利用非负矩阵分解算法将分解成基矩阵W和第一系数矩阵。非负矩阵分解（Non-negative Matrix Factorization，NMF）是一种无监督学习算法，旨在将一个维的非负矩阵分解为两个非负矩阵（一个维的非负矩阵和一个维的非负矩阵）的乘积。由于约束所有矩阵元素非负值，使得分解结果具有可解释性，更符合生物学场景。同时，非负矩阵分解还具有捕获数据关键模式、去除噪声、识别显著特征等作用。

在本申请的一些可实施方案中，步骤S104具体如下：

将待注释的空间组学数据的特征计数矩阵中的所有非标志特征集中的特征剔除，得到空间标志特征矩阵，基于所述空间标志特征计数矩阵利用所述基矩阵获得第二系数矩阵，利用所述第二系数矩阵表征每个空间位置上各种细胞类型的权重值。

其中，，表示空间位置数量，

在本申请的一些具体实施方案中，利用以下公式计算第二系数矩阵，

，

其中，为维的矩阵，表示每个空间位置上各种细胞类型的权重值。

在本申请的一些可实施方案中，步骤S104中，利用负最小二乘法计算第二系数矩阵。

在本申请中，所述组学数据选自表观基因组学数据、转录组学数据、蛋白质组学数据中的一种。

在本申请的一些可实施方案中，所述组学数据为转录组学数据，所述特征为基因，相应地，所述特征计数为基因表达计数。

本申请第二方面提供一种计算机设备，包括：存储器，用于存储计算机程序；处理器，用于执行所述计算机程序时实现如本申请第一方面任一所述基于元细胞的空间组学数据细胞类型注释方法的步骤。

本申请第三方面提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如本申请第一方面任一所述基于元细胞的空间组学数据细胞类型注释方法的步骤。

相对于现有技术，本申请具有以下有益效果：

本申请的方法、设备和介质，通过将单细胞聚类生成元细胞，并基于元细胞进行细胞类型特征计算和非负矩阵分解，计算空间组学中各空间位置的细胞类型权利，进一步获取各个空间位置的细胞类型注释结果，大幅度减少了分析过程的计算量，提高了分析效率，增加了注释的准确性和稳定性。

应当理解，本部分所描述的内容并非旨在标识本申请的实施例的关键或重要特征，也不用于限制本申请的范围。本申请的其它特征将通过以下的说明书而变得容易理解。

附图说明

通过参考附图阅读下文的详细描述，本申请示例性实施方式的上述以及其他目的、特征和优点将变得易于理解。在附图中，以示例性而非限制性的方式示出了本申请的若干实施方式，其中：

图1示出了本申请实施例1中基于元细胞的空间转录组学数据细胞类型注释方法的实现流程示意图；

图2示出了本申请实施例2中星型胶质细胞（astrocyte）权重值在空间位置上的分布；

图3示出了本申请实施例2中内皮细胞（endothelial cell）权重值在空间位置上的分布；

图4示出了本申请实施例2中神经元细胞（neuron）权重值在空间位置上的分布；

图5示出了本申请实施例2小胶质细胞（microglial cell）权重值在空间位置上的分布；

图6示出了本申请实施例2中少突胶质细胞（oligodendrocyte）权重值在空间位置上的分布；

图7示出了本申请实施例2中每个空间位置的细胞类型鉴定结果。

具体实施方式

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

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

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

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

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

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

实施例1 基于元细胞的空间转录组细胞类型注释方法

本实施例提供一种基于元细胞的空间转录组细胞类型注释方法，流程示意图如图1所示。

主要包括步骤S101至步骤S104，详细说明如下：

S101，将已经包含具有细胞类型注释的单细胞转录组数据按细胞类型逐个拆分，将同一细胞类型的细胞聚类并生成元细胞。通过如下步骤S1011至S1014实现：

S1011，选取同一细胞类型（如T细胞）的单细胞转录组数据，使用seurat软件的NormalizeData、FindVariableFeatures、ScaleData、RunPCA函数进行对原始的细胞-基因表达计数矩阵（如表1所示）进行归一化、标准化、高变基因计算、PCA降维处理，获得细胞-PCA降维矩阵（表2）。

表1 原始的细胞-基因表达计数矩阵（部分）

表2 细胞-PCA降维矩阵（部分）

S1012，基于细胞-PCA降维矩阵，计算矩阵中每个细胞间的欧氏距离。基于细胞间的欧氏距离，使用k-means算法将细胞聚成10类，并使用KNN（k-nearest neighbor，k最邻近）算法计算每个细胞最相似（即欧式距离最小）的k个细胞（k为超参数，可按细胞数调整，设定为大于1的任意整数，通常设置为25）。

S1013，生成元细胞。设置预期的元细胞数量N（N取值100~1000，取该细胞类型的细胞数量的1%，但当该细胞类型细胞数量少于10000时，N设为100；当该细胞类型细胞数量大于100000时，N设为1000），在矩阵中随机选取N个细胞（N个细胞互不重叠，并须在步骤S1012获得的10个聚类中均有分布）。在N个细胞中，逐一将每个细胞与各自最相似的k个细胞的原始基因表达计数矩阵合并，如表3所示。

表3 随机选中的细胞及其最邻近的k个细胞的原始表达谱（部分）

利用以下公式计算随机选中的细胞中各基因的平均表达计数：

其中，表示随机选中的细胞中第i个基因的平均表达计数；表示随机选中的细胞中第i个基因的表达计数，表示随机选中的细胞第j个最相似细胞中第i个基因的表达计数，j=1~k。

由此，N个细胞生成N个基因平均表达计数矩阵，每个基因平均表达计数矩阵即为元细胞基因表达计数矩阵，即生成N个元细胞基因表达计数矩阵，如表4所示。

表4 N个元细胞基因表达计数矩阵（部分）

S1014，对每种细胞类型的单细胞转录组数据重复步骤S1011至S1013，即每个细胞类型均获得N个元细胞基因表达计数矩阵，合并之后则获得本实施例单细胞转录组数据对应的完整的元细胞基因表达计数矩阵，即元细胞count矩阵。

S102，对元细胞count矩阵进行归一化和测序深度矫正后，计算各细胞类型的标志基因。

使用seurat软件的NormalizeData函数对矩阵进行归一化和测序深度矫正，将处理后的矩阵称为元细胞data矩阵。使用seurat软件的FindAllMarkers函数和Wilcox差异检验方法计算元细胞data矩阵中每种细胞类型每个基因的avg_log₂FC（log₂处理后的差异倍数）和p值（即显著性水平）——本质上是逐一将单个细胞类型的细胞与剩余所有细胞类型的细胞的基因进行Wilcox秩和检验，以获得每个细胞类型特异性高表达的基因。在每个细胞类型中，将基因按p值从小到大进行排序，对于p值相同的基因，则按avg_log₂FC从大到小排序。每个细胞类型仅保留前10个标志基因（如表5所示），将所有细胞类型的标志基因合并生成标志基因集。

表5 各细胞类型TOP10标志基因（部分）

S103，剔除元细胞data矩阵中所有非标志基因集的基因，称为元细胞标志基因data矩阵，即一个维的矩阵，其中m表示标志基因数量，表示元细胞数量。使用非负矩阵分解算法将矩阵分解成一个维的基矩阵W和维的系数矩阵，即。其中，i表示细胞类型的种类。

S104，将空间转录组数据的空间位置-基因表达计数矩阵用seurat软件的NormalizeData函数进行归一化和测序深度矫正处理，同样剔除所有非标志基因集的基因，生成空间标志基因data矩阵，即一个维的矩阵，其中表示空间位置数量。使用非负最小二乘法，计算中的系数矩阵。其中为维的矩阵，表示每个空间位置上各种细胞类型的权重值。对于每个空间位置，将权重值最高的细胞类型作为该空间位置的细胞类型注释。

实施例2 基于元细胞的空间转录组细胞类型注释方法的应用

本实施例提供一种基于元细胞的空间转录组细胞类型注释方法的应用。使用单细胞转录组数据对小鼠全脑的高清空间转录组数据进行细胞类型注释，具体方法如实施例1所示。

空间转录组原始数据如表6所示，随后使用seurat软件的NormalizeData函数进行归一化和测序深度矫正处理，同样剔除所有非标志基因集的基因，生成空间标志基因data矩阵，即矩阵。如步骤S104所示，使用非负最小二乘法计算每个空间位置上各种细胞类型的权重值（如表7和图2~图6），并将权重值最高的细胞类型作为该空间位置的细胞类型注释（如表8和图7）。

表6 空间转录组原始的空间位置-基因表达计数矩阵

表7 每个空间位置上各种细胞类型的权重值

表8 每个空间位置的细胞类型的鉴定结果

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

Claims

1.一种基于元细胞的空间组学数据细胞类型注释方法，其特征在于，包括以下步骤：

S104，基于待注释的空间组学数据的特征计数矩阵，利用所述基矩阵，计算每个空间位置对应的细胞类型权重值，取权重值最高的细胞类型作为该空间位置的细胞类型注释，

其中，步骤S103具体如下：

利用非负矩阵分解算法将所述元细胞特征计数矩阵中所有非标志特征集中的特征剔除，得到元细胞标志特征矩阵，，m表示标志特征数量，表示元细胞数量，

将分解成一个维的基矩阵W和维的第一系数矩阵，即，其中，i表示细胞类型的种类，

步骤S104具体如下：

将待注释的空间组学数据的特征计数矩阵中的所有非标志特征集中的特征剔除，得到空间标志特征计数矩阵，其中，，表示空间位置数量，

基于所述空间标志特征计数矩阵，利用以下公式计算第二系数矩阵，

，

2.根据权利要求1所述的一种基于元细胞的空间组学数据细胞类型注释方法，其特征在于，步骤S101通过以下步骤实现：

S1013，随机选取N个细胞，所述N个细胞互不重叠并且在所有聚类中均有分布，其中N=100~1000，对于所述N个细胞中的任意一个细胞的任意一个特征，计算该特征在该细胞及与该细胞最相似的k个细胞中的平均计数，由此得到N个特征平均计数矩阵，

3.根据权利要求1所述的一种基于元细胞的空间组学数据细胞类型注释方法，其特征在于，步骤S102中，计算所述元细胞特征计数矩阵中每种细胞类型每个特征相对于其他细胞类型的平均倍数和p值，在每个细胞类型中，将特征按p值从小到大进行排序，对于p值相同的特征，则按平均差异从大到小排序，每个细胞类型仅保留前10个特征作为标志特征。

4.根据权利要求1~3任一所述的一种基于元细胞的空间组学数据细胞类型注释方法，其特征在于，所述组学数据选自表观基因组学数据、转录组学数据、蛋白质组学数据中的一种。

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

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

处理器，用于执行所述计算机程序时实现如权利要求1-4任一所述的方法的步骤。

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

所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1-4任一所述的方法的步骤。