整合单细胞转录组数据的空间转录组解卷积

Spatial Transcriptome Deconvolution Integrating Single-Cell Transcriptome Data

摘要: 空间转录组技术是理解组织结构和功能的关键。然而，许多空间转录组技术没有单细胞分辨率，其基因组数据反映了细胞群的平均表达水平，在此基础上的差异、聚类和关联分析等数据分析问题实际是在细胞群水平上进行的，忽略了不同细胞类型之间的差异，导致分析结果出现偏差甚至错误。本文开发了一种解卷积算法stDRPC，该算法首先从单细胞转录组数据中得到细胞类型特异性表达信息矩阵，将该矩阵作为参考基矩阵，然后利用稳健偏相关性算法对空间转录组数据进行基于参考基矩阵的解卷积，从而得到各细胞类型所占比例。小鼠脑皮层及人胰腺导管腺癌(PDAC)数据分析表明所提算法在准确性及生物学意义上都优于常用算法。

Abstract: Spatial transcriptome technology is crucial for understanding tissue structure and function. However, many spatial transcriptome techniques do not have single-cell resolution, and their genomic data reflects the average expression level of cell populations. On this basis, the data analysis problems such as differential, clustering and association analysis are actually carried out at the cell population level, ignoring the differences among different cell types, resulting in biases and errors. This article develops a deconvolution algorithm called stDRPC, which first obtains cell type specific expression information matrix from single-cell transcriptome data, uses this matrix as a reference-based matrix, and then uses a robust partial correlation algorithm to deconvolve the spatial transcriptome data based on the reference-based matrix to obtain cell type proportions. Data analyses of mouse cortical and human pancreatic ductal adenocarcinoma (PDAC) showed that the proposed algorithm is superior to commonly used algorithms in terms of accuracy and biological significance.

1. 引言

空间转录组技术(spatial transcriptomics technologies)使我们能够利用空间定位信息对许多组织位置进行基因表达谱分析，从而能够表征空间组织结构[1]。然而，大多数技术(如基于测序的技术)的空间分辨率有限，它们收集组织位置的基因表达测量值是细胞类型群体的混合物[2]-[4]。因此，对每个空间位置进行细胞类型解卷积成为解开细胞类型空间定位和表征复杂组织结构的重要任务[5]。

借鉴Bulk RNA-seq的解卷积方法，并考虑到空间转录组数据的特点，国内外许多研究者提出了一系列空间转录组数据的解卷积方法，根据算法的建模原理，这些方法可分为两类：基于概率模型的方法，主要包括stereoscope [6]、cell2location [7]、RCTD [8]、CARD [9]、MuSic [10]、Seurat [11]等。其中stereoscope [6]利用相同组织的scRNA-seq作为参考信息，在计数数据均服从负二项分布的假设下构建了单细胞表达谱映射到空间转录组数据的概率模型，然后使用最大后验估计求解模型参数。cell2location [7]构建了层次贝叶斯模型，使用负二项回归求解空间转录组数据中的各种细分细胞类型。RCTD [8]同样借助注释后的单细胞数据作为参考，建立了空间转录组数据关于单细胞转录组数据的“泊松–对数正态混合模型”，以考虑两种数据类型的平台差异。CARD [9]通过结合单细胞RNA测序(scRNA-seq)数据和空间相关性信息，提高了组织中细胞类型解卷积的准确性。即使在参考数据不匹配的情况下，它也能精确识别细胞类型，并能够生成高分辨率的组织空间图。CARD计算效率高，结果比其他方法更加准确，已被应用于多种数据集的分析，包括癌症研究。MuSiC [10]通过权重基因的跨主体和跨细胞一致性，实现了更精确的细胞类型比例估计。研究表明，该方法在具有密切相关细胞类型的组织中表现出色，且其准确性优于其他现有方法，广泛应用于胰腺和肾脏等器官的研究，帮助揭示疾病机制。Seurat [11]通过字典学习将不同模态的数据转换为共享的特征空间，使得不同模态的数据可以在一个统一的框架内对齐，从而提高计算效率并支持大规模的单细胞数据集整合。文章展示了该方法在基因表达、染色质可及性、蛋白质水平等多种模态上的成功应用。研究表明，该方法在具有密切相关细胞类型的组织中表现出色，且其准确性优于其他现有方法，广泛应用于胰腺和肾脏等器官的研究，帮助揭示疾病机制。基于非负矩阵分解和最小二乘模型的方法，例如spatialDWLS [12]和SPOTlight [13]。其中spatialDWLS [12]首先利用富集分析识别每个spot中可能包含的细胞类型，然后利用改进的阻尼加权最小二乘方法(DWLS)方法推断每个spot中细胞类型的比例。SPOTlight [13]利用NMF识别scRNA-seq参考数据中的细胞类型特异性主题表达谱，然后利用NNLS算法识别spot主题表达谱，得到细胞类型比例矩阵。

