CN116189785B 基于空间转录组学数据特征提取的空间域识别方法 （西安电子科技大学）

(19)国家知识产权局(12)发明专利(10)授权公告号CN116189785B(65)同一申请的已公布的文献号(73)专利权人西安电子科技大学地址710071陕西省西安市太白南路2号(72)发明人贾松卫崔议文兰猛(74)专利代理机构陕西电子工业专利中心专利代理师王品华(54)发明名称基于空间转录组学数据特征提取的空间域识别方法本发明公开了一种基于空间转录组学数据特征提取的空间域识别方法，主要解决现有技术对空间转录组数据特征提取存在过拟合及空间域识别精度低的问题。其实现方案为：对空间转录组中测量得到的基因表达数据和空间信息进行预处理；基于基因表达特征矩阵和空间信息构建基因相似性网络和空间邻域网络；对基因相似性网络和空间邻域网络进行数据增强；构建特征提取模型，并将增强的数据输入到模型中计算对比损失和重构损失；依据计算损失训练模型，将未增强的数据输入到训练好的模型中得到低维嵌入；对低维嵌入进行聚类完成对空间域识别。本发明避免了在特征提取过程中的过拟合，提高开始开始行训练结束网络21.一种基于空间转录组学数据特征提取的空间域识别方法，其特征在于，包括如下步(1)使用空间转录组测序技术对所需组织切片中每一个像素点同时测量其基因表达值和空间位置坐标，得到包含像素点-基因表达矩阵和组织切片中每个像素点空间位置的空间转录组数据；(2)对空间转录组数据的基因表达矩阵进行预处理：(2a)删除空间转录组数据中基因表达值少于三个像素点的表达基因；(2b)对删除后的数据进行数值归一化，使得各个细胞的计数和为所有细胞的中位数，再对归一化后的数据进行对数转换，并将其标准化为零均值和单位方差；(2c)对标准化后的数据进行主成分分析PCA,提取前n个主成分，生成基因表达的特征矩阵X;(3)构建空间邻域网络：(3a)基于空间坐标信息计算组织切片中每个像素点之间在空间位置上的欧式距离d:(3b)基于空间坐标计算的欧式距离d选取每个像素点的前k个最近邻，构建表征空间信息的邻接矩阵A;(3c)将步骤(2)生成的基因表达特征矩阵X作为节点属性特征矩阵；(3d)基于表征空间信息的邻接矩阵A和节点属性特征矩阵X,构成空间邻域网络G₁(A,(4)构建基因表达相似性网络：(4a)基于步骤(2)生成的基因表达特征矩阵X,计算组织切片中每个像素点基因表达值之间的欧式距离d′;(4b)基于基因表达值计算的欧式距d′,选取每个像素点的前k个最近邻，构建表征基因表达相似性的邻接矩阵B;(4c)基于表征基因表达相似性的邻接矩阵B和节点属性特征矩阵X,构成基因表达相似性网络G₂(B,X);(5a)对空间邻居网络中的边和节点属性特征依照符合伯努利分布的给定边遮掩概率p和节点特征遮掩概率p进行遮掩，得到增强后的空间邻居网络G₁(A₁,X₁);(5b)对基因表达相似性网络中的边和节点属性特征依照符合伯努利分布的给定边遮掩概率p′和节点特征遮掩概率p′进行遮掩，得到增强后基因表达相似性网络G₂(B₁,X₂);(6)构建由编码器f(·)分别与解码器h(·)和投影器g(·)级联组成的空间转录组数据的特征提取模型，并使用对比损失Ln和重构损失Lrecon的加权和作为损失函数L;(7)对空间转录组数据的特征提取模型进行训练：(7a)将数据增强后的空间邻居网络G₁(A₁,X₁)的邻接矩阵A1和节点属性特征矩阵X₁以及基因表达相似性网络G₂(B₁,X₂)的邻接矩阵B₁和节点属性特征矩阵X₂输入到空间转录组特征提取模型中，由编码器生成低维嵌入Z₁和Z₂,解码器生成重构基因表达特征矩阵X₁和2;(7b)计算低维嵌入Z₁和Z₂的对比损失和重构基因表达特征矩阵X₁和₂与节点属性特征矩阵X的重构损失，依据计算得到的损失对编码器和解码器的网络参数进行更新，直到损失3(8)将未进行数据增强的空间邻域网络的邻接矩阵A和节点属性特征矩阵X输入到(7b)和j分别代表空间邻域网络中的两个节点，若节点i包括在节点j的前k个基于空间坐标计算4(5a1)根据空间邻域网络的邻接矩阵A中的每个元素a;,依照伯努利分布采样一个边遮掩矩阵R∈{0,1}N×N,其表示如下：(5a2)将空间邻域网络的邻接矩阵A与(5a1)中生成的采样矩阵R按照元素相乘得到增(5a3)依照伯努利分布B(1-p)采样一个随机向量，生成一个维度与基因特征向量相同5的节点特征遮掩向量m,其中，p是空间邻域网络中每个节点特征向量中的值被删除的概(5a4)将空间邻居网络的节点属性特征矩阵X与(5a3)生成的节点特征遮掩向量m按照X₁=[x₁°m;x₂°m;…;xi°m;…个边遮掩矩阵R∈{0,1}×,表示如下：络中每条边被删除的概率，i和j分别代表基因表达相似性网络中的两个节点，i,j=(5b2)将基因表达相似性网络的邻接矩阵B与(5b1)中生成的采样矩阵R按照元素相乘(5b3)依照伯努利分布B(1-p)采样一个随机向量，生成一个维度与基因特征向量相同的节点特征遮掩向量m,其中p′是基因表达相似性网络中每个节点特征向量中的值被删除(5b4)将基因表达相似性网络的节点属性特征矩阵X与(5a3)生成的节点特征遮掩向量m按照元素相乘得到增强后基因表达相似性网络的节点属性特征矩阵X₂:X₂=[x₁°m;x₂°m;…;x;°m;…6其中，[;]代表着拼接操作，x;为X₂中的第i行，代表基因表达相似性网络中节点i上其基因特征向量。9.根据权利要求1所述的方法，其特征在于，步骤(6)构建的空间转录组数据特征提取所述编码器f(·),由输入图神经网络GCN层和两层隐藏GCN层级联组成，其输入维度为转录组基因特征数、第一隐藏GCN层为256维、第二隐藏GCN层为128维，在每层GCN之间使用PRelu函数作为激活函数；所述解码器h(·),由输入全连接层和三层隐藏全连接层级联组成，其输入全连接层为128维、第一隐藏全连接层为128维、第二隐藏全连接层为256维、第三隐藏全连接层维度为转录组基因特征数，在每层全连接层之间使用Relu函数作为激活函数；所述投影器g(·),由输入全连接层和隐藏全连接层级联组成，其输入全连接层为128所述损失函数，为对比损失和重构损失的加权和，其表示为：损失，L代表对比损失。10.根据权利要求1所述的方法，其特征在于，步骤(7a)中由编码器生成低维嵌入，解码器重构基因表达特征矩阵，实现如下：(7a1)编码器f(·)提取G₁(A₁,X)的节点低维嵌入Z₁和G₂(B₁,X₂)的节点低维嵌入Z₂:Z₁=f(X₁,A₁)=GCk+1(GC(X₁,Z₂=f(X₂,B₁)=GCk+1(GC(X₂,其中，GC(·)表示编码器的第k层网络，X₁和A₁分别为空间邻域网络的节点特征矩阵和邻接矩阵，X₂和B₁分别为基因表达相似性网络中的节点特征矩阵和邻接矩阵，k=1;(7a2)将(7a1)得到的两个低维嵌入Z₁和Z₂分别作为解码器h(·)的输入，得到G₁(A₁,X)重构的基因表达特征矩阵X₁和G₂(B₁,X₂)重构的基因表达特征矩阵8₂:11.根据权利要求1所述的方法，其特征在于，步骤(7b)中计算对比损失和重构损失，实现如下：(7b1)将步骤(7a)中生成的低维嵌入分别输入到投影器g(·),分别得到Z₁和Z₂对比损(7b2)根据(7b1)的结果计算每个节点i与其他节点k的对比损失1(z′,z′2i):7在G₂(B₁,X₂)作为输入时投影器的输出，i,k=1....n,n为组织切片中所有像素点的个数；(7b3)根据(7b2)得到的每个节点的对比损失，计算整个模型的对比损失Lon:(7b4)根据步骤(7a)中生成的重构基因特征矩阵X₁和₂,计算重构损失Lrecon:像素点i由G₁(A₁,X₁)重构的基因特征；xzi是8₂第i行的向量，其代表像素点i由G₂(B₁,X₂)重构的基因特征；(7b5)根据对比损失Lon和重构损失Lrecn计算损失函数L:12.根据权利要求1所述的方法，其特征在于，步骤(9)使用Leiden聚类算法对训练后得到的联合低维嵌入Z进行聚类，实现如下：(9a)根据步骤(8)中提取到的联合低维嵌入Z计算每个像素点的邻居，构建邻域图，保(9b)通过UMAP算法将联合低维嵌入Z进行降维得到维度降低后的嵌入Z′;(9c)根据(9a)中的邻域标签1’和(9b)中的嵌入Z′(9d)将聚类标签1和低维嵌入Z′进行UMAP可视化，在组织切片上对每个像素点依照聚类标签1进行染色，将相同颜色的像素点视为一个域，即实现空间域的识别。8技术领域[0001]本发明属于数据挖掘技术领域，具体涉及一种空间域识别方法，可用于为探究生物发育及疾病治疗提供参考数据。背景技术[0002]在组织切片中，一些区域具有相似性空间基因表达谱，形成组织中特异的结构或者亚结构。这些区域由于细胞类型的组成和基因表达的不同具有不同的功能分区，进而形成具有特定生物学意义结构的空间域。空间域的识别对于研究组织结构和细胞-细胞相互作用的影响至关重要。[0003]单细胞转录组测序技术scRNA-seq可以用来提供高分辨率的基因表达谱，然而由于制备样品时无法保留空间位置信息，故对于下游分析带来了局限性。空间转录组测序技术，包括基于原位杂交的成像技术和基于空间条形码的原位测序技术，其同时提供了基因表达谱和空间位置信息，该空间信息对于理解健康组织发育和疾病肿瘤微环境的生物学意义至关重要。因此空间转录组数据的提出有助于更好的描述细胞的空间组织方式。通过聚类的方式针对空间转录组数据挖掘具有相似表达模式的区域，以解读细胞的空间组织方式，即识别空间域是空间转录组学的最重要的工作之一。[0004]传统的聚类算法，如：Louvain和K-means不能有效利用可用的空间信息，的结果在组织切片中无法连续识别有明显分层结构的组织区域，无法为下游分析提供准确的参考，因此需要有针对空间转录组的数据开发同时利用基因表达谱和空间位置坐标的空间聚类方法。[0005]2021年JianHu等人在NatureMethods上提出了一个名为SpaGCN的深度学习算法，通过图卷积网络集成基因表达、空间位置和组织学图像。该算法首先结合空间位置和组织学图像构建了一个表示各点之间关系的图，再利用一个图卷积层来聚集来自邻近点的基因表达信息；然后采用无监督迭代聚类算法，利用聚合表达式矩阵对点进行聚类。[0006]2021年EdwardZhao等人在NatureBiotechnology上提出了一个名称为BayeSpace的算法。BayesSpace算法通过对基因表达矩阵的低维表示进行建模，通过贝叶斯统计方法在先验算法中引入空间邻居结构鼓励相邻像素点属于同一簇，从而实现空间聚[0007]2022年ShihuaZhang等人在NatureCommunications上提出了一种新的基于图注意力自编码器的框架STAGATE,其利用图注意力自编码器，在嵌入空间信息的同时通过注意力机制自动学习节点间边的权重，考虑了空间域边界像素点的空间相似性。[0008]2022年ChangXu等人在NucleicAcidsResearch上提出了一种深度神经网络框架DeepST,其用神经网络提取组织学图像特征，并与基因表达和空间位置创建空间增强的基因表达矩阵，使用图卷积网络和去噪自编码器联合生成增强ST数据的潜在表示。[0009]以上这些算法均存在以下不足：[0010]一是由于加入了组织学特征，因而在提高聚类精度的同时，增加了模型的复杂度，9[0011]二是同时一些算法因为过分重视空间信息从而过矫正了基因表达特征，造成了聚类的过拟合，故无法识别出一些精细的区域，不能对识别出的空间域进行生物学功能的精准分析。[0012]三是不稳定，每次复现的结果差异较大，并且只在基于原位测序的空间转录组测序手段测量的数据集上表现好，而在基于成像的数据集上表现很差，不能对空间转录组的数据集进行广泛的分析。发明内容[0013]本发明的目的在于针对上述现有技术的不足，提出一种基于空间转录组学数据特征提取的空间域识别方法，以提取空间转录组数据中基因表达谱和空间位置信息的联合特征，提高在空间转录组学基于测序和基于成像这两种测序手段上的泛化能力，完成对空间域生物学功能的精准分析。[0014]本发明的技术方案是：对空间转录组中测量得到的基因表达数据和空间信息进行预处理；基于基因表达特征矩阵和空间信息构建基因相似性网络和空间邻域网络；对基因相似性网络和空间邻域网络进行数据增强；构建特征提取模型，并将增强的数据输入到模型中计算对比损失和重构损失；依据计算损失训练模型，将未增强的数据输入到训练好的模型中得到低维嵌入；对低维嵌入进行聚类完成对空间域识别。其实现步骤包括如下：[0015](1)使用空间转录组测序技术对所需组织切片中每一个像素点同时测量其基因表达值和空间位置坐标，得到包含像素点-基因表达矩阵和组织切片中每个像素点空间位置的空间转录组数据；[0016](2)对空间转录组数据的基因表达矩阵进行预处理：[0017](2a)删除空间转录组数据中基因表达值少于三个像素点的表达基因；[0018](2b)对删除后的数据进行数值归一化，使得各个细胞的计数和为所有细胞的中位数，再对归一化后的数据进行对数转换，并将其标准化为零均值和单位方差；[0019](2c)对标准化后的数据进行主成分分析PCA,提取前n个主成分，生成基因表达的特征矩阵X;[0021](3a)基于空间坐标信息计算组织切片中每个像素点之间在空间位置上的欧式距离d;[0022](3b)基于空间坐标计算的欧式距离d选取每个像素点的前k个最近邻，构建表征空间信息的邻接矩阵A;[0023](3c)将步骤(2)生成的基因表达特征矩阵X作为节点属性特征矩阵；[0024](3d)基于表征空间信息的邻接矩阵A和节点属性特征矩阵X,构成空间邻域网络G₁[0025](4)构建基因表达相似性网络：[0026](4a)基于步骤(2)生成的基因表达特征矩阵X,计算组织切片中每个像素点基因表达值之间的欧式距离d′;[0027](4b)基于基因表达值之间的欧式距d′,选取每个像素点的前k个最近邻，构建表征基因表达相似性的邻接矩阵B;[0028](4c)基于表征基因表达相似性的邻接矩阵B和节点属性特征矩阵X,构成基因表达[0030](5a)对空间邻居网络中的边和节点属性特征依照符合伯努利分布的给定边遮掩概率p和节点特征遮掩概率p进行遮掩，得到增强后的空间邻居网络G₁(A₁,X);[0031](5b)对基因表达相似性网络中的边和节点属性特征依照符合伯努利分布的给定边遮掩概率p和节点特征遮掩概率p进行遮掩，得到增强后基因表达相似性网络G₂(B₁,X₂);[0032](6)构建由编码器f(·)、并行的解码器h(·)和投影器g(·)级联组成的空间转录组数据的特征提取模型，并使用对比损失Lon和重构损失Lrec的加权和作为损失函数L;[0033](7)对空间转录组数据的特征提取模型进行训练：[0034](7a)将增强后的空间邻居网络G₁(A₁,X₁)的邻接矩阵A₁和节点属性特征矩阵X₁以及基因表达相似性网络G₂(B₁,X₂)的邻接矩阵B₁和节点属性特征矩阵X₂输入到空间转录组特征提取模型中，由编码器生成低维嵌入Z₁和Z₂,解码器生成重构基因表达特征矩阵X₁和X₂;[0035](7b)计算低维嵌入Z₁和Z₂的对比损失和重构基因表达特征矩阵X₁和X₂与节点属性特征矩阵X的重构损失，根据计算的损失对编码器和解码器的网络参数进行更新，直到损失函数L收敛，得到训练好的空间转录组特征提取模型；[0036](8)将未进行数据增强的空间邻域网络的邻接矩阵A和节点属性特征矩阵X输入到(7b)中训练好的空间转录组特征提取模型中，得到包含空间信息和基因表达的联合低维嵌入Z;[0037](9)使用Leiden聚类算法对得到的联合低维嵌入Z进行聚类，在组织切片上得到具有一致性基因表达的区域，即空间域。[0038]本发明与现有技术相比，具有如下优点：[0039]1)本发明由于结合空间转录组数据的空间信息和基因表达谱构建空间邻域网络和基因表达相似性网络，与现有方法相比，可以更好的平衡空间信息和基因表达信息，防止对基因表达谱的过拟合，提高了空间域识别的准确性和鲁棒性。[0040]2)本发明由于只使用了基因表达谱和空间信息，没有加入组织学图像的特征，因此提升了模型的效率，既减少了运行的时间，又可以应对日后产生的更大的数据集的挑战。[0041]3)本发明由于引入对比损失训练低维嵌入，使样本中相似的样本更相似，不相似的样本更远离，很好的拟合了空间聚类的问题，与现有技术相比，提高了在空间转录组两种测序手段生成的数据集上的泛化能力。[0042]4)本发明由于设计了编码器、解码器级联的模型架构，同时考虑对比损失和重构损失，因而会生成去噪后的数据，更好的保留了原始样本中的生物意义。附图说明[0043]图1为本发明的实现流程图；[0044]图2为本发明中进行数据增强的示意图；[0045]图3为本发明中构建的空间转录组数据的特征提取模型图；[0048]现有的空间转录组数据包括基于原位杂交的成像技术和基于空间条形码的原位测序技术，其中成像技术包括STAPmap、MERFISH,原位测序技术包括spatial类背外侧前额皮质层151673切片的空间转录组数据集为例，该数据集中包含3639个像素[0051]1.1)获取10xVisium的空间转录组测序的人类背外侧前额皮质层151673切片的空间转录组数据集中的像素点-基因表达矩阵数据，删除空间转录组数据中基因表达值中[0052]1.2)将过滤后的转录组数据进行中值归一化，即每列数据除以该列数据的中位[0053]1.3)对标准化后的数据进行主成分分析PCA,提取前300个主成分，生成基因表达[0057]2.1)获取10xVisium的空间转录组测序的人类背外侧前额皮质层151673切片的间信息的邻接矩阵A:[0062]其中，为空间邻域网络的邻接矩阵A中第i行第j列的元[0064]2.4)基于表征空间信息的邻[0066]3.1)基于步骤1中生成的基因表达特征矩阵X,计算组织切片中每个像素点基因表[0068]其中x;k和xik分别是像素点i和像素点j基因表达特征向量的第k维度的值，k=[0072]3.3)基于表征基因表达相似性的邻接矩阵B和节点属性特征矩阵X,构成基因表达节点属性特征矩阵X及基因表达相似性网络的邻接矩阵B和节点属性特征矩阵X进行数据增[0076]4.1)根据空间邻域网络的邻接矩阵A中的每个元素a;;,依照伯努利分布采样一个边遮掩矩阵R∈{0,1}N×N:[0089]4.6)将基因表达相似性网络的邻接矩阵B与4.5)中生成的采样矩阵R按照元素相乘，b,为基因表达相似性网络的邻接矩阵B中第i行第j列的元素，R;为采样矩阵R中第i行[0092]4.7)依照伯努利分布B(1相同的节点特征遮掩向量m,其中p'=0.3是基因表达相似性网络中每个节点特征向量中[0093]4.8)将基因表达相似性网络的节点属性特征矩阵X与4量m按照元素相乘得到增强后基因表达相似性网络的节点属代表重构损失，L代表对比损失。[0105]步骤6,对空间转录组数据的特征提取模型进行训练。[0106]6.1)通过编码器f(·)分别提取G₁(A₁,X₁)的节点低维嵌入Z₁和G₂(B₁,X₂)的节点低维嵌入Z₂:阵和邻接矩阵，X₂和B₁分别为基因表达相似性网络中的节点特征矩阵和邻接矩阵，k=1;[0110]6.2)将步骤6.1)得到的两个低维嵌入Z和Z₂分别作为解码器h(·)的输入，得到G₁(A₁,X,)重构的基因表达特征矩阵X₁和G₂(B₁,X₂)重构的基因表达特征矩阵X₂:[0113]6.3)将步骤6.1)中生成的两个低维嵌入Z₁和Z₂分别输入到投影器g(·),分别得到Z₁的对比损失低维嵌入Z′和Z₂的对比损失低维嵌入Z′2:[0116]6.4)根据步骤6.3)的结果计算每个节点i与其他节点k的对比损失1(z′1,z′2;):素点i在G₂(B₁,X₂)作为输入时投影器的输出，i,k=1n,n=3639为组织切片中所有像素点的个数；[0119]6.5)根据步骤6.4)得到的每个节点的对比损失，计算整个模型的对比损失Lcon:[0121]6.6)根据步骤6.2)中生成的重构基因特征矩阵X₁和X₂,计算重构损失Lrecon:代表像素点i由G₁(A₁,X,)重构的基因特征；Xzi是X₂第i行的向量，其代表像素点由G₂(B₁,X₂)重构的基因特征；[0124]6.7)根据对比损失L和重构损失Lrecon计算损失函数L:[0126]6.8)根据步骤6.7)得到的损失函数L对编码器和解码器的网络参数进行更新，直到损失函数L收敛，得到训练好的空间转录组特征提取模型。[0127]步骤7将未进行数据增强的空间邻域网络的邻接矩阵A和节点属性特征矩阵X输入到步骤6训练好的空间转录组特征提取模型中，得到包含空间信息和基因表达的联合低维嵌入Z;[0128]步骤8使用Leiden聚类算法对步骤7得到的联合低维嵌入进行聚类。[0129]8.1)根据步骤7中提取到的联合低维嵌入Z计算每个像素点的邻居，构建邻域图，[0131]8.3)根据步骤8.1)得到的邻域标签1’和步骤8.2)得到的低维后嵌入Z′通过Leiden算法得到聚类标签1;[0132]8.4)将聚类标签1和低维嵌入Z′进行UMAP可视化，在组织切片上对每个像素点依照聚类标签1进行染色，将相同颜色的像素点视为一个域，即实现空间域的识别。[0133]以下结合模拟实验，对本发明的技术效果进行描述。[0135]仿真实验的计算机硬件CPU为IntelCore[0138]仿真1:用本发明和现有的6种方法SEDR、STAGATE、DeepST、scanpy、stlearn、SpaGCN在基于原位测序的10xVisium人类背外侧前额皮质层DLPFC的12个切片的空间转录组数据集和基于成像的STARmap小鼠视觉皮层的空间转录组数据集这两种空间转录组测序手段生成的数据集上进行空间聚类，并使用调整兰德指数ARI作为评价各方法的空间聚类结果的评价指标，结果如表1:[0139][0142]SEDR,LingS,HRepresentationofSpatialTranscriptomics[J].bioRxiv,2021.resolvedtranscriptomicswithanadap.Naturecommunications,2022,13(1):1-12.spatialtranscriptomicsbydeeplearning[J].NucleicAcidsResearc[0145]Scanpy,WolfFA,AngererP,TheisFJ.SCANPY:large-sgeneexpressiondataanalysis[J].Genometissuemorphologyandgeneexpressio[0147]SpaGCN,LiM,HuJ,LiX,etal.SpaGCN