CN114091603A - 一种空间转录组细胞聚类、分析方法 - Google Patents

Abstract

一种空间转录组细胞聚类方法，包括步骤：对于空间转路组中每个细胞基因表达进行预处理；根据所述空间转路组的细胞坐标生成邻接矩阵A，获得空间转录组细胞的图结构表示，由细胞特征矩阵X表示细胞基因表达，将邻接矩阵A和细胞特征矩阵X输入经过训练的图卷积神经网络模型DGI；所述图卷积神经网络模型DGI输出具有空间信息的结点特征表示；对所述结点特征表示采用降维、聚类算法处理后，识别、获得所述空间转录组细胞类型。

Description

一种空间转录组细胞聚类、分析方法

技术领域

本发明涉及生物信息学领域，特别涉及一种空间转录组细胞聚类、分析方法。

背景技术

空间转录组技术是生物信息学领域近年来的重大突破。该技术可以同时测量大量细胞的空间位置和细胞内的转录组计数。通过这一方法获得的数据集包含了大量细胞的空间坐标信息，同时每个细胞都有高维属性信息。目前对该类数据的分析大都基于传统概率统计模型，难以处理海量全局空间信息和计数数据，并依赖特定模型先验。其中，细胞类型发现这一任务是众多数据分析和下游研究的必经步骤。

发明内容

本发明实施例之一，一种基于图卷积神经网络的空间转录组细胞聚类方法对于空间转路组中每个细胞基因表达进行预处理；

根据所述空间转路组的细胞坐标生成邻接矩阵A，获得空间转录组细胞的图结构表示，由细胞特征矩阵X表示细胞基因表达，

将邻接矩阵A和细胞特征矩阵X输入经过训练的图卷积神经网络模型DGI；

所述图卷积神经网络模型DGI输出具有空间信息的结点特征表示；

对所述结点特征表示采用降维、聚类算法处理后，识别、获得所述空间转录组细胞类型。

所述的对所述结点特征表示采用降维处理，是指采用PCA算法降低所述结点特征的维度。

使用K-means++算法对PCA降维后的结点进行聚类，从而发现不同的细胞类型。

附图说明

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

图1为根据本发明实施例之一空间转录组细胞聚类分析方法流程图。

图2为根据本发明实施例的一个实例的结果示意图，图中展示的是各类细胞所相关的GO term及其对应的显著性指标p值。

具体实施方式

现有细胞分类方法普遍借鉴scRNA-Seq的策略，即没有考虑空间坐标所能带来的信息。而图神经网络是近年机器学习领域的一个重要研究热点，也显示出了处理特定结构数据的强大能力，有较好的潜力应用于空间转录组数据，从而完成带有空间信息的细胞类型发现。

本发明针对现有细胞聚类方法未充分利用细胞的空间分布信息的缺陷，依赖于特定先验知识，无法满足新细胞类型发现的问题，提出一种基于图神经网络的细胞聚类方法，将机器学习技术与空间转录组领域的知识相结合，从而使用聚类的方法，发现新的细胞亚型。

根据一个或者多个实施例，一种基于图网络机器学习的针对空间转录组数据的细胞聚类方法。首先对每个细胞中的基因表达进行预处理，来得到节点特征。同时，根据细胞坐标计算得到邻接矩阵，从而将空间转录组数据转化为图结构的表示形式，每个细胞即为图上的一个结点，空间上邻近的细胞之间被一条边连接。进而将图卷积神经网络应用到这一图结构的数据上，实现各个节点的特征嵌入。之后再进行降维及聚类，最终实现带有空间结构的细胞类型发现。为了验证本方法能够发现不同特征的细胞类型，本发明分别使用了邻居丰度分析和差异表达基因分析、Gene Otology(GO)terms分析来检验不同细胞聚类在空间分布特征和生物功能上的差异。

所述的基因表达预处理，通过以下方式得到：首先根据均值排除掉表达的基因，再根据方差排除低差异的基因，然后对于从多个批次下采集到的数据，用Scanorama来去除批次效应，最后进行数据的归一化，是每个细胞中剩余的基因总表达量之和为10000。

所述的邻接矩阵是指：根据每个细胞的空间坐标信息，计算出各个细胞之间的的欧氏距离，根据距离构建邻接矩阵，具体包括：

步骤1)对于同一个批次下的细胞组，第i个细胞二维空间坐标为v_i＝(x_i,y_i)，第j个细胞的坐标为v_j＝(x_j,y_j)，则这两个细胞间的欧式距离为d_ij＝‖v_i-v_j‖，该细胞组的距离矩阵为

