摘要

空间分辨转录组学技术能够全面测量完整组织背景下的基因表达模式。然而,现有技术存在分辨率低或测序深度浅的问题。在此,我们提出了一种基于深度学习的方法 DIST,它能对未测量位置的基因表达谱进行推算,并通过自我监督学习和迁移学习来增强原始测量点和推算点的基因表达。我们评估了 DIST 在估算、聚类、差异表达分析和功能富集分析方面的性能。结果表明,DIST 可以准确地估算基因表达量,提高低质量数据的基因表达量,帮助检测更多有生物学意义的差异表达基因和通路,从而更深入地了解生物学过程。

关键词:空间转录组学;归因;去噪;自监督学习;迁移学习

引言

空间转录组学(ST)是近年来的一项技术创新,它可以测量组织中具有空间信息的基因表达。了解转录本的相对位置对于理解生物、生理和病理过程至关重要,因此,空间转录组学技术已被应用于多个生物学领域,如肿瘤微环境[1, 2]、胚胎发育[3, 4]和神经学[5, 6]。ST 技术主要可分为两类:(1) 基于下一代测序(NGS)的方法,如 ST [7] 和 10X Genomics Visium [8],在测序前将位置信息编码到转录本上;(2) 基于成像的方法,包括基于原位测序的方法、FISSEQ[9]和 STARmap[10]为代表的原位测序方法,其中转录本在组织中进行扩增和测序;以及 MerFISH [11] 和 seqFISH [12]为代表的原位杂交方法,其中成像探针在组织中依次杂交[13]。

这些技术通常需要在分辨率和基因通量之间做出权衡。基于成像的方法可提高分辨率和灵敏度,但通量较低,限制了其探索全转录组相互作用和发现新序列的潜力。基于 NGS 的方法无偏见,因为它们原则上可以对整个转录组进行采样,但分辨率和灵敏度较低 [13-15]。最近,基于 NGS 的方法的分辨率迅速提高,例如,DBiT-seq [16] 的分辨率达到了 10 μm,Stereo-seq [17] 的分辨率达到了纳米级(220 nm)。虽然这些基于 NGS 的方法能测量捕获位置(称为点)上千个基因的表达水平,但它们的空间覆盖范围不完整,限制了它们在研究详细表达模式方面的作用。ST 的斑点直径为 100 μm,中心到中心的距离为 200 μm;Visium 的斑点直径为 55 μm,中心到中心的距离为 100 μm。由于扩散距离为 1 μm,DBiT-seq 的理论中心到中心距离可小至 2 μm。Stereo-seq 的标准 DNA 纳米球芯片的斑点直径为 ∼220 nm,中心到中心的距离为 500 或 715 nm。

测量点的直径和密度都限制了目前基于 NGS 方法的空间分辨率,从多细胞到亚细胞都有 [14]。为了获得高分辨率(HR)空间表达谱,人们开发了新的方法。例如,BayesSpace [18] 采用贝叶斯方法,利用空间邻域的信息提高分辨率,但 BayesSpace 通过推断子点的表达来提高分辨率,无法预测未测量位置的表达;XFuse [15] 使用由空间表达数据和相应组织学图像训练的深度生成模型来推断超分辨率表达图;HisToGene [19] 是一种深度学习方法,采用 Vision Transformer 从组织学图像中预测基因表达SpatialPCA [20] 利用改进的概率主成分(PC)分析,提取归一化和缩放 ST 数据的低维表示,由于数据的生成性,可以在未测量的点上归纳空间 PC。但是,由于 NGS 数据的生成性,它无法解释原始计数中存在的均方差关系,从而导致潜在的功率损失[21]。在此,我们主要关注基于 NGS 的非密集空间覆盖 ST 技术。考虑到深度学习在 ST 方面的成功应用[22-24]、基因表达与空间信息和图像之间的本质联系以及受基于学习的图像处理的启发[25, 26],我们提出了 DIST 方法,该方法仅利用空间基因表达数据对新的和未测量的点进行基因表达谱推算,从而得到分辨率更高的精细空间图谱。与 XFuse 和 HisToGene 不同的是,DIST 利用深度学习增强了 ST 数据,不需要组织学图像等任何附加信息。此外,尽管基于 NGS 的 ST 技术有所改进,但包括扩增偏差在内的各种技术因素,尤其是较低的 RNA 捕获率,都会导致测序过程中出现大量噪音[27]。由于测序相对较浅,相对 HR 技术面临的挑战更大。由于扩增和丢失事件造成的噪声可能会阻碍分析并破坏潜在的生物信号。为了解决这个问题,DIST 可以通过从高质量数据中学习,对估算的空间基因表达进行去噪处理。总之,我们将 DIST 解释为空间转录组学的去噪和归因。我们通过模拟 STARmap 数据集和 Stereo-seq 数据集来说明 DIST 的优势,并将其应用于分别用 ST 和 10X Visium 平台获得的两个真实世界数据集。结果表明,与传统的插值方法相比,DIST 能更准确地完成估算任务,同时揭示更多的生物信号。