本文将空间转录组数据与单细胞转录组数据相整合，构建统计模型。首先从单细胞转录组数据中得到细胞类型特异性表达信息矩阵，将该矩阵作为参考基矩阵，然后利用稳健偏相关性算法对空间转录组数据进行基于有参考基矩阵的解卷积，从而得到各细胞类型的比例矩阵，并通过小鼠脑皮层和人胰腺导管腺癌(PDAC)数据分析来表明所提算法的优越性。

2. 整合单细胞转录组数据的空间转录组解卷积模型(stDRPC)构建

2.1. stDRPC算法概述

stDPRC算法的输入为：同一组织的空间转录组数据和单细胞转录组数据(scRNA-seq)，输出为：每个空间位置(spot)上每个细胞类型所占比例。

stDPRC算法中scRNA-seq数据是参考数据，假设其由K个细胞类型和G个细胞类型特异性的基因组成。X表示N个空间位置上测量的同一组G个信息基因上的基因表达矩阵。W表示细胞类型特异性表达矩阵，维数为，其通常被称为参考基矩阵，可由scRNA-seq计算求出。stDPRC算法构建的统计模型为：

.(2.1)

其中，H表示维数为的细胞类型比例矩阵，每个元素非负，H的每一行表示每个spot上K个细胞类型的比例，行和为1。E是一个的残差矩阵，每个元素独立且服从正态分布。我们的目标是根据已知的X和求得的W来估计H。

Figure 1. Flowchart of the stDPRC algorithm

图1. stDPRC算法流程图

图1显示了stDRPC的工作流程，其中包括两个主要步骤：第一步是利用scRNA-seq来构建参考基矩阵W，第二步是将空间转录组矩阵X与参考基矩阵W整合，利用稳健偏相关算法进行解卷积，求解细胞组成比例矩阵H。具体步骤如下。

2.2 利用已知的scRNA-seq数据求解参考基矩阵W

假设某组织的scRNA-seq数据由K个细胞类型组成，受试者j在基因g上mRNA分子的总数为，那么，其中表示受试者j的细胞c中基因g的mRNA分子的数量，是样本j中属于细胞类型k的细胞索引集，表示该集合中的细胞总数()。受试者j中基因g的相对丰度为

,(2.2)

可以计算得到

$S{c}_{jg}={\displaystyle {\sum }_{k=1}^K{m}_j^k{S}_j^k{\theta }_{jg}^k=m_j{\displaystyle {\sum }_{k=1}^K{p}_j^k{S}_j^k{\theta }_{jg}^k}}$ .(2.3)

其中，$S_j^k={\displaystyle {\sum }_{c\in C_j^k}{\displaystyle {\sum }_{g^{\prime }}^{G}S{c}_{j{g}^{\prime }c}}}/m_j^k$ 是细胞类型K的全部细胞的平均量，$m_j={\displaystyle {\sum }_{k=1}^K{m}_j^k}$ 是组织中mRNA分子的总数，$p_j^k=m_j^k/m_j$ 是细胞类型k的比例。假设$Y_{jg}=S{c}_{jg}/{\displaystyle {\sum }_{g^{\prime }}^{G}S{c}_{j{g}^{\prime }}}$ 为样本j中基因g的相对丰度，则

$Y_{jg}=S{c}_{jg}/{\displaystyle {\sum }_{g^{\prime }=1}^{G}S{c}_{j{g}^{\prime }}}$ .(2.4)

因此，下式成立

$\left[\begin{array}{c}{Y}_{j1}\\ ⋮\\ Y_{jG}\end{array}\right]\infty \left[\begin{array}{ccc}{\theta }_{j1}^1& \cdots & {\theta }_{j1}^k\\ ⋮& \ddots & ⋮\\ {\theta }_{jG}^1& \dots & {\theta }_{jG}^k\end{array}\right]\cdot \left[\begin{array}{ccc}{S}_j^1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & S_j^k\end{array}\right]\cdot \left[\begin{array}{c}{P}_j^1\\ ⋮\\ P_j^k\end{array}\right]$ .(2.5)

根据已知的scRNA-seq数据信息，求得细胞大小因子$S_j^k$ 和细胞类型特异性相对丰度${\theta }_{jg}^k$ ，细胞类型特异性参考基矩阵W为：