N为细胞的个数；

步骤2)基于上述所得的距离矩阵

通过设置合适的阈值d_thres可以得到其邻接矩阵A₀＝[a_ij]^N×N，其中：

步骤3)对邻接矩阵进行进一步调整，引入权重因子λ实现对空间信息权重的调节，调整后的邻接矩阵为A＝λ*I+(1-λ)A₀,I为一个对角矩阵：I∈R^N×N。

所述的图神经网络模型，主要使用了Deep Graph Infomax(DGI)算法的思想，DGI的基本思路是通过GCN层作为编码器来获得原始网络节点的表示向量，接下来对图网络结构随机打乱，获得打乱后的节点表示作为负样本，最后优化判别器的损失函数，使得算法可以区分正负样本是否来自原网络。DGI的目标通过最大化全局特征和局部特征之间的互信息，来学习一个编码器E来学习节点的表示。

所述的降维方法为主成分分析(PCA)降维。

所述的聚类方法为K-means++，将所有细胞向指定的K个类聚类。

所述的邻居丰度分析，首先构造一个邻居分布矩阵C＝[c_ij]^K×K，其中c_ij为第i类细胞周围的邻居中出现第j类细胞的次数。邻居的定义为：若存在一条边将两个结点直接相连，则两节点对应的细胞互为邻居细胞。之间然后对邻居分布矩阵进行归一化，

其中n_i为第i类细胞的个数。

所述的差异表达基因分析，本方法使用曼-惠特尼u检验来找到各个聚类中显著高表达的top200个基因，对于显著性指标p值小于0.5的基因进行去除。

所述的GO term分析，本方法使用的是各个聚类中筛选出的200个高表达基因，将其上传至GO term分析官方网站上，获得由高到低的一系列相关生物功能。其中背景基因选择为本数据集中的全部基因。

为了实现上述方法，涉及到的相关计算分析系统包括对空间转录组的图结构数据建立，以及建立基于DGI的空间转录组数据特征嵌入、降维及聚类和聚类结果分析模块。

本发明实施例提出了新的空间转录组数据的处理方法，即，

将图卷积神网络技术应用于细胞类型的发现，在聚类任务上完成了端到端的数据驱动的机器学习任务，在多个数据集上都有了新的细胞亚型发现。

根据一个或者多个实施例，如图1所示，一种基于图卷积神经网络的针对空间转录组的细胞聚类方法，包括如下步骤：

步骤1)首先进行单细胞基因表达的预处理，通过以下方式得到：

步骤1.1)计算每个基因在所有细胞中的平均表达量，低于某一阈值的基因将被排除；

步骤1.2)计算剩余的各个基因在所有细胞中的表达量的标准差，低于某一阈值的基因将被排除；

步骤1.3)针对不同批次下获得的空间转录组数据，该空间转路组数据包含：细胞空间坐标和各细胞内的基因表达量，使用Scanorama去除批次效应。批次效应(batcheffect)简单说来就是因为实验做了几个批次导致的实验误差。误差可能来源于实验环境、实验人员、试验时间的改变。因此需要对采集到的数据进行处理去除批次效应。这里采用的Scanorama源自：Hie B,Bryson B,Berger B.Efficient integration of heterogeneoussingle-cell transcriptomes using Scanorama[J].Nature biotechnology,2019,37(6):685-691。

步骤1.4)对每个细胞剩余的基因进行归一化，以细胞i为例，其第j个基因的归一化方法为：

count_ij为前一个步骤处理后得到的各基因表达量，x_ij即为所得到的细胞i的第j维特征,可得特征矩阵X＝[x_ij]^N×N。这里，1<＝i<＝N,1<＝j<＝N,i不等于j。N作为细胞总个数，由该数据集内所决定。该数据集内包含多少个细胞，N就是细胞数量。

步骤2)提取根据同一批次下细胞们的坐标信息，计算出每对细胞之间的欧氏距离，根据距离构建邻接矩阵，具体包括：

步骤2.1)对于同一个批次下的细胞组，第i个细胞二维空间坐标为v_i＝(x_i,y_i)，第j个细胞的坐标为v_j＝(x_j,_yj)，则这两个细胞间的欧式距离为d_ij＝‖v_i-v_j‖，该细胞组的距离矩阵为

L为细胞的个数；

步骤2.2)基于上述所得的距离矩阵

通过设置合适的阈值d_thres可以得到其邻接矩阵A₀＝[a_ij]^L×L，其中：