方法

DIST 概述

图 1.DIST 的结构。(A) 矩阵(上)和蜂巢(下)排列平台上的归因和向下取样草图。估算:在每两个相邻点之间预测一个新点。向下采样:以一个点为间隔定期提取非相邻点。

DIST 主要侧重于基于阵列的 ST 技术。这些技术的测量点按一定的模式排列,一般包括两类:矩阵排列(如 ST)和蜂巢排列(如 Visium)。DIST 在每两个相邻的测量点之间预测一个未测量的点(图 1A)。DIST 的整体思路是通过自我监督学习来推算未测量斑点的基因表达。由于我们没有用于训练的低分辨率(LR)和高分辨率(HR)基因表达对,因此我们对原始基因表达图进行下采样,创建人工 LR 基因表达图来训练模型。如图 1B 和 C 所示,DIST 首先以监督的方式训练一个插值深度神经网络--变异学习网络(VLN)[26],使用下采样 LR 基因表达作为输入,原始基因表达作为 HR 标签。输出结果是与原始数据分辨率相同的估算基因表达图谱。然后在预测阶段(或估算阶段),将所有点的原始基因表达量输入训练有素的模型,以输出 HR 基因表达量DIST 首先学习如何从向下采样的 LR 基因表达图和原始基因表达图中推算未测量的基因表达,然后应用学习到的规则在原始数据上推算未测量的基因表达。这种结构允许 DIST 通过各种方式增强 ST 数据。除了通过自我监督学习对未测量点进行归因外,DIST 还可以通过在高质量 ST 数据上进行迁移学习来改进低质量 ST 数据的归因,并通过在训练集输入中引入合成噪声来对归因表达进行去噪

(B) 自我监督学习框架。下采样:对原始基因表达图进行下采样,创建人工低分辨率基因表达图,形成 LR-HR 基因表达对。训练过程:使用下采样 LR 基因表达作为输入,原始基因表达作为标签,以监督方式训练插值深度神经网络 VLN。预测过程:将原始基因表达图谱输入训练好的模型,以输出更高分辨率的基因表达。

DIST 方法

仅使用带有空间坐标的基因表达矩阵来完成任何数据空间的估算,如基因计数、归一化基因计数和其他嵌入。本文主要分析基因计数矩阵。

我们考虑一项 ST 研究,收集 m 个点上 n 个基因的基因表达矩阵X_{m\times n }。这些点具有已知的空间坐标,例如第 i 个点的坐标为 (x_{i }, y_{i }) ∈ Z^{2}所有点的坐标合并成一个坐标矩阵 C_{m\times 2 },其第 1 行和 X_{m\times n } 的第 1 行指向同一个点。X_{m\times n } 的第 j 列(用 Xj 表示)是一个 m 向量,代表 m 个点中第 j 个基因的表达量。X^{j} 的值可以通过相应的空间坐标进行协调,并转换为形状为(u,v)的矩阵I^{j} 。具体来说,I^{j} 的第 x 行和第 y 列上的值就是 X^{j}的第 i 个元素。组织内矩阵 M 的形状与 I^{j}相同,只包含 0 和 1,我们用它来表示是否存在斑点。例如,如果 M 的第 p 行和第 q 列上的值等于 1,则存在相应坐标(p,q)的斑点;如果等于 0,则不存在斑点。这里,I^{j} 可视为表达图,注意图集 I = {I^{1} , --- , I^{n} }。我们用函数 T (·) 表示从空间基因表达到图的转换。函数

在将基因表达转化为表达图时,由于组织切片的图形不规则以及过滤点的质量较差,大多会出现缺失值。我们不考虑没有细胞的组织外点。因此,我们的方法侧重于组织内的斑点,并恢复组织切片内部的细节。

(C) VLN 的网络架构。

DIST 通过基于学习的插值法完成空间表达估算DIST 的插值部分是一个 VLN,用于估计 HR 值与其最近的 LR 值之间的变化值[26](图 1C)。我们在由基因表达矩阵 Xtrain 和坐标矩阵 Ctrain 构建的数据集 I_{train} = {I_{train}^{1}, --- , I_{train}^{ntrain} } 上训练 VLN,并在预测或测试阶段使用由基因表达矩阵 X_{test}和坐标矩阵 C_{test} 构建的 I_{test} = {I_{test}^{1}, --- , I_{test}^{ntest} }:

其中 n_{train} 和 n_{test} 分别表示训练数据集和测试数据集的基因编号,I_{train}^{j}M_{train}的形状相同(utrain、vtrain),I_{test}^{j}M_{test} 的形状相同(u_{test}v_{test})。