$W=\left[\begin{array}{ccc}{\theta }_{j1}^1& \cdots & {\theta }_{j1}^k\\ ⋮& \ddots & ⋮\\ {\theta }_{jG}^1& \dots & {\theta }_{jG}^k\end{array}\right]\cdot \left[\begin{array}{ccc}{S}_j^1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & S_j^k\end{array}\right]$ .(2.6)

2.3. 将已知的X和求得的W相结合，利用稳健偏相关算法进行解卷积

首先，对已知的空间转录组数据X进行预处理，去掉在所有spot上表达量低的基因，对剩余的基因在spot上的矩阵进行处理，去掉批次效应。将参考基矩阵W和空间转录组矩阵X进行匹配，使其有共同的基因。其次，对每个spot，将空间转录组矩阵X在该spot上的表达量作为因变量，将参考基矩阵W作为自变量，利用Robust Linear Model (RLM)进行解卷积(rlm函数)，将回归模型输出的系数进行标准化处理，使其值大于0，且和为1。最后，将每个spot上的估计值作为行构成矩阵，该矩阵即为最终所求细胞类型比例矩阵H。

2.4. 算法评价指标

本文选取皮尔逊相关系数(PCC)、结构相似性指数(SSIM)、Jensen-Shannon散度(JSD)、均方根误差(RMSE)四个指标来评价算法的性能。高的PCC和SSIM，低的JSD和RMSE表明算法有更高的准确性。四个指标具体计算公式如下：

(1) PCC的计算公式为：

$\text{PCC}=\frac{{\displaystyle \sum \left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}}{\sqrt{{\displaystyle \sum {\left(X_i-\overline{X}\right)}^2}{\displaystyle \sum {\left(Y_i-\overline{Y}\right)}^2}}},$(2.7)

其中，$X_i$ 和$Y_i$ 分别表示真实比例和估计比例构成的向量，$\overline{X}$ 和$\overline{Y}$ 分别是它们对应的均值，PCC值越高，表明算法精确度越高。

(2) SSIM指标可用来衡量结构相似性，计算公式如下：

$\text{SSIM}\left(X,Y\right)=\frac{\left(2{\mu }_X{\mu }_Y+C_1\right)\left(2{\sigma }_{XY}+C_2\right)}{\left({\mu }_X^2+{\mu }_Y^2+C_1\right)\left({\sigma }_X^2+{\sigma }_Y^2+C_2\right)},$ (2.8)

其中X和Y分别为真实比例和估计比例构成的向量；$u_X$ 和$u_Y$ 分别为真实比例和估计比例的均值；${\sigma }_X^2$ 和${\sigma }_Y^2$ 分别为真实比例和估计比例的方差；${\sigma }_{XY}$ 是两者的协方差；常数$C_1$ 和$C_2$ 用于避免分母为零，通常取小值(例如0.01和0.03)。它衡量了细胞类型真实比例与估计比例在局部结构上的相似性。SSIM值越高，说明估计比例的空间分布特征与真实比例越接近，算法的准确性越好。

(3) JSD是基于信息熵的距离度量，用来衡量真实比例和估计比例之间的概率分布差异。公式为：

$\text{JSD}\left(P,Q\right)=\frac{1}{2}{D}_{\text{KL}}\left(P\left|\right|M\right)+\frac{1}{2}{D}_{\text{KL}}\left(Q\left|\right|M\right),$ (2.9)

其中，P和Q分别为真实比例和估计比例的概率分布，$M=1/2\left(P+Q\right)$ ，而$D_{KL}$ 是Kullback-Leibler散度。JSD值越小，表示真实比例与估计比例的概率分布越接近，算法的准确性越好。

(4) RMSE用于评估估计比例与真实比例之间的总体误差。计算公式为：

$\text{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(Y_i-{\widehat{Y}}_i\right)}^2}},$(2.10)

其中$Y_i$ 为真实比例，${\widehat{Y}}_i$ 为估计比例，N是样本数。RMSE值越小，算法的准确性越好。

3. 小鼠脑皮层及人胰腺导管腺癌数据实证分析

3.1. 数据介绍

小鼠脑皮层数据。我们从GEO数据库获取了使用Smart-seq技术生成的scRNA-seq数据和seqFISH+的小鼠脑皮层空间转录组数据。小鼠脑皮层的scRNA-seq数据共含有14,249个细胞，我们去除了RNA数少于200个的细胞。小鼠脑皮层的空间转录组数据含有524个spot，根据这524个spot的位置信息，我们绘制了网格图，每个网格包含1到18个细胞，最终得到一个含有72个spot的空间转录组数据。由于每个网格(即spot)中包含几个细胞类型、每个细胞属于哪个细类型都是已知的。因此，该spot中某个细胞类型所占比例可用属于该细胞类型的所有细胞数量之和除以所有细胞数量之和求得。即在小鼠脑皮层的空间转录组数据上，细胞类型组成比例矩阵H是已知的。因此，该数据可用作金标准来评价算法的性能[14]。