步骤2.3)对邻接矩阵进行进一步调整，引入权重因子λ实现对空间信息权重的调节，调整后的邻接矩阵为A＝λ*I+(1-λ)A₀，I为一个对角矩阵：I∈R^N×N。本实例中λ＝0.8。

步骤3)经过前两步的处理，空间转录组数据已被转化为了图结构的表示，由特征矩阵X表示细胞本身基因表达，由邻接矩阵A表示其空间邻接信息。将他们输入图卷积神经网络模型DGI，从而将细胞的基因表达同空间分布进行信息融合，得到带有空间信息的节点嵌入，具体方法为：

实验中，DGI模型包含4个图卷积层，在第l个图卷积层内，节点特征按以下公式更新：

其中H^l和H^l+1分别是第l层的输入和输出，W^(l)是该层用于特征变换的权重矩阵。

是加入了self-loops的邻接矩阵,即

I∈R^N*N,

最后一个卷积层的输出将被输入PReLU激活函数:

其中a是一个可学习的参数。

DGI的目标是学习一个编码器E来学习节点的表示：H＝E(X,A)。这里还需要引入一个读取函数S来从局部特征H提取全局特征s:s＝S(H)。对于这些局部和全局特征，DGI中训练了一个判别器D来评估每一对(h_i,s)之间的匹配程度。D(h_i,s)越高，说明更多的全局信息被包含在了当前的局部特征里。为了训练这个判别器，需要引入一个腐蚀函数(corruptionfunction)C来生成负样本。在腐蚀函数C中，原始的图结构将被随机打乱，从而获得负样本的邻接矩阵

进而基于

和X，可以获得负样本上的局部特征

最终的损失函数为:

通过最大化h_i和s之间的互信息，DGI的输出即为学习了空间信息之后的节点表示，同时节点维度也将得到降低。

步骤4)针对DGI输出的结点表示，使用PCA以进一步降低节点特征的维度。

步骤5)使用K-means++算法对PCA降维后的节点进行聚类，从而发现不同的细胞类型。在用K-means++聚类时，通过采用多次随机初始化聚类中心的办法，确保了聚类结果的鲁棒性。

步骤6)邻居丰度分析，首先构造一个邻居分布矩阵C＝[c_ij]^K×K，其中c_ij为第i类细胞周围的邻居中出现第j类细胞的次数。邻居的定义为：若存在一条边将两个结点直接相连，则两节点对应的细胞互为邻居细胞。之间然后对邻居分布矩阵进行归一化，

其中n_i为第i类细胞的个数。

步骤7)差异表达基因分析，本方法使用曼-惠特尼u检验来找到各个聚类中显著高表达的top200个基因，对于显著性指标p值小于0.5的基因进行去除。

步骤8)GO term分析，本方法使用的是各个聚类中筛选出的200个高表达基因，将其上传至GO term分析官方网站上，获得由高到低的一系列相关生物功能。其中背景基因选择为本数据集中的全部基因。

本实施例在MERFISH数据集上进行了实验，该数据集包含三个批次的细胞，分别为645,400和323个，每个细胞提供了10050种基因的表达量。通过本发明实施例中的方法对其进行处理和分析，共得到了4个具有显著生物功能差别的聚类，可以很好地对应到细胞周期的4个不同阶段。

首先，计算邻接矩阵和节点特征，得到图结构的数据。将数据输入到DGI模型中进行训练，epoch设置为5000，输出维度为256。之后将模型的输出作为结合了空间信息之后的细胞特征，使用PCA将细胞特征进一步降低到30维，并在降维后的空间对所有细胞进行聚类，类别数K设置为5。

在聚类结果上，由于其中一类只包含两个细胞，本方法只对其余4类做了进一步的分析。在空间分布上，4类细胞都主要和自己同类型的细胞相互邻近。如表1所示。

表1.邻居丰度矩阵

聚类标签	C0	C1	C3	C4
					C0	0.705208	0.086957	0.151457	0.055901
C1	0.072481	0.810832	0.093190	0.023497
					C3	0.127258	0.093938	0.768366	0.010438
C4	0.131019	0.066069	0.029115	0.772676

在差异基因表达分析和GO term分析上，4类细胞展示出了不同的生物功能，和细胞周期中的四个阶段：G1，S，G2，M可以很好的对应(如图2所示)。相比之下，现有的方法都无法直接通过聚类直接区分出不同阶段的细胞，同时各分类结果之间GO term存在严重的重复，意味着先前方法得到的聚类结果并未在生物意义上很好的将细胞区分开来。如表2所示。

表2.各聚类中top30 GO term中的平均重复率