在训练阶段,DIST 创建 n_{train}对 的 LR 和 HR 表达图 \left \{ {(L_{train}^{j},I_{train}^{j} )}\right \}_{j=1}^{train}来训练 VLN,其中 L_{train}^{j}是通过向下采样从I_{train}^{j}中提取的(图 1A)。这一过程可表示如下:

我们将下采样函数记为 D (·),因此函数上

L_{train}^{j} 的形状为([u_{train}/2] , [v_{train}/2 ]),其中[·]表示四舍五入操作。

F (·) 表示从输入的 LR L_{train}^{j}恢复 HR I_{train}^{j} 的网络。我们采用以下均方误差来训练 VLN:

我们使用 ADAM 优化器,通过小批量梯度下降来最小化损失

请注意,该网络可用于不同规模的输入。在预测或测试阶段,将测试集I_{test}输入训练有素的网络:

H_{test}^{j}的形状等于 (2u_{test},2v_{test})

接下来,将 HR H_{test}^{j}转换为估算的基因表达矩阵

其中,h (·) 从原始空间信息中学习组织边界,并恢复归纳基因表达的边界。h (M_{test}) (p, q) = 1,当且仅当位置(p, q)指向一个原始点或位于两个相邻原始点之间。Y_{test} = (Y_{test}^{1}, --- , Y_{test}^{ntest} ) 是估算的 X_{test}D_{test}Y_{test}对应的坐标矩阵。估算输出 Y_{test} 可能会对与 X_{test} 相同的点上的表达式进行修改,因为点数的增加会改变尺度,而微小的修改可以增强一致性。但表达量的变化不会太大,基因间的皮尔逊相关系数(PCC)大多接近 1。

如果 I_{train}I_{test}来自相同的 ST 数据,那么学习过程就是自监督学习;如果不是,那么就是迁移学习。

VLN 的结构

(C) VLN 的网络架构。
迭代
沿通道轴堆叠

如图 1C 所示,VLN 采用递归卷积结构,以监督的方式执行渐进恢复路径[26]。假设有一个 HR 表达图 y,它被分成四个部分:x_{tl}(= x)、x_{tr}x_{bl}x_{br},这四个子图由从每 2 × 2 个非重叠斑块中提取的左上、右上、左下和右下值组成。我们假设 x 的形状为 (u, v)。

首先,网络提取 LR 特征,并通过递归卷积层对其进行增强。假设有 K 个递归层,每个递归层包含 L 个卷积层。让f_{in}^{k} 表示第 k 个递归子网络的输入。子网络的输出f_{out}^{k}按如下方式逐步更新:

其中,f_{l}^{k} 是第 l 层的特征图,W_{in}^{k}W_{l-1}^{k}b_{in}^{k}b_{l-1}^{k}是卷积层的滤波器参数和偏置

在这

我们必须连接输入 x 和 K 个递归子网络

其次,通过对增强特征的最后一次卷积重建变化图,计算方法如下:

其中,x_{tl}x_{tr}x_{bl}x_{br} 分别为每个 2 × 2 非重叠片段中相对于输入 x 的左上、右上、左下和右下的估计差值,W、b 分别表示重构层的滤波器参数和偏置。

第三,网络会将每个非重叠斑块中的值变化和相应的左上角值结合起来,从而得到

形状为 (4, u, v)。

最后,通过位置感知上采样层重建 HR 表达图。输出结果用 \hat{xup} 表示,计算公式如下:

其中,"\left \lceil \right \rceil"表示上限操作,"p "表示行索引,"q "表示列索引,"r "表示指向 "x_{tl}"、"x_{tr}"、"x_{bl} 或"x_{br}"的分量索引。这里,所有参数 = {W, b} ∪ {W_{in}^{k}, W_{2}^{k} --- , W_{L-2}^{k}, b_{in}^{k}, b_{2}^{k}, --- , b_{L-2}^{k} }_{k=1}^{K},其中滤波器参数初始化为正态分布,偏置初始化为 0。

在训练阶段,默认 K = 2,L = 5。前端层的学习率设置为 0.001,重建层的学习率设置为 0.00001。在大多数情况下,损失会在 200 个历时内收敛。

引入合成噪声

DIST 可以通过在训练过程中引入合成噪声来对估算的表达式进行去噪处理。按照 SAVER [28],我们在基因表达中引入了如下噪声:对于基因 jon spot i,我们将原始计数作为 X_{ij }处理,而噪声值 X_{ij } 则从泊松分布中提取生成,Y_{ij } ∼ Poisson (λX_{ij } ),其中 λ 是效率损失。在随后的实验中,我们让 λ = 0.6 为 Visium 人类浸润性导管癌(IDC)数据[29] 引入噪声。