人胰腺导管腺癌数据。我们从GEO数据库下载胰腺导管腺癌(PDAC)空间转录组数据和scRNA-seq数据作为stDRPC算法的输入。PDAC的空间转录组数据含有428个spot和25,753个基因。scRNA-seq数据含有1926个细胞。由于在该数据上细胞类型比例矩阵H是未知的，因此对该数据我们从生物学意义方面来进行分析。

RCTD、CARD、SpatialDWLS、MuSic、Seurat是经典的空间转录组解卷积算法。因此，下文我们将stDRPC与这5种算法进行比较。对小鼠脑皮层数据，用PCC、JSD、SSIM、RMSE这4个指标来评价算法的性能[15] [16]。对人胰腺导管腺癌数据，从生物学意义方面来进行对比分析。

3.2. 小鼠脑皮层数据实证分析

小鼠脑皮层数据共有13种细胞类型，进行网格化后，每一种细胞都被唯一细胞类型注释，并且落在唯一的一个网格内。根据图2(A)，最终得到一个含有72个spot的空间转录组数据。由于该数据上细胞类型比例矩阵H是已知的，因此我们将stDRPC算法与RCTD、CARD、SpatialDWLS、MuSic、Seurat进行比较，图2(B)显示了6种算法在4个指标下的结果。stDRPC估计的细胞类型所占比例与真实的细胞类型比例之间的PCC值高达0.86，而RCTD为0.853、CARD为0.857、spatial DWLS为0.864、MuSic为0.89、Seurat为0.83。stDRPC的SSIM为0.68，仅次于MuSic，但比其它4种方法的值都大。stDRPC的JSD为0.16，仅高于MuSic，但比其它4种都小。stDRPC的RMSE为0.10，仅高于MuSic，但比其它4种都小。显然，整体来看stDPRC在四个指标上都取得了很好的准确性。

Figure 2. (A) Grid plots of mouse cortex datasets; (B) Performance comparison of six algorithms under four indicators

图2. (A) 小鼠脑皮层数据的网格化图；(B) 四种指标下六种算法的性能比较

3.3. 人胰腺导管腺癌数据实证分析

图3(A)展示了使用stDRPC与其他五种常用反卷积方法(Seurat、MuSiC、RCTD、SpatialDWLS和CARD)估计的细胞类型比例的相关性分析结果。横坐标为不同的反卷积方法，纵坐标显示了与匹配的单细胞RNA测序(scRNA-seq)数据中的细胞类型比例的相关系数。该相关系数用于衡量不同方法估计的空间转录组数据与真实单细胞测序数据中细胞类型比例的一致性。结果表明，stDRPC方法在相关性分析中表现最佳，相关系数达到0.87，显著优于其他反卷积算法，如CARD (0.79)、SpatialDWLS (0.55)、RCTD (0.53)、MuSiC (0.16)和Seurat (0.13)。这表明stDRPC在反卷积任务中的性能更为优越，能够更准确地预测不同细胞类型在空间转录组数据中的分布。此外，图中使用的PDAC (胰腺导管腺癌)数据集包含多个组织区域，包括癌症、胰腺、导管和间质区域。通过将预测的细胞类型分布与基于H&E染色的组织学家注释对比，stDRPC能够准确地将胰腺和肿瘤细胞定位到相应的组织区域，从而验证了其空间定位能力。

图3(B)展示了六种算法在每个空间位置上20种细胞类型比例组成的空间散点饼图。根据stDRPC估计的细胞类型组成矩阵的第一个主成分(优势细胞类型)，可以清楚地将癌症和非癌症区域、导管和间质区域以及胰腺和导管区域之间的大体区域分离，层次清晰。相比其他方法，效果更为显著。

Figure 3. Results of real data deconvolution for PDAC: (A) Comparison of average correlation coefficients between RPCD and five other deconvolution methods. (B) Spatial scatter pie charts of cell type composition at each spatial location for six methods. (C) Expression of specific marker genes for selected cell types deconvoluted by RPCD across six cell types. (D) Proportion of each cell type under cancerous and normal cell types