在本申请所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另外，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口、装置或单元的间接耦合或通信连接，也可以是电的，机械的或其它的形式连接。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (9)

1.一种空间转录组细胞聚类方法，其特征在于，包括以下步骤：

对于空间转路组中每个细胞基因表达进行预处理；

根据所述空间转路组的细胞坐标生成邻接矩阵A，获得空间转录组细胞的图结构表示，由细胞特征矩阵X表示细胞基因表达，

将邻接矩阵A和细胞特征矩阵X输入经过训练的图卷积神经网络模型DGI；

所述图卷积神经网络模型DGI输出具有空间信息的结点特征表示；

对所述结点特征表示采用降维、聚类算法处理后，识别、获得所述空间转录组细胞类型。

2.根据权利要求1所述的空间转录组细胞聚类方法，其特征在于，所述的对所述结点特征表示采用降维处理，是指采用PCA算法降低所述结点特征的维度。

3.根据权利要求2所述的空间转录组细胞聚类方法，其特征在于，使用K-means++算法对PCA降维后的结点进行聚类，从而发现不同的细胞类型。

4.根据权利要求1所述的空间转录组细胞聚类方法，其特征在于，所述的对于空间转路组中每个细胞基因表达的预处理具体包括以下步骤：

步骤1.1，计算每个基因在所有细胞中的平均表达量，低于第一设定阈值的基因将被排除，高于等于第一设定阈值的基因被纳入第一基因集；

步骤1.2，计算第一基因集中各个基因在所有细胞中的表达量的标准差，低于第二设定阈值的基因将被排除，高于等于第二设定阈值的基因被纳入第二基因集；

步骤1.3，针对不同批次下获得的空间转录组数据，使用Scanorama去除批次效应，该空间转路组数据包含细胞空间坐标和各细胞内的基因表达量；

步骤1.4，根据第二基因集内的基因，对每个细胞的基因进行归一化，即，细胞i，其第j个基因归一化为：

count_ij为各基因表达量，x_ij即为所得到的细胞i的第j维特征,可得细胞特征矩阵X＝[x_ij]^N×N，1<＝i<＝N,1<＝j<＝N,i不等于j，N为细胞总数。

5.根据权利要求1所述的空间转录组细胞聚类方法，其特征在于，所述的邻接矩阵A生成方法包括：

步骤2.1，对于同一个批次下的细胞组，第i个细胞二维空间坐标为v_i＝(x_i,y_i)，第j个细胞的坐标为v_j＝(x_j,y_j)，则这两个细胞间的欧式距离为d_ij＝‖v_i-v_j‖，该细胞组的距离矩阵为

N为细胞的个数；

步骤2.2，根据距离矩阵

设置阈值d_thres，得到其邻接矩阵A₀＝[a_ij]^N×N，

其中，

步骤2.3，引入权重因子λ，获得邻接矩阵为A＝λ*I+(1-λ)A₀,I为对角矩阵，I∈R^N×N。

6.根据权利要求1所述的空间转录组细胞聚类方法，其特征在于，所述图卷积神经网络模型DGI模型包含4个图卷积层，

在第l个图卷积层内，节点特征为：

其中H^l和H^l+1分别是第l层的输入和输出，W^(l)是该层用于特征变换的权重矩阵，

是加入了self-loops的邻接矩阵,即

I∈R^N*N,

最后一个卷积层的输入为PReLU激活函数，

其中a是一个可学习的参数。

7.一种空间转录组细胞分析方法，其特征在于，在采用如权利要求1所述的空间转录组细胞聚类方法后，进一步采用邻居丰度分析、差异表达基因分析、GeneOtology(GO)terms分析来检验不同细胞聚类在空间分布特征和生物功能上的差异。

8.一种空间转录组细胞聚类装置，其特征在于，所述装置包括存储器；以及耦合到所述存储器的处理器，该处理器被配置为执行存储在所述存储器中的指令，所述处理器执行以下操作：

对于空间转路组中每个细胞基因表达进行预处理；

根据所述空间转路组的细胞坐标生成邻接矩阵A，获得空间转录组细胞的图结构表示，由细胞特征矩阵X表示细胞基因表达，

将邻接矩阵A和细胞特征矩阵X输入经过训练的图卷积神经网络模型DGI；

所述图卷积神经网络模型DGI输出具有空间信息的结点特征表示；

对所述结点特征表示采用降维、聚类算法处理后，识别、获得所述空间转录组细胞类型。

9.一种存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时，实现如权利要求1至7中任一所述的方法。

Images