数据分析

我们按照质量控制、归一化、聚类、差异表达分析和通路富集等常规分析流程[30]探索 IDC 的生物学功能。我们采用对数转换归一化方法,选择了 2000 个高变异基因,并应用 Louvain [31] 方法对斑点进行聚类。由于现有的识别重要差异表达(DE)基因的方法大多适用于整数的唯一分子标识符(UMI)计数,因此我们在差异表达分析前对估算的计数进行了四舍五入。我们使用 t 检验来寻找 DE 基因,并根据调整后 P 值和 log2-foldchange 的阈值来选择显著的 DE 基因,这些值列举如下:

(i) 原始 IDC 与估算 IDC 的比较:调整后 P 值小于 0.05 且 log2-foldchange > 1;

(ii) 去噪 IDC 与非去噪 IDC 的比较:调整后 P 值小于 1×10-6 且 log2-foldchange > 0.6。富集分析在 Metascape [32] 网站上完成。

结果

DIST 准确估算了空间基因表达

为了严格评估 DIST 的估算性能,我们从 STARmap 数据集中创建了一个模拟 ST 数据,在该数据集中,903 个基因的 RNA 分子被测量并定位在小鼠胎盘组织上[33]。我们利用这一分子分辨率数据,将所有 RNA 投射到伪点(240 × 240 像素)上,从而模拟出单细胞分辨率的点状数据。240 × 240 像素范围内的总 RNA 计数被视为该模拟伪斑的基因表达量。因此,共有 7465 个伪点的 RNA 计数超过 5 个。伪点的数量与 ClusterMap [33] 估算的细胞数大致相当。我们将该模拟单像素分辨率数据集作为地面实况,并对其进行下采样,以创建模拟粗分辨率 ST 数据。然后,我们将 DIST 与其他插值方法进行了比较,包括近邻插值法(NN)、线性双曲插值法(Linear)、立方样条插值法(Cubic)、新边缘定向插值法(NEDI)[34]、随机森林快速图像插值法(FIRF)[35]和空间PCA。我们没有将 DIST 与贝叶斯空间(BayesSpace)、XFuse 和 HisToGene 等用于 ST 归因的先进方法进行比较,因为贝叶斯空间将一个点分成几个子点,无法归因未测量位置的基因表达;XFuse 和 HisToGene 需要组织学图像,而模拟无法提供。由于模拟数据中的基因数量比整个转录组的基因数量少一个数量级以上,我们用 IDC 数据对 DIST 进行了预训练,然后在模拟数据本身构建的训练集上对模型进行了微调。

(A) 基于粗分辨率模拟数据、单细胞分辨率地面实况以及使用 DIST、线性、立方、NEDI 和 NN 估算的数据,具有不同空间模式的多个基因的空间表达。一致的颜色范围按基因-地面实况的最小值和最大值缩放。

图 2A 和补充图 S1 显示了具有不同空间模式的几个基因的地面实况、模拟粗分辨率表达和使用 DIST、Linear、Cubic、NEDI 和 NN 估算的表达。由于空间覆盖和质量控制不完全,模拟表达存在空缺位置。DIST 和上述插值方法都能填补空缺表达,确保组织切片的完整性和连续性。直观地说,DIST 比其他插值方法具有更清晰的基因表达谱和更精细的描述,尽管其他方法也能达到与 DIST 相同的分辨率。NN 使插值表达图不连续,容易出现锯齿现象。线性和立方都能提高平滑度,但会导致模糊。这些插值方法仅利用周围位置的点来推算未测量点的表达,而没有利用其他点的有效知识。而 DIST 则考虑了所有点的空间表达及其内部统计信息。此外,DIST 是在整个基因集上训练的,利用了所有基因和斑点的优势,而不是像其他方法那样只利用一个基因和周围的斑点。

我们计算了 903 个基因的地面实况与所有点的每个估算表达之间的 PCC,结果如图 2B 所示。与表现最好的插值法相比,DIST 获得了更高的基因间 PCC 中位数(中位数:DIST = 0.523,Linear = 0.389)。DIST 的 PCC 第 25 百分位数(0.498)甚至高于插值法的最大第 75 百分位数(Linear = 0.477)。我们计算了 SpatialPCA 在比例嵌入空间中的 PCC,因为 SpatialPCA 只能获得估算的比例数据。需要注意的是,SpatialPCA 的目的是降维,而不是预测未测量点的真实表达值。FIRF 的性能不佳,主要是因为它是为图像插值设计的,并且是在图像数据集上训练的。