图3. PDAC的真实数据反卷积结果：(A) RPCD与五种反卷积方法平均相关系数比较。(B) 六种方法每个空间位置上细胞类型组成的空间散点饼图。(C) RPCD反卷积出六种细胞类型，选定细胞类型的特异性标记基因的表达量。(D) 每种细胞在癌症细胞类型和正常细胞类型下的比例

图3(C)第一行分别显示stDRPC估计出的六个细胞类型(Cancer clone A、Cancer clone B、Ductal terminal、Ductal centroacinar、Ductal MHC Class II、Ductal high hypoxic)在所有spot上估计出的比例，第二行显示对应细胞类型特异性标记基因的表达水平。显然，上亚区由癌症克隆A细胞主导，标记基因TM4SF1富集。下亚区由标记基因S100A4富集的癌症克隆B细胞主导。S100A4是早期胰腺癌的预后标志物，其空间富集提示癌亚区域底部可能是早期癌症区域。TM4SF1在stDRPC的迁移和侵袭中扮演重要角色，其空间富集提示癌上亚区可能是具有转移能力的晚期癌症区域。事实上，stDRPC也检测到癌症上部亚区域富含成纤维细胞和成纤维细胞标记基因CD248，成纤维细胞与晚期肿瘤淋巴结转移密切相关。

由于癌症和非癌症区域的分离，图3(D)展示了采用双侧Wilcoxon秩和检验比较stDRPC推断的癌症区域(n = 273)与非癌症区域(n = 155)的细胞类型比例，我们可以看出，Cancer clone A和Cancer clone B细胞高度富集于癌症区(秩和检验得到$p=2\times {10}^{-6}$ )。stDRPC还揭示了两个巨噬细胞亚群在癌症区域和非癌症区域之间的不同分布，这是癌症组织区域化的一个关键功能特征。

stDRPC可将细胞类型定位到特定组织区域，该结果与细胞类型相应标记基因的表达模式一致。stDRPC捕获了导管中心腺泡和末端导管细胞，图4(A)显示这两种细胞无论是在肿瘤区域还是在非肿瘤区域都广泛存在。此前，没有任何一种方法能够捕获导管中心腺泡和末端导管细胞的预期空间定位，其他算法推断腺泡细胞要么在胰腺区缺失，要么从胰腺区向外扩散到间质区和癌区。但stDRPC推断出腺泡细胞主要富集于正常胰腺组织区。

Figure 4. (A) stDRPC captures ductal centroacinar cells and terminal duct cells; (B) Correlation of cell type proportions at different spatial locations inferred by stDRPC

图4. (A) stDRPC捕获了导管中心腺泡和末端导管细胞；(B) stDRPC推断的不同空间位置间细胞类型比例的相关性

图4(B)展示了stDRPC推断的细胞类型之间在不同空间位置上的细胞类型比例相关性矩阵。矩阵的每一行和每一列分别代表不同的细胞类型，颜色反映了这些细胞类型之间的相关系数(右侧的颜色条表示了相关系数的范围，红色为正相关，绿色为负相关)。相关系数用于评估在空间转录组数据中，不同细胞类型在相邻或相似的空间区域内是否存在相互关系或共定位趋势。从图中可以看到，细胞类型之间的相关性存在显著的差异。例如，导管高缺氧细胞与癌症细胞(如“Ductal high hypoxic”和“Cancer clone A/B”)之间表现出较强的正相关性，表明它们可能共定位于同一肿瘤微环境中，形成缺氧和营养不良的肿瘤环境。类似地，内皮细胞(“Endothelial cells”)与成纤维细胞(“Fibroblasts”)也表现出较强的相关性，暗示它们在胰腺癌的基质区域中可能存在相互作用，协同促进肿瘤进展。这个相关性矩阵基于stDRPC推断出的细胞类型空间分布，揭示了不同细胞类型在胰腺癌(PDAC)组织中的空间分布及其相互作用关系[17]。

4. 结论

本文针对空间转录组数据的解卷积问题，开发了统计算法stDRPC，该算法将单细胞转录组数据与空间转录组数据相融合，从单细胞转录组数据中得到细胞类型特异性表达信息矩阵，利用稳健偏相关性算法对空间转录组数据进行基于参考基矩阵的分解，从而得到各细胞类型所占比例。小鼠脑皮层数据分析表明所提算法在准确性上优于常用经典算法。人类胰腺导管腺癌(PDAC)真实肿瘤样本分析发现，stDPRC能准确定位不同细胞类型在不同组织区域的分布，这有助于识别不同空间定位的多种细胞类型和分子标记。

基金项目

研究得到了江西省自然科学基金(20212BAB202001)的支持。

NOTES

^*通讯作者。