(B) 使用 DIST、Linear、Cubic、NEDI、NN、SpatialPCA 和 FIRF 计算的地面实况与估算表达之间的基因相关性。方框表示第 25、50 和 75 百分位数。胡须表示所有非离群值的范围,定义为距离铰链 1.5 个四分位数间范围内的观测值。(C) 从立体序列成年小鼠半脑创建的模拟。使用 DIST_transfer(迁移学习)、DIST_self(自我监督学习)、Linear、Cubic、NEDI、NN、XFuse、SpatialPCA 和 FIRF 计算基本真实表达和估算表达之间的基因相关性。方框中的解释与(B)一致。

接下来,我们从立体测序数据中创建了另一个模拟 ST 数据,该数据提供了成年小鼠半脑的亚细胞分辨率空间表达[17]。与上述过程类似,我们计算了地面实况与 DIST、Linear、Cubic、NEDI、NN、XFuse、SpatialPCA 和 FIRF 估算表达之间的 PCC(图 2C)。对于 DIST,我们分别评估了来自自我监督学习和基于 IDC 数据训练的迁移学习的估算值。此外,我们还利用组织切片的核酸染色图像运行了 XFuse。XFuse 可以按像素预测表达,分辨率与染色图像一样高。为了进行比较,我们将像素按其位置融合成斑点。图 2C 显示,DIST 获得了最高的基因表达中位数,而迁移学习提高了预测的准确性。与 STARmap 模拟数据相同,DIST 的 PCC 是所有方法中最高的。由于测序深度较浅,Stereo-seq 模拟数据中 DIST 的 PCC 中位数(0.465)低于 STARmap 数据。但在使用高质量数据进行迁移学习后,Stereo-seq 数据的 PCC 中值提高到了 0.518,与 STARmap 数据相当接近。

探索人类黑色素瘤 ST 数据的空间模式

接下来,我们测试了 DIST 是否有助于发现 ST 的空间模式。我们将基于扩散模型的转录本特征空间模式识别方法 Sepal 应用于使用 DIST 估算前后的人类黑色素瘤 ST 数据集[7]。按照 Sepal 的教程,首先,我们将基因表达谱按空间结构程度从明显模式到随机模式进行排序;然后,将具有明显空间模式的排名靠前的基因分成模式家族,同一模式家族的成员表现出相似的空间组织;最后,根据基因本体论(Gene Ontology:生物过程(GO:BP)数据库[36]进行功能富集分析。

图 3A 显示了一些在黑色素瘤估算数据中排名靠前但在原始黑色素瘤数据中排名靠后的转录本特征。CSPG4 与黑色素瘤肿瘤的形成和某些黑色素瘤的不良预后有关[37]。DLL3 在许多转移性黑色素瘤中都有表达,靶向该基因可能是一种很有前景的治疗策略,可防止炎症加重黑色素瘤的进展[38]。CD37 几乎只在免疫系统细胞中表达,尤其是在成熟的 B 细胞中[39]。如图 3A 所示,这些基因的表达具有独特的空间结构。较低的分辨率可能导致这些基因在原始黑色素瘤数据中的排名相对较低。而 DIST 估算的 HR 数据使 Sepal 更有说服力,大大提高了这些基因的排名。

图 3.探索 ST 人类黑色素瘤数据的空间模式。(A) 三个与疾病相关的基因在估算数据(下图)中的排名较高,因为它们具有独特的空间模式,但在原始数据(上图)中的排名较低。每个转录本档案的页眉给出了相关基因的名称和排名。(B) 原文中病理学家对苏木精和伊红染色组织图像的注释,其中黑色:黑色素瘤,红色:基质,黄色:淋巴组织。

接下来,我们将排名前 150 位的基因划分为三个主要的模式家族,并进行了富集分析。从每个模式族的代表性主题(补充图 S2)、独特的标记基因(补充图 S3 和 S4)和生物学过程来看,这些模式族与组织学注释(图 3B)有明显的联系,包括黑色素瘤、淋巴细胞和基质。然后,我们比较了原始 ST 数据和估算 ST 数据中具有相同空间模式的这些基因组。推算的 ST 数据中排名最高的 150 个基因共分为 87 个黑色素瘤相关基因、37 个淋巴相关基因和 25 个基质相关基因,而原始 ST 数据中排名最高的 150 个基因共分为 108 个黑色素瘤相关基因、31 个淋巴相关基因和 8 个基质相关基因(图 3C)。就三个家族而言,尽管黑色素瘤相关基因较少,但估算数据比原始数据富集了更多的 GO:BP 术语(图 3C)。此外,淋巴家族富集了 124 个新通路,其中多个过程与免疫反应和淋巴细胞活化有关。基质家族则有更多与生长因子有关的通路(图 3D)。我们还选择了 100 或 200 个排名靠前的基因重复上述分析步骤,结果与 150 个排名靠前的基因一致(附图 S5)。这表明 DIST 估算数据中的排名靠前的基因更具有区域特异性和生物学意义,尤其是在淋巴区域。与黑色素瘤区域相比,淋巴区域和基质区域相对较小。由于淋巴区域和基质区域较小,与黑色素瘤等较大区域相比,淋巴相关基因和基质相关基因的排名要低得多。DIST 提高了分辨率,从而使 Sepal 更有能力找到淋巴和基质相关基因

(B) 原文中病理学家对苏木精和伊红染色组织图像的注释,其中黑色:黑色素瘤,红色:基质,黄色:淋巴组织。(C) 每个家族的基因数量(左)和有意义(P 值小于 0.05)的 GO:BP 术语(右)。绿色符号数字来自原始数据;橙色符号数字来自估算数据。(D) 在淋巴(上图)和基质(下图)家族中,仅通过估算数据确定的前 10 个有意义的 GO:BP 术语。

探索 Visium 人类 IDC 数据的生物学功能

我们进一步分析了在 Visium 平台上制备的人类 IDC 样本,随后进行了聚类、差异表达分析和通路富集等分析。我们感兴趣的是,与原始 ST 数据相比,DIST 能带来哪些新发现。

(A、B)使用两种策略对估算的 Visium 人类 IDC 进行卢万聚类显示。
(A) 基于自我监督学习策略。
(B)基于在高质量数据(Visium 小鼠矢状脑后部)上训练的迁移学习模型。

我们尝试了两种策略来训练我们的模型。首先,我们使用了基本的自监督学习方案,训练集和测试集来自相同的 ST 数据。结果表明,新的估算斑点不能很好地与簇上的原始斑点合并,导致了类似批量效应的现象(图 4A)。这可能是因为 IDC 数据质量不高,DIST 无法得到很好的训练。因此,我们接下来尝试了迁移学习方案,在质量更好的数据上进行训练。在这里,我们使用 Visium 小鼠后脑矢状面数据[40]来构建训练集。小鼠后脑矢状面数据每个点的平均读数为 87 128 [40],而 IDC 的平均读数仅为 40 795 [29]。估算的聚类结果在视觉上非常平滑,与原始数据的聚类结果一致(图 4B)。这表明,DIST 可以利用从高质量数据中学到的信息,改善低质量数据的估算结果。高质量数据的噪声较小,因此迁移学习可以使 DIST 避免过多噪声,从而提高估算的准确性。以下分析基于迁移学习的结果。

(C) 用 LISI 测量原始(橙色)和估算(蓝色)IDC 在获得平滑和连续空间域方面的聚类性能。LISI 分数越低,表明空间结构越均匀和连续。(D) 原始表达的域。(E)(B)的合并版,与(D)用相同颜色匹配的域。

首先,我们比较了 DIST 归因前后 IDC 的聚类结果(图 4B 和 D 及补充图 S6)。结果表明,DIST 估算表达的聚类结果更为平滑。我们利用卢万算法,使用不同的 PCs 数和分辨率对原始表达和估算表达进行了聚类。从直观上看,在估算的 IDC 上检测到的空间域更连续、更平滑(图 4B 和补充图 S6),较低的局部逆辛普森指数(LISIs)[41]也证实了这一点(图 4C)。在每个参数组合中,基于原始表达的 LISIs 中值都大于基于估算表达的 LISIs 中值。基于原始表达的 LISIs 的最小中值为 1.709,甚至大于基于估算表达的 LISIs 的最大中值 1.510。图 4B 和 D 显示了使用 15 个 PCs 和 0.8 分辨率的聚类结果,该结果被用于进一步的差异表达分析和功能富集分析。

局部逆辛普森指数(LISI):定义了单个细胞附近的数据集的有效数量。LISI值越高,说明邻域的数据集越多,批处理效应越小。要减少批处理效应的影响。

图 5.探索 IDC 数据的生物学功能。(A) 每个域中有意义(调整后 P 值小于 0.05 且 log2-foldchange > 1)的 DE 基因数量。绿色表示来自原始数据的数字;橙色表示来自估算数据的数字;灰色表示两者的重叠。(B) 与单细胞参考文献比较的成对和重叠的显著 DE 基因数量。(C) 病理学家对免疫荧光组织图像的注释(来自文献[18]),其中红色:浸润性癌,黄色:原位癌,绿色:良性增生,灰色:未分类肿瘤。

差异表达分析结果表明,与原始数据相比,从估算表达中发现了更多有意义的标记。为了便于比较原始 IDC 和估算 IDC,我们手动合并了两个域,在两个聚类结果之间形成一对一的映射(图 4D 和 E)。在几乎所有空间域中,估算的 IDC 都能比原始数据发现更多重要的 DE 基因(图 5A)。在每个域中,重叠的 DE 基因数量都接近原始转录组数据,这表明 DIST 可以在保留原始数据信息的同时,帮助识别更多重要的 DE 基因。此外,我们还分析了与之匹配的 IDC 单细胞 RNA-seq 数据[42]作为参考,验证了更多的重要 DE 基因并不是由假阳性引起的(图 5B)。我们认为单细胞 RNA-seq 的 DE 结果更可信,因为它具有较深的测序深度和单细胞分辨率。在这里,我们使用 Seurat 标签转移法将每个细胞映射到空间域中[43]。在癌与免疫域的六次比较中,估算的 IDC 有更多显著的 DE 基因,而且其与单细胞参考的重叠也比原始数据多得多,这表明 DIST 可以发现更多真正的 DE 基因。

根据组织切片的免疫荧光成像和组织病理学注释[18](图 5C),这 12 个域可分为浸润癌(标号:0,4,7,8)、原位癌(标号:3,9)、良性增生(标号:2)和免疫反应相关域(标号:1,5,6,10,11)。在疾病相关领域,从估算数据而非原始数据中的 DE 基因中发现了几个与癌症相关的基因。例如,域 0 中的 EZH2 [44] 和 SOCS1 [45],以及域 4 中的 AKT1 [46],都是众所周知的 IDC 标志 [47]。在免疫相关域 10 中,新发现的 DE 基因包括 ADA,这是一种能促进 CD4+ T 细胞分化和增殖的酶:IGHG1和IGHG4,预计参与激活免疫反应和吞噬作用[47]。

(D) 前 10 个重要的 GO:仅在域 0(上图)和域 10(下图)中通过估算数据确定的 BP 术语。

为了进一步探索生物功能,我们根据 GO: BP 数据库对这些结构域进行了功能富集分析。推算表达带来了一些新发现。例如,域 0 中富集了与染色体组织和 DNA 代谢过程相关的通路,这些通路与癌细胞的生长有关(图 5D)。在域 10 中,包括 B 细胞、T 细胞、白细胞活化和免疫反应正调控在内的多个过程被显著富集(图 5D)。

DIST 对 Visium 人类 IDC 数据进行去噪

转移学习可以提高估算的准确性,但无法降低估算表达的噪声。在本节中,我们将展示如何使用 DIST 通过增加一个额外步骤来对估算表达式进行去噪。

由于单个点内的 RNA 测序量较低,低质量 ST 数据通常具有噪声和稀疏性,导致生物信号被掩盖。对于这些嘈杂和稀疏的数据,DIST 不仅可以通过迁移学习改进估算,还可以对估算表达进行去噪处理,以增加 UMI 计数的深度并降低丢失率。去噪可以通过额外的步骤来实现,即在训练过程中为输入引入合成噪声。LR 输入比其 HR 地面实况有更多的噪声,因此 DIST 可以学习如何同时对低质量数据进行估算和去噪

我们通过与上一节的非变性估算进行比较,显示了基于讨论的 IDC 的降噪效果。我们比较了去噪表达与非去噪表达之间每个点的 UMI 计数总数以及每个基因 UMI 计数为零的点的百分比。结果如图 6A 所示。去噪表达的总计数平均值从 7237 降至 12 367,去噪表达的丢失百分比平均值从 6% 降至 2%。从 UMI 总计数和脱落百分比来看,去噪表达的质量明显提高

图 6.DIST 可同时对低质量的 IDC 进行归因和去噪。(A) 总计数(每个斑点的计数总数)(左)和丢失百分比(每个基因未出现在斑点中的百分比)(右)的小提琴图,分别基于非去重估算的 IDC(浅绿色)和去重估算的 IDC(浅黄色)。(B) 去染估算的合并聚类结果(补充图 S7),除域 1 和域 6 外,与图 4E 相同颜色的域。 (C) 每个域中显著(P 值小于 1×10-6 且 log2-foldchange > 0.6)DE 基因的数量。绿色表示来自非变性数据的数字;橙色表示来自变性数据的数字;灰色表示两者的重叠。

接下来,我们用与非变性数据相同的方法分析了变性数据,包括聚类、差异表达分析和通路富集。去噪数据的聚类结果与非去噪数据相似(图 4E 和 6B)。在每个匹配域中,我们都从去噪数据中发现了比非去噪数据更多的 DE 基因(图 6C)。在疾病相关领域,去噪数据的 DE 基因中包含了几个著名的标记,而非去噪数据则没有。这些基因的例子有域 0 中的 NCOA4 [48] 和 SOCS1 [45],域 4 中的 ADGRB1 [49],它们都是 IDC [47] 的著名标记物;域 3 中的 HOXB13 [50] 和 KLK10 [51],它们是导管原位癌 [47] 的两个标记物。此外,我们还发现了一些富集在上调 DE 基因中的 GO 术语,这些术语只能在去噪数据中检测到,而不能在非去噪数据中检测到。这些富集的术语与癌症和免疫反应的生物学功能高度相关。例如,域 0 中富集了与蛋白质和磷酸盐相关的通路(图 6D)。而与免疫反应相关的多个过程在域 10 中被显著富集(图 6D)。总之,去噪可以降低推算表达的噪声,提高识别重要 DE 基因和生物信号的能力

(D) 前 10 个重要的 GO:BP 术语仅由领域 0(上图)和领域 10(下图)的去噪数据确定。

此外,我们还使用了不同的效率损失 λ 来研究不同的合成噪声水平对去噪性能的影响。我们选择了效率损失 λ = 0.2、0.4、0.6 和 0.8,结果如补充图 S8 所示。从聚类结果来看,效率损失越小,域越平滑,细节越少。因此,效率损失宜大于 0.5,以避免原始细节的缺失。

讨论

我们开发了一种基于深度学习的方法 DIST,该方法通过自我监督学习和迁移学习来增强 ST 数据。通过自我监督学习,DIST 可以准确推算未测量区域的基因表达水平,并在总计数和丢失百分比方面提高数据质量。此外,迁移学习还能让 DIST 借用其他高质量数据的信息来改进推算的基因表达水平。值得注意的是,高质量的信息数据不一定要与目标 ST 数据来自相同的组织或物种,因为我们通过在高质量的小鼠大脑数据上训练的模型增强了质量较差的人类 IDC 数据。我们还使用了不同物种或组织的数据作为来源,对 IDC 数据进行迁移学习,结果表明(补充图 S9),来自同一物种的数据源能带来更好的性能,但差异很小。即使来自不同的物种或组织,DIST 仍能从源数据中学到足够的知识来提高预测性能

在空间分辨转录组学研究中,确定显示空间表达模式的基因是描述复杂组织中 ST 图谱的重要一步 [21]。当试图利用全转录组技术描绘组织内的生物过程和通路时,对具有独特空间模式的基因表达谱特别感兴趣[52]。由于统计能力较低,利用 LR 转录组学数据很难在小规模区域内找到具有空间模式的生物学意义基因。通过分析 ST 人类黑色素瘤数据,我们发现 DIST 使我们能在小区域中检测到更多具有明显空间模式的基因,并使小区域特异性基因(如淋巴和茎基质相关基因)的排名比原始数据高得多。除 LR 外,测序深度较浅或丢弃事件等噪音也是导致探究生物学意义的能力下降的另一个原因。在 Visium IDC 数据分析中,我们发现 DIST 可以通过去噪改善 UMI 总计数并减少丢弃,使用 DIST 增强数据,我们可以发现更多 DE 基因,从而找到更多重要的富集通路。

为了量化 DIST 的时间和内存成本,我们记录了 NN、DIST、Linear、Cubic、NEDI、SpatialPCA 和 XFuse 在相同 ST 数据上的时间和内存使用情况。结果如补充表 S1 所示。DIST 是第二快的方法,只需要小于 2 分钟就能完成整个过程。

虽然我们的评估是在 ST 和 Visium 平台上进行的,但 DIST 也应适用于其他平台,如 Stereo-seq[17]。我们将 DIST 应用于 E9.5 小鼠胚胎的立体序列数据,结果如补充图 S10 所示。与 IDC 数据相同,DIST 显著提高了总计数,并发现了比原始数据更多的重要 DE 基因。此外,我们的工作主要针对二维组织切片及其基因表达,而 DIST 可以扩展到三维(3D)ST 重构,为生物探索提供更全面的视角。为了重建三维表达,现有方法通常会对多个相邻组织切片的 ST 数据进行配准和整合 [53]。由于组织切片之间没有匹配的组织学图像,一些归因方法无法轻松应用于三维表达。然而,DIST 只需要基因表达和斑点坐标,因此只需对网络参数和 VLN 输入的维度稍作修改,就能应用于三维表达。总之,DIST 提供了全面的数据增强功能,包括归因和去噪,可以成为 ST 数据分析的有用工具。

要点

- DIST 提供全面的数据增强功能,包括估算和去噪,可作为空间转录组(ST)数据分析的有用工具。

- DIST 利用自学增强 ST 数据,不需要组织学图像等任何附加信息。

- DIST 利用迁移学习从其他高质量数据中借用信息,从而显著提高了 ST 数据的质量。

- DIST 在整个基因集上进行训练,利用了所有基因和斑点的优势,而不是只利用一个基因和周围的斑点。

- DIST 提高了聚类的准确性,有助于为 ST 数据识别出更多有生物学意义的差异表达基因和通路,从而更深入地了解生物学过程。

Logo

DAMO开发者矩阵,由阿里巴巴达摩院和中国互联网协会联合发起,致力于探讨最前沿的技术趋势与应用成果,搭建高质量的交流与分享平台,推动技术创新与产业应用链接,围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