文献阅读-《Benchmarking spatial and s

2022-07-02 本文已影响0人 whitebird

文献名称：用于转录本空间分布预测和细胞类型反卷积的单细胞转录组和单细胞空间组学的多种整合方法的基准测试
文献链接：https://www.nature.com/articles/s41592-022-01480-9
文献代码：https://github.com/QuKunLab/SpatialBenchmarking
参考资料：https://www.zhihu.com/column/c_1340036890003693568

摘要

空间转录组技术大大提高了我们检测组织中RNA转录本空间分布的能力，但表征空间中单个细胞的全转录组水平仍然具有挑战性。为了满足这一需求，研究人员开发了整合方法，将空间转录组数据与单细胞RNA-seq数据相结合，以预测未检测到的基因的空间表达分布和/或对组织切片中的spot进行细胞类型去卷积。然而，迄今为止，还没有专门的研究对这些整合方法进行基准分析来衡量它们的性能。本文我们使用45个配对数据集（包括空间转录组数据和 scRNA-seq 数据）和 32 个模拟数据集对 16 种整合方法进行基准测试。我们发现Tangram、gimVI 和 SpaGE在预测RNA转录物的空间分布方面优于其他整合方法，而Cell2location、SpatialDWLS 和 RCTD是用于spot细胞类型去卷积的最佳方法。我们提供了一个基准分析流程来帮助研究人员选择最佳的整合方法来处理他们的数据集。

问题1：为什么要进行空间转录组数据与单细胞RNA-seq数据进行整合？
我的理解：他山之石，可以攻玉。单细胞转录组数据包含每个细胞全基因组信息，但丢失了细胞的空间位置信息。通过单细胞转录组技术，来弥补当前空间转录组技术的缺陷；

基于图像的空间转录组技术，优点是空间单细胞分辨率，缺点是检测基因的数目有限，利用配对的单细胞转录组信息，可以预测未检测基因在组织上的空间表达分布；“弥补”检测基因数目的限制；

2.基于测序技术的空间转录组技术，优点是检测全转录组（基因数目不受限），缺点是空间上不是单细胞分辨率，spot包含多个细胞，通过配对的单细胞转录组数据，对组织切片中的spot进行细胞类型去卷积，预测spot细分细胞类型；

空转的问题和解谜：
实验技术是我们探索生命世界的工具，算法工具都是在弥补实验技术的不足，或者说弥补不到，尽量去找近似解。算法是技术的补充和延伸。
一个新的实验技术的产生，它能解决重要的问题，比如说10X空转，它可以捕获组织spot的空间位置信息，也知道spot的转录表达的信息。
问题：
但这也有一个致命的问题，相比单细胞转录组技术，空转技术，为了捕获空间的信息，丢掉了单细胞分辨率的特性，一个spot包含多个细胞（1-10），也就是一个spot是多细胞的混杂，也是所有下游分析，多样本整合，找差异基因等都绕不开的问题。
解题：
有问题就会有吸引很多解谜者。空转技术没有实质性突破，空转+单细胞转录组联合的算法研究将持续成为热点。

问题2：怎么理解45 个配对数据集（包括空间转录组学和 scRNA-seq 数据）？
答：怎么理解配对数据？什么样的数据可以称为配对数据？
组织样本上，空间位置尽量的近，一部分样本做单细胞转录组，一部分做单细胞空间转录组；相当于彼此的近似解，细胞类型的组成应该是一致的；
粗略找一两篇篇文献具体看下，有两种类型：
1.空转样本+单细胞转录组公共数据集；

2.组织空转样本+组织单细胞样本

问题3：16种单细胞转录组+空间转录组的整合方法的介绍
答：1.预测RNA转录物的空间分布的算法有8种；2.用于spot细胞类型去卷积的算法有12种；其中4个算法(Tangram, Seurat, SpaOTsc, novoSpaRc)是通用的，均适用于两种功能；总共16种整合分析方法。

前言

空间转录组技术使我们能够检测样本空间中的RNA转录信息，这些方法已被用于研究各种组织和器官中基因在空间上的表达分布，包括大脑、心脏、胰腺和皮肤。

空间转录组技术包括以下两种：一种基于原位杂交（situ hybridization）和基于图像荧光显微镜（fluorescence microscopy）的空间转录组学方法——包括 seqFISH、osmFISH6和 MERFISH——以高分辨率和高准确度检测基因的空间表达分布，但是它们在检测总基因的数目是有限的。另一种，基于二代测序（seq-based）的空间转录组技术，例如 ST、10X Visium 和 Slide-seq，可以从空间中的spot捕获全转录组规模的基因表达水平，但每个spot（半径 10-100 μm）可能包含多个细胞，这限制了这些方法的空间分辨率。这些空间转录组技术的局限性阻碍了它们在空间中以单细胞分辨率捕获全转录组规模数据的能力。

问题4：单细胞空间转录组的spot孔径越小，会带来哪些弊端？
答：我的理解是：基于测序的空间转录组技术，spot孔径越小，也不等价于单细胞分辨率；spot跟单个细胞(cell)还是不能一一匹配，另外，spot越小会导致那些跨spot的细胞越多，这类细胞应该不利于下游分析；另外，实验端也是有问题的，切片会弥漫，孔径越小，弥漫程度越高，实验不好操作；

预测基因的空间表达分布算法
为了突破空间转录组技术的局限性，生物信息学家提出并开发了各种整合方法来结合空间转录组学和单细胞 RNA-seq (scRNA-seq) 数据。例如，gimVI 采用深度整合模型（deep generative model）来推断未检测到基因的可能空间表达分布； SpaGE使用域适应算法PRECISE 和 k-最近邻回归来预测未检测到基因的空间表达分布； Tangram使用非凸优化和深度学习框架来学习scRNA-seq 数据的空间对齐； Seurat应用典型相关分析（cca）将空间和 scRNA-seq 数据嵌入到一个共同的潜在空间中，并将 scRNA-seq 数据中的细胞映射到空间转录组数据的spot上； LIGER 使用综合非负矩阵分解和共享因子邻域图来预测空间中的基因表达水平； novoSpaRc和 SpaOTsc 各自使用最优传输方法在 scRNA-seq 数据的基础上构建细胞的空间度量； stPlus结合了自动编码器和加权 k 最近邻方法来预测空间基因表达。这些整合方法使研究人员能够预测未检测到基因的空间表达分布。

预测组织切片中spot的细胞类型组成方法
此外，Seurat、Tangram、novoSpaRc 和 SpaOTsc 能够将 scRNA-seq 数据中的细胞分配到组织切片中的空间位置；这对于提高使用空间转录组技术（如 ST 或 10X Visium）生成的空间转录组数据的分辨率非常有用。此外，Cell2location使用 scRNA-seq 数据中细胞亚群的基因表达特征来估计每个spot的每种细胞类型的丰度； RCTD应用从scRNA-seq数据获得的细胞类型图谱，通过监督学习来分解混合的细胞类型； SpatialDWLS采用加权最小二乘法推断细胞类型组成； Stereoscope利用基于模型的概率方法和 scRNA-seq 数据对空间数据中的细胞混合物进行去卷积；SPOTlight 将种子非负矩阵分解应用于spot的去卷积；DSTG使用基于图的卷积网络对空间转录组数据进行去卷积；STRIDE使用从 scRNA-seq 数据训练的主题轮廓从空间混合物中分解细胞类型； DestVI采用变分推理和潜变量模型来描绘细胞类型比例。这些整合方法使研究人员能够预测组织切片中spot的细胞类型组成。

整合方法的评估
这些整合方法的出现无疑加深了我们对空间转录组数据及相关生物学和病理过程的理解。然而，据我们所知，没有专门的研究全面评估这些整合方法在预测转录本空间分布或组织切片中spot的细胞类型去卷积方面的性能。在这里，我们使用多种指标系统地对16种整合方法的性能进行基准测试，这些方法可以预测未检测到的转录本的空间分布，或组织切片中spot的细胞类型组成（图1a），基于对45个包含空间转录组学数据和scRNA序列数据的配对数据集以及32个模拟数据集的处理（图1b）。我们评估了每种整合方法在预测转录本空间分布上的准确性，包括从原始数据集中下采样(抽样)的稀疏空间转录组数据(sparse spatial transcriptomics data)。我们还根据数据集的模拟评估组织切片中spot的细胞类型去卷积积分方法的准确性，其中每个spot可以包含多个不同类型的细胞。最后，我们评估了每种整合方法消耗的计算资源。我们的研究结果可以帮助研究人员为他们的数据集选择合适的整合方法，另外，这也提出了一些有趣的问题，即各种处理和数据集特定属性是如何影响这些工具在空间转录组学研究中的整合性能。

图 1 | 基准测试框架和配对测试数据集的重要特征；a-用于比较配对空间转录组学和 scRNA-seq 数据集的整合方法性能的基准测试分析流程示意图。我们使用 16 种整合方法来结合空间和单细胞转录组学数据，然后比较它们在 (1) 预测 RNA 转录本的空间分布； (2) 组织切片spot的去卷积后细胞类型；这些方面的性能。我们还评估了整合方法消耗的计算资源。

b-本研究中使用的 45 个配对数据集和 32 个模拟数据集的信息：每个数据集都包含同一组织的空间转录组数据和 scRNA-seq 数据。数据源的详细信息见方法和补充表 1。

问题：13种空间转录组技术的特点和应用场景？、

问题：32 个模拟数据集是如何生成的？
可具体查看SimulatedData.ipynb的代码

def generate_dataset_data(scRNADir1, scRNADir2, celltype_key,reverse = False):
    scrna, spatial = check_dataset(scRNADir1, scRNADir2, top_genes=200, plot=False,reverse=reverse,check=False)
    if not reverse:
        output_dir = scRNADir1
    else:
        output_dir = scRNADir2
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
    generate_spatial_data(spatial,celltype_key,5,3,3,1.5,save_path = output_dir + '/SimulatedSpatial.h5ad')
    if len(scrna.obs) > 10000:
        scrna = scrna[np.random.permutation(scrna.obs_names)[:10000]].copy()
    scrna.obs = scrna.obs.loc[:,['celltype_final']]
    celltype_counts = scrna.obs.celltype_final.value_counts()
    celltype_drop = celltype_counts.index[celltype_counts < 2]
    print(f'Drop celltype {list(celltype_drop)} contain less 2 sample')
    scrna = scrna[~scrna.obs.celltype_final.isin(celltype_drop),]
    sc.pp.normalize_total(scrna, target_sum=1e4)
    sc.pp.log1p(scrna)
    sc.tl.rank_genes_groups(scrna,groupby='celltype_final')
    scrna.write_h5ad(output_dir + '/Simulated_scRNA.h5ad')

结论

1-基准测试框架和测试数据集

为了评估 16 种整合方法的性能，我们从已发表的文献研究中收集了45组配对的空间转录组和scRNA-seq数据集（图 1 和补充表 1）。空间转录组数据集由13种空间转录组技术产生，包括FISH、osmFISH、seqFISH、MERFISH、STARmap、ISS、EXseq、BaristaSeq、ST、10X Visium、Slide seq、seq scope和HDST，scRNA seq数据集由Drop-seq、Smart-seq和10X Chromium平台获得。我们设计了一个分析流程来评估整合空间转库组和单细胞转录组数据集方法的性能（图1a）。在scRNA-seq数据集的预处理过程中，我们剔除了表达基因少于200 个的细胞。对于空间转录组数据集，我们使用 2 个标准生成了“ground truth”：对于检测到的基因数目少于 1,000 个的样本，我们使用了所有的基因；对于检测到的基因数目超过1,000 个的样品，我们使用 1,000 个高可变基因（根据每个基因表达的变异系数进行评估；详见Methods）。

此外，我们采用了RCTD 和 Stereoscope文献中提出的模拟数据方法，并从16 个配对的 scRNA-seq 数据集（补充表 2 和 3）中生成了 32 个模拟的 10X Visium 数据集。一个模拟spot包含从 scRNA-seq 数据集（方法）中随机采样的 5-15 个细胞，每个spot的基因表达值为该spot中所有细胞的总和。

收集数据集后，我们首先评估了8种整合方法的性能，包括 Tangram、gimVI、SpaGE、Seurat、SpaOTsc、novoSpaRc、LIGER 和 stPlus，在预测空间转录组数据集中未检测到基因的空间表达分布。我们使用收集的45个配对数据集，来评估这些整合方法预测 RNA空间分布的准确性。然后，我们对空间转录组数据进行下采样，以测试具有稀疏表达矩阵的数据集的整合方法的性能。

除了预测基因的空间表达分布外，Tangram、Seurat、SpaOTsc 和 novoSpaRc 还可以将 scRNA-seq 数据中的细胞分配到组织切片中的空间位置。此外，通过结合空间转录组数据和 scRNA-seq 数据，使用 Cell2location、SpatialDWLS、RCTD、Stereoscope、DestVI、STRIDE、SPOTlight 和 DSTG 来预测组织切片中spot的细胞类型组成。所有 12 种整合方法都能够对使用 10X Visium 或 ST 平台生成的空间转录组数据集中的spot细胞类型进行去卷积。为了比较这些整合方法在细胞类型去卷积中的性能，我们使用数据集4和10作为基础来模拟代表低空间分辨率数据集的“网格”，并且我们从 scRNA-seq 数据中模拟了32个数据集作为ground truth （方法）。简而言之，在模拟的低分辨率数据集中，每个网格“spot”包含 1-18 个细胞，类似于 10X Visium 或 ST 方法生成的空间转录组数据集。最后，我们评估了每种整合方法消耗的计算资源。

2-预测基因的空间表达分布的算法

我们对45个配对数据集使用十倍交叉验证（详见方法）来评估每种整合方法在预测RNA转录本空间分布方面的准确性。我们通过计算空间转录组数据集的ground truth中基因的表达量与每种整合方法预测的结果中相同基因的表达向量之间的Pearson相关系数（PCC）来量化每种整合方法的预测性能（方法）。我们首先检查了已知标记基因空间分布的预测结果。例如，Lein等人报道Igsf21和Rprm在皮质的L5/L6层高度表达。与数据集 4（seqFISH+；Smart-seq；小鼠皮层）的 ground truth 相比，Tangram 在预测 Igsf21 的空间分布方面表现最好（PCC = 0.79），gimVI、SpaGE 和 Seurat 紧随其后（PCC=0.77、0.71和0.70）（图 2a）。对于Rprm的空间表达分布，由SpaGE和Seurat生成的结果具有最高的PCC值（PCC=0.79），其次是SpaOTsc、gimVI、Tangram和LIGER（PCC=0.78、0.71、0.66、0.65）（图2b）。

图 2 | 比较能够预测基因空间表达分布的8种整合方法的准确性
图 2- a/b：数据集 4（seqFISH+；Smart-seq；小鼠皮层）中 Igsf21 (a) 和 Rprm (b) 的空间分布，包括ground truth和每种整合方法的预测结果；

我们还检查了数据集 42（ST；10X Chromium；人类鳞状细胞癌）中 COL17A1 基因空间分布的预测结果。 COL17A1 是鳞状细胞癌基底细胞的已知标记基因。 Tangram、gimVI、novoSpaRc和SpaGE成功预测COL17A1在鳞癌基底细胞中高表达；值得注意的是，这四种整合方法的 PCC 值分别为 0.86（Tangram）、0.84（gimVI）、0.76（novoSpaRc）和 0.70（SpaGE），高于其他整合方法的最佳结果（Seurat，0.48； SpaOTsc，0.40；LIGER，0.31；stPlus，0.27）。

为了进一步量化每种整合方法的预测精度，我们采用了除 PCC 之外的三个指标：（1）结构相似度指数（SSIM），它结合了均值、方差和协方差来衡量预测结果与 ground truth 之间的相似性；（2）均方根误差（RMSE），预测分布与ground truth之间的绝对误差； (3) Jensen-Shannon 散度 (JS)，它使用相对信息熵来衡量两个分布之间的差异。对于一个基因，较高的 PCC/SSIM 或较低的 RMSE/JS 值表示更好的预测准确性。我们还通过汇总四个指标（方法）来定义准确度分数（accuracy score ，AS），以简化对每种整合方法的准确度评估（AS 值越高表示性能越好）。

以数据集 4 (seqFISH+; Smart-seq; mouse cortex) 作为基于图像的空间转录组技术示例，Tangram、gimVI 和 SpaGE 明显优于其他整合方法。具体来说，我们发现 Tangram、gimVI 和 SpaGE 的平均 PCC/SSIM 分别为 0.54/0.45、0.52/0.43 和 0.49/0.39，高于其他 5 种方法的 PCC/SSIM 值，以及平均 RMSE/JS 这三种方法中分别为 0.94/0.18、0.97/0.19 和 0.99/0.21，低于其他方法的平均 RMSE/JS（图 2c）。此外，Tangram、gimVI 和 SpaGE 预测的平均 AS 分别为 1.0、0.875 和 0.75，高于其他方法（图 2d）。我们还计算了数据集 42 中所有转录本的预测结果的 PCC、SSIM、RMSE、JS 和 AS 值（作为基于 seq 的空间转录组学方法的示例），并发现基于这些指标， Tangram 和 gimVI 优于其他整合方法（扩展数据图 1b，c）。

为了系统地评估8种整合方法对未检测到基因的空间表达分布的预测准确性，我们确定了所有 45 个配对数据集的预测结果的 PCC、SSIM、RMSE、JS 和 AS 值（图 2e，扩展数据）图 2）。 Tangram、gimVI 和 SpaGE 预测的平均 AS 分别为 0.96、0.84 和 0.69，均超过 Seurat (0.50)、SpaOTsc (0.55)、LIGER (0.25)、novoSpaRc (0.47)、和 stPlus (0.31)。请注意，当我们分别评估基于图像的数据集、基于 seq 的数据集和 32 个模拟数据集时，Tangram仍然是表现最好的整合方法，其次是 gimVI 和 SpaGE（扩展数据图 3a-c）。由于 10X Visium、seqFISH、MERFISH 和 Slide-seq 已发布的数据集超过 3 个，我们在处理使用这4种空间转录组技术获得的数据时，进一步比较了8种整合方法的 AS得分（扩展数据图 3d-g）。我们发现，对于从 10X Visium、seqFISH 和 MERFISH 平台生成的数据，Tangram、gimVI 和 SpaGE优于其他整合方法，并且 Tangram 和 gimVI 是处理 Slide-seq 数据集的顶级方法。

几种整合方法（例如，Seurat、LIGER、SpaGE 和 stPlus）在整合之前默认对空间转录组数据进行标准化。在这里，我们测试四种输入表达矩阵方案：（1）空间数据的原始表达矩阵和scRNA-seq数据的原始表达矩阵（R-R）； (2)空间数据的归一化表达矩阵和scRNA-seq数据的原始表达矩阵(N-R)； (3)空间数据的原始表达矩阵和scRNA-seq数据的归一化表达矩阵(R-N)； (4)空间数据的归一化表达矩阵和scRNA-seq数据的归一化表达矩阵(N-N)。

有趣的是，对于 28 对基于seq的数据集，由 Tangram、gimVI、SpaGE、Seurat、SpaOTsc 和 LIGER 生成的基因空间分布，在使用 R-R 和 R-N 输入方案时，比使用 N-R 或 N-N 输入方案时具有显着更高的 PCC 值，并且在 28 个配对数据集中的 16 个中观察到了这种趋势（P 值 < 0.01，配对 t 检验）（扩展数据图4和5a)；对于 SpaGE、Seurat、SpaOTsc 和 novoSpaRc，在 28 个配对数据集中的 19 个，使用 R-R 输入方案的结果的 PCC 值高于其他输入方案的结果（P < 0.01），并且在 28 个配对数据集中的 18 个中， stPlus 使用 R-R 输入方案时，生成的结果具有更高的 PCC 值比 N-N 输入方案（P < 0.05）。

对于 15 对基于图像的数据集（扩展数据图 4 和 5b），由 Tangram、gimVI、SpaGE 和 Seurat 生成的基因空间分布，使用 R-R 或 R-N 输入方案时，比使用 N-R 或 N-N 输入方案时，具有更高的 PCC 值（15 个数据集中的 11 个，P 值 < 0.05）；使用 R-R 输入方案时，SpaGE、Seurat 和 LIGER 的 PCC 值高于使用其他输入方案时的 PCC 值（15 个数据集中的 11 个，P < 0.05）；使用 R-R 输入方案时 SpaOTsc 的 PCC 值高于使用 NN 输入方案时的 PCC 值（15 个数据集中的 12 个，P < 0.05）。尽管如此，应该强调的是，无论使用什么输入方案，Tangram总是优于其他集成方法（扩展数据图 5c-f）。

3-矩阵稀疏性的影响

值得注意的是，对于数据集 12、13、40 和 44，所有8种整合方法在预测基因的空间表达分布方面的准确性都较低（即平均 PCC/SSIM<0.3，扩展数据图 2）。我们通过计算4个指标（PCC、SSIM、RMSE 和 JS）之间的相关系数来研究这些数据集的整合方法明显较差的性能，并考虑了空间转录组数据集的几个特征，包括表达矩阵（空间转录组和 scRNA-seq 数据的稀疏性定义为表达矩阵中零元素的百分比）、检测到的基因数、检测到的spot数和每个spot的基因数。最终，我们发现 8 种方法的 JS 值均随着表达矩阵稀疏性的增加而线性增加（P 值 < 1 × 10−6 ，决定系数（R2 ）≥0.50）（扩展数据图 6 ）。

为了进一步表征矩阵稀疏性的影响，我们接下来在输入一个非常稀疏的空间表达矩阵（从稀疏度低于 0.7 的高质量数据集下采样）时评估每种整合方法的性能。具体来说，我们检查了从 >100 个spot和捕获 >1,000 个基因的高质量空间转录组数据集。为了模拟具有“高稀疏性”的基因表达矩阵，我们采用 Splatter 和 Scuttle对原始表达矩阵中的非零元素进行不同程度的下采样（方法）。然后，我们使用原始和下采样的表达矩阵作为8种整合方法的输入。

首先，我们评估基因表达矩阵稀疏性对预测已知标记基因空间分布的影响。据Drew等人报道，Cplx1 在皮质的 L5 层高度表达。查看数据集 4 的原始数据和下采样数据（下采样率 = 0.8）中的 Cplx1，我们针对原始数据和下采样数据，观察到由 Tangram、gimVI 和 SpaGE 预测的 Cplx1 基因的空间分布（扩展数据图 7a）。

然后，我们通过计算原始数据和下采样数据，超过PCC 阈值 0.5 的基因比例来评估每种整合方法的性能，我们将其视为“稳健性得分”（RS）。针对数据集 4，Tangram 具有最高的 RS 值（0.60），其次是 gimVI（0.55），然后是 SpaGE（0.51）（图 3a）。此外，我们注意到（1）RS 值随着下采样率的增加而降低；（2）Tangram、gimVI 和 SpaGE 的 RS 值始终高于 Seurat、SpaOTsc、LIGER、novoSpaRc 和 stPlus（图 3b)。包括来自 19 个数据集的下采样数据的组合分析再次突出了 Tangram、gimVI 和 SpaGE 的强大性能：即使下采样率达到 0.8，这三种方法的平均 RS 值仍然 > 0.50（图 3）。总之，Tangram、gimVI 和 SpaGE在预测高度稀疏数据集的基因空间表达分布方面优于其他整合方法。

4-细胞类型去卷积算法的性能

使用 10X Visium 和 ST 等空间转录组技术时，会遇到的一个常见问题是组织切片中的每个spot可能包含多个细胞，因此不可能正确分配每个spot的细胞类型组成。如上所述，Seurat、SpaOTsc、Tangram 和 novoSpaRc 能够将来自scRNA-seq数据的每个细胞类型分配给空间转录组数据的spot，这意味着它们可以用来解卷积每个spot的细胞类型。为解决此问题，研究者开发了Cell2location、SpatialDWLS、RCTD、Stereoscope、DestVI、STRIDE、SPOTlight 和DSTG。

为了比较 12 种整合方法在预测spot的细胞类型组成方面的性能，我们通过“网格化”一个没有此问题的数据集（数据集 10 ，使用 STARmap 获得；Smart-seq；小鼠视觉皮层）来模拟 ST 和 10X Visium 数据集遇到的这种“多细胞spot问题”。已经报告了数据集 10 中每个spot的细胞类型组成，并且可以在模拟每个spot中具有潜在模糊细胞类型分配的数据集时用作基本事实（图 4a 和方法）。原始数据集 10 捕获了 1,549 个细胞，对应15 种细胞类型。网格化后，模拟数据有 189 个spot，每个点包含 1-18 个细胞。我们绘制了 L4 兴奋性神经元的位置，发现 RCTD 和 Stereoscope 在 PCC 值 (0.87) 方面表现更好，其次是 Tangram (0.85)、Cell2loacation (0.83)、STRIDE (0.80)、SPOTlight (0.79)、Seurat ( 0.76）、SpaOTsc（0.74）和 DSTG（0.71）（图 4b）。然后，我们使用 PCC、SSIM、RMSE、JS 和 AS 指标来量化 12 种整合方法在预测网格数据集 10 中spot的细胞类型组成的准确性（图 4c 和扩展数据图 8a）。 RCTD 的 AS 得分最高 (0.94)，其次是 Stereoscope (0.92)。

为了比较 12 种整合方法在预测spot的细胞类型组成方面的性能，我们通过“网格化”一个单细胞分辨率空间转录组数据集，没有spot问题（数据集 10 ，使用 STARmap 获得；Smart-seq；小鼠视觉皮层）来模拟 ST 和 10X Visium 数据集遇到的这种“多细胞spot问题”。已经报告了数据集 10 中每个spot的细胞类型组成，并且可以在模拟每个spot中具有潜在模糊细胞类型分配的数据集时用作基本事实（图 4a 和方法）。原始数据集 10 捕获了 1,549 个细胞，对应15 种细胞类型。网格化后，模拟数据有 189 个spot，每个spot 包含 1-18 个细胞。我们绘制了 L4 兴奋性神经元的位置，发现 RCTD 和 Stereoscope 在 PCC 值 (0.87) 方面表现更好，其次是 Tangram (0.85)、Cell2loacation (0.83)、STRIDE (0.80)、SPOTlight (0.79)、Seurat ( 0.76）、SpaOTsc（0.74）和 DSTG（0.71）（图 4b）。然后，我们使用 PCC、SSIM、RMSE、JS 和 AS 指标来量化 12 种整合方法在预测网格数据集 10 中spot的细胞类型组成的准确性（图 4c 和扩展数据图 8a）。 RCTD 的 AS 得分最高 (0.94)，其次是 Stereoscope (0.92)。

我们还对包含 14 种细胞类型的 524 个细胞的数据集 4（seqFISH+；Smart-seq；小鼠皮层）进行了相同的分析。在“网格化”之后，模拟数据集有 72 个spot（扩展数据图 8b）。使用 L5/6 兴奋性神经元位置的基本事实，我们发现 SpatialDWLS、RCTD、Tangram、Cell2location 和 Stereoscope 的 L5/6 分配的 PCC 值为 0.88、0.86、0.85、0.83 和 0.81 兴奋性神经元，高于其他整合方法（扩展数据图 8c）。此外，在数据集 4 的所有细胞类型的预测结果中，SpatialDWLS、Tangram 和 RCTD 的 AS 值排名前 1、2 和 3 位（分别为 1.0、0.92 和 0.83），其次是 Cell2location（0.67）和 Stereoscope （0.65）（图 4d）。

我们在从 scRNA-seq 数据集（补充表 2 和 3 和方法）合成的 32 个模拟数据集，进一步量化了这些整合方法中spot的细胞类型去卷积中的性能。由于这些 scRNA-seq 数据集中每个细胞的细胞类型信息已被数据源论文报道，因此可以从其包含的细胞中推断出模拟spot的细胞类型组成。请注意，novoSpaRc 和 SpaOTsc 需要每个spot的空间位置信息，因此被排除在外，因为模拟数据集中没有空间位置信息。我们使用 32 个模拟数据集作为基本事实来评估其余 10 种集成方法（包括 Seurat、Tangram、Cell2location、SpatialDWLS、RCTD、Stereoscope、DestVI、STRIDE、SPOTlight 和 DSTG）在对spot中的细胞类型进行去卷积时的性能。我们发现 Cell2location、SpatialDWLS 和 STRIDE 的平均 PCC 和 SSIM 值分别为 0.83/0.75、0.78/0.71 和 0.83/0.69，高于其他整合方法，以及 Cell2location 的平均 RMSE 和 JS 值， SpatialDWLS、RCTD 和 STRIDE 分别为 0.08/0.33、0.10/0.32、0.096/0.37 和 0.11/0.37，低于其他整合方法（扩展数据图 8e）。我们还使用4个指标的聚合分数（即 AS 分数）来对这些整合方法在预测spot的细胞类型组成方面的性能进行排名，我们发现 Cell2location、SpatialDWLS、RCTD 和 STRIDE 优于其他整合方法（图 4e）。

5-计算资源

我们使用所有45个配对数据集来比较8种整合方法消耗的计算资源，这些方法可以预测未检测到的转录本的空间分布（补充表 4）。我们使用相同配置的CPU服务器（2.2GHz, 45 MB L3 cache, 144 CPU cores）来测试每种方法。我们知道gimVI和Tangram可以支持 GPU 处理；然而，当处理最大的数据集 40（空间转录组数据中的 19,522 个spot和 scRNA-seq 数据中的 26,252 个细胞）时，这两种整合方法在我们的GPU服务器（具有 12 GB 内存的 NVIDIA Tesla K80）上报告了内存错误。值得注意的是，Seurat 和 LIGER 处理每个数据集的CPU时间不到 10 分钟，而Tangram和LIGER消耗的内存不到 32 GB。

然后，我们评估了各种数据属性（包括 scRNA-seq 数据中的细胞数、空间数据中的spot数以及用于训练的基因数）对这8种整合方法消耗的计算资源的影响。通过对数据集40中的细胞数和spot数以及数据集6中的训练基因数进行下采样，我们发现在预测未检测到的转录本细胞空间分布的8 种整合方法中，Seurat始终是计算效率最高的方法（扩展数据图 9a-c）。

为了比较对spot细胞类型进行去卷积的 10 种整合方法消耗的计算资源，我们使用了一个包含 10,000 个细胞、20,000 个spot和 56 种细胞类型的大型模拟数据集（详见方法）。对于这个数据集，Cell2locations在我们的GPU平台上报告了内存错误。 Seurat 和 angram占用的 CPU 时间不到 30 分钟，而 Stereoscope、Tangram 和 DestVI 消耗的内存不到 8 GB（扩展数据图 9d）。然后，我们评估了 scRNA-seq 数据中的细胞数、空间数据中的spot数和细胞类型数对这 10 种整合方法所消耗的计算时间的影响，发现 Tangram 和 Seurat 是处理spot细胞类型去卷积的前两种最有效的方法（扩展数据图 9e-g）。

问题：为什么此处没有用12个整合方法，而是10个整合方法？
答：此处是10个整合方法，而不是12个，缺少的是novoSpaRc，SpaOTsc；
原因是32个模拟数据集，SpaOTsc和novoSpaRc被剔除，因为此两个算法都需要spot的空间位置信息，而在我们的模拟数据中是没有此此信息的。

讨论

在这项研究中，我们对整合空间转录组数据和单细胞转录组数据的16种整合方法的性能进行了基准测试。我们发现 Tangram、gimVI 和 SpaGE 在预测基因的空间表达分布方面优于其他整合方法，而Cell2location、SpatialDWLS 和 RCTD在组织切片中spot的细胞类型去卷积方面优于其他整合方法。我们的研究帮助研究人员选择合适的工具并优化数据分析工作流程，能够准确有效地整合空间转录组学数据和scRNA-seq 数据。我们还提供了一个基准分析流程（https://github.com/QuKunLab/SpatialBenchmarking）和一个指导表（下表），总结了所有基准方法的属性和性能，以指导研究人员选择与其数据组合相匹配的合适工具。

基于概率模型结合负二项式或泊松分布构建的方法，例如 gimVI、Cell2location 和 RCTD，通常在预测转录本的空间分布或spot去卷积的细胞类型方面表现更好。深度学习算法也被应用于几种整合方法中，其中Tangram是预测未检测到的转录本空间分布的最佳方法之一。从技术上讲，Tangram在模型中采用了非凸优化，并且在损失函数中只选择了scRNA-seq观测值的最优子集。这些策略的组合可能有助于提高工具的预测能力。

我们的比较分析中一个观察结果是，空间转录组表达矩阵的稀疏性严重影响了8种整合方法在预测RNA转录本空间分布的性能。有多种策略可用于解决空间转录组学表达矩阵的稀疏性问题：1.研究人员可以增加测序深度、2.筛选spot和保留严格质控的基因，以减少质控后表达矩阵的稀疏性，或3.考虑应用缺值填补的算法（例如，SAVER、MAGIC和 WEDGE）来插补基因表达式矩阵中的零数值。

空间转录组学的另一个潜在应用是预测空间上彼此接近的两种细胞类型之间的配体-受体相互作用。为此任务开发了许多分析工具，例如 SpaOTsc、Giotto、CellChat、NicheNet、ICELLNET 和 SingleCellSignalR。 然而，不同方法的结果存在巨大差异，使得信息比较变得困难。例如，只有一小部分 (<5%) 的预测配体-受体相互作用在 > 3 种方法 (https://github.com/QuKunLab/SpatialBenchmarking/tree/main/FigureData) 的结果是共享。因此，未来基准分析可能依赖于更多的细胞间相互作用的实验验证。

不同细胞通讯分析结果比较
比较6种预测细胞间相互作用的方法（SpaOTsc、Giotto、CellChat、ICELLNET、NicheNet和SingleCellSignal）预测的配体-受体相互作用。
a：这些方法在数据集 4 中 L2/3 eNeuron 和 L4 eNeuron 之间共享的配体-受体相互作用的数量；
b-d：通过不同方法对数据集 1 中每个细胞类型对发现的配体-受体相互作用的数量
c：数据集4（ seqFISH；Fluidigm c1；小鼠原肠胚形成）
b：数据集1（seqFISH+；Smart-seq；小鼠皮层）；
d：数据集 41（Slide-seqV2；Drop-seq；小鼠海马区）；
对于所有方法，我们使用空间转录组学数据的表达量作为输入。

问题：为什么不用单细胞细胞通讯的主流软件cellphonedb？
猜测是cellphonedb没有针对小鼠配体-受体数据库；

问题：为什么这些单细胞细胞通讯的分析结果差异很大？
软件是底层数据库+算法的结合；软件自带的配体-受体数据库差异很大，

scRNA-seq 技术（Drop-seq、Smart-seq 和 10X Chromium）和空间测序技术（FISH、osmFISH、seqFISH、MERFISH、STARmap、ISS、EXseq、BaristaSeq、ST、10X Visium、Slide-seq、seq scope和HDST）的广泛多样性和快速发展状态使该基准分析的任务复杂化。我们采用3种不同的视角来克服测序技术多样性带来的挑战：（1）我们将空间转录组数据集分为两类（基于seq的技术和基于图像的技术）； (2) 我们添加了迄今为止已发布三个以上数据集的四种空间转录组学技术（MERFISH、seqFISH、Slide-seq 和 10X Visium）的整合方法的独立比较； (3) 我们使用了几个内在参数（例如，捕获的基因数量、捕获的spot数和表达矩阵的稀疏性）来表征不同空间转录组学技术生成的数据集。还有一些方面（即不同数量的基因和不同的空间组织）可能会影响整合方法的性能和用户的期望。然而，根据目前的数据集，我们发现这些整合方法的性能排名几乎不受生成这些数据的空间转录组技术的影响。

空间转录组技术的进步，例如升级版本的10X Visium 和 BGI Stereo-seq，可以检测直径远小于细胞大小的spot的转录组信息。然而，在接下来的几年内，使用这些技术时，每个spot仍可能不完全对应于单个细胞。此外，考虑到有大量公开可用的空间转录组数据集对不同的研究团体有很高的价值，以及目前正在进行的使用空间转录组技术的密集研究工作，我们认为，整合空间转录组和scRNA-seq数据以确定细胞或未检测到转录本的空间分布将引起人们极大的兴趣。
注：实验技术和算法处于此消彼长的局势，技术的不完善加速算法的研究热度。

总之，本研究对空间转录组数据和scRNA-seq数据的多种整合方法进行独立性比较，以确定细胞类型去卷积或未检测基因的空间分布。我们的结果对于想要分析其空间转录组数据的生物学家或想要改进最先进技术的方法开发人员来说也是很有用的资源。

分析方法

1-数据集的预处理

我们通过以下步骤对每个数据集进行了预处理： (1) 去除低质量细胞。对于 scRNA-seq 数据，我们使用 Seurat 和参数“min.features = 200”来去除捕获少于 200 个基因的细胞。 (2) 表达矩阵的归一化。对于空间转录组学数据集，我们测试了非归一化和归一化的表达矩阵，用于每种整合方法的输入。为了规范化表达式矩阵，我们使用了以下等式：

image.png

其中，Cij 表示spot j 中基因 i 的原始UMI计数，Dij 表示spot j 中基因 i 的标准化UMI计数，N¯ 是每个细胞检测到的转录本的中位数(N¯ is the median number of detected transcripts per cell)。对于 scRNA-seq 数据集，我们使用函数“NormalizeData”和 Seurat 3.2 中的默认参数对其表达矩阵进行归一化。 (3)选择高度可变的基因。对于检测到超过 1,000 个基因的空间转录组数据集，我们使用以下等式计算每个基因的变异系数：
CVi = σi/ui;
其中 CVi 是基因 i 的变异系数； σi 是基因 i 在所有spot的空间分布的标准差s.d. 。ui 是基因 i 在所有spot的平均表达量。我们将 CVi 值最高的 1,000 个基因确定为高可变基因，并使用它们与相应 scRNA-seq 数据中检测到的基因的共享基因来构建每个数据集的ground truth。对于检测到的基因少于 1,000 个的空间转录组数据集，我们使用在空间转录组和 scRNA-seq 数据中检测到的基因来构建每个数据集的ground truth。

2-整合方法的参数设置

我们使用十倍交叉验证评估了8种整合方法的性能，这些方法可以预测未检测到的转录本的空间分布。对于空间数据中的一组基因，我们将基因分成10份，迭代使用9份基因进行整合（即用于训练的参考基因集）；剩下的一部分基因用于预测。每个整合方法的参数如下所述设置：

3-基准指标

我们构建了一个通用分析流程来评估针对45 对数据集不同整合方法的性能。在分析流程中，我们使用以下5个指标来评估每种整合方法。
1) PCC
皮尔逊相关系数(Pearson Correlation Coefficient, pcc)可以度量两个随机变量的线性相关程度，pcc的取值范围为[-1,1]。
PCC 值使用以下方程计算：

其中，∼xi 和 xi 分别是基因 i 在 ground truth 和预测结果中的空间表达值； ui 和 ~ui 分别是基因 i 在 ground truth 和预测结果中的平均表达值； σi 和 ∼σi 分别在ground truth和预测结果中基因i的空间表达值得标准差。对于一个基因，较高的 PCC 值表示更好的预测准确性。

2) SSIM
图像质量评估算法SSIM（结构相似性）
SSIM的全称为structural similarity index，即为结构相似性，是⼀种衡量两幅图像相似度的指标。该指标⾸先由德州⼤学奥斯丁分校的图像和视频⼯程实验室(Laboratory for Image and Video Engineering)提出。⽽如果两幅图像是压缩前和压缩后的图像，那么SSIM算法就可以⽤来评估压缩后的图像质量。
SSIM如何表征相似性：
我们首先对表达矩阵进行了如下缩放（scale），使得每个基因的表达值都在 0 到 1 之间：

其中 xij 表示基因 i 在spotj 中的表达值，M 是spot数目。然后我们使用缩放的基因表达值，利用以下等式计算每个基因的SSIM值：

其中，ui、~ui、σi 和 ~σi 的定义与计算 PCC 值的定义相似（但针对缩放的基因表达值）； C1 和 C2 分别为 0.01 和 0.03； cov(xi, ∼xi) 是ground truth中基因 i 的表达值（即 x′ i ）与预测结果的表达值（即 ∼x′ i ）之间的协方差。对于一个基因，较高的 SSIM 值表示更好的预测准确度。

3) RMSE
均方根误差是预测值与真实值偏差的平方与观测次数n比值的平方根。
我们首先计算了所有spot的每个基因空间表达值的 z-score，然后使用以下等式计算 RMSE：

其中zij和~zij分别是ground truth和预测结果中spot j中基因i的空间表达的z分数。对于一个基因，较低的 RMSE 值表示更好的预测准确性。

4) JS
JS散度（Jensen-Shannon）用来描述两个分布的相似性。
JS 使用相对信息熵（即 Kullback-Leibler 散度）来确定两个分布之间的差异。我们首先计算每个基因的空间分布概率如下：

其中 xij 表示基因 i 在spot j 的表达值，M 是spot数目，Pij 是基因 i 在spot j 的分布概率。然后我们使用以下等式计算每个基因的 JS 值：

其中，Pi 和 ∼Pi 分别是基因 i 在 ground truth 和预测结果中的空间分布概率向量，KL(ai||bi) 是两个概率分布 ai 和 bi 以及 aij 和 bij 之间的 Kullback-Leibler 散度分别是点 j 中基因 i 的预测概率和真实概率。对于一个基因，较低的 JS 值表示更好的预测精度。

5) AS
我们通过聚合 PCC、SSIM、RMSE 和 JS 来定义 AS，以评估每个数据集的整合方法的相对准确性。对于一个数据集，我们计算了每种整合方法预测的所有基因的平均 PCC、SSIM、RMSE 和 JS。然后我们将不同整合算法的PCC和SSIM值按升序排序得到RANKPCC和RANKSSIM；具有最高 PCC/SSIM 值的方法将具有 RANKPCC/SSIM = N，具有最低 PCC/SSIM 值的方法将具有 RANKPCC/SSIM = 1。我们还对不同整合算法的RMSE和JS值进行了降序排序，得到RANKRMSE和RANKJS； RMSE/JS 值最高的方法的 RANKRMSE/JS = 1，RMSE/JS 值最低的方法的 RANKRMSE/JS = N。最后，我们计算了 RANKPCC、RANKSSIM、RANKRMSE 和 RANKJS获取各个整合方式的AS值，如下：

对于一个数据集，具有最高AS值的算法在整合方法评估中具有最好的性能。

4-模拟“多细胞spot”数据集

为了构造每个spot包含已知细胞类型组成的多细胞数据集，我们对单细胞分辨率的空间转录组数据集 dataset4 和 dataset10 进行网格化操作，，模拟每个spot具有潜在模糊细胞类型的数据集。对于 dataset4(seqFISH+; Smart-seq; mouse cortex)，我们定义了一个 500×500 像素 (~51.5μm) 的正方形作为一个spot区域来网格化 seqFISH+ 切片，此步骤参考 SpatialDWLS 文章中引入的粗粒度操作（https://www.srcesq.com/articles/10.1186/s13059-021-02362-7）。我们将一个网格中所有细胞的基因表达值累加来模拟一个可能包含多个细胞的spot，并以该网格的中心代表spot的位置。 dataset4 共模拟出 72 个spot，我们计算了每个spot的细胞类型的百分比作为 ground truth。对于dataset10 (STARmap; Smart-seq; mouse visual cortex)，我们使用 Seurat 对细胞进行聚类，并使用标记基因 79,80 注释每个聚类的细胞类型。我们使用标记基因 Slc17a7 和 Gad1 分别注释兴奋性神经元和抑制性神经元。 L2/3、L4、L5 和 L6 兴奋性神经元（eL2/3、eL4、eL5、eL6）分别由标记基因 Nov、Rorb、Sulf2 和 Pcp4 注释。此外，VIP、SST 和 PV 抑制神经元分别由标记基因 Vip、Sst 和 Pvalb 注释。小胶质细胞、星形胶质细胞、少突胶质细胞、平滑肌和内皮细胞分别由标记基因 Pdgfra、Aqp4、Enpp2、Mgp 和 Bsg 注释。然后我们使用一个 750 像素的窗口来网格化 STARmap切片。我们将一个网格中所有细胞的基因表达值相加来模拟一个可能包含多个细胞的spot，并将网格的中心作为spot的位置。dataset10 的模拟数据共有 189 个spot，我们计算了每个spot中细胞类型的百分比作为 ground truth。

我们还使用 PCC、SSIM、RMSE 和 JS 来评估 Seurat、SpaOTsc、Tangram 和 novoSpaRc 在将细胞分配到组织切片中的空间位置时的准确性。我们首先计算了每个spot中各种类型细胞的比例。然后我们将每个点的细胞类型比例引入方程式3-方程式9，计算PCC、SSIM、RMSE和JS值，量化预测结果与ground truth之间的相似度（PCC/SSIM）或差异（RMSE/JS）。最后，我们使用双边 Mann-Whitney U 检验来计算不同方法之间预测精度差异的统计显着性。

5-使用scRNA-seq数据集模拟空间数据集

对于32个模拟数据集，我们参考RCTD和Stereoscope算法文章介绍的方法，设计了空间转录组模拟数据的生成流程。对于每个模拟spot，我们首先在 5-15 范围内对均匀分布的细胞数 (Nc) 进行采样，并在 2-6 范围内对均匀分布的细胞类型数 (Nt) 进行采样。然后我们假设这些细胞类型在spot中具有相等的分布可能性P=1/Nt，并从scRNA-seq数据的每个细胞类型中随机分配细胞到spot。为了获得每个空间位置的基因表达值，我们将一个spot的所有细胞的基因表达值相加。参考构建 RCTD 中使用的模拟数据集的方法，我们使用 Scuttle (http://bioconductor.org/packages/release/bioc/html/scuttle.html) 将每个spot的计数下采样到 10%原始值。我们可以通过计算细胞类型对应的细胞数量来获得每个spot的细胞类型的百分比。对于用于评估每种整合方法效率的大型模拟数据集，我们使用与上述相同的算法，但将spot数设置为 20000，细胞数设置为 10000，细胞类型数设置为 56。

6-模拟高稀疏性的数据集

我们对45个配对数据集中的19个（即datasets 4、7、18、19、22、23、24、25、27、28、29、30、31、32、33、36、37、38、42 ) 来测试基因表达矩阵的稀疏性对每种整合方法的影响。为了模拟具有高稀疏度的数据集，我们应用Scuttle包（http://www.bioconductor.org/packages/release/bioc/html/scuttle.html）对数据集的空间表达矩阵进行下采样，我们还使用Splatter包对数据集进行下采样（下采样率 = 0.2、0.4、0.6 和 0.8）。我们以不同的速率对每个数据集进行 10 次下采样，以避免随机选择导致的错误。为了量化基因表达矩阵稀疏性对每种整合方法的影响，我们计算了从原始数据和下采样数据预测基因的空间分布PCC 值均大于 0.5 的百分比，这被定义为稳健性得分。

问题：Scuttle包和Splatter包的下采样有什么区别？
Scuttle包的下采样是downsampleMatrix函数，Splatter下采样是kersplatSimCellMeans，针对SingleCellExperiment with counts matrix的数据进行下采样；
https://bioconductor.org/packages/release/bioc/manuals/splatter/man/splatter.pdf
http://www.bioconductor.org/packages/release/bioc/manuals/scuttle/man/scuttle.pdf

小结：
此处构建高稀疏性的数据，是评估算法的稳健性。

7-计算平台

CPU服务器配置为：4个Intel Xeon E78860v4 CPU（2.2GHz，45 MB 三级缓存，共 144 个 CPU 内核）和 1 TB 内存（DDR4 2,400MHz）的计算集群，对16种整合方法进行CPU测试。
GPU服务器配置：Intel Xeon E5-2680v4 CPU（2.4GHz，35 MB 三级缓存，总共 14 个 CPU 内核）、128 GB 内存和 NVIDIA Tesla K80 GPU（12 GB 内存，共 2496 个 CUDA 核心）服务器，对gimVI和 Tangram方法进行GPU 测试。

计算资源测试包括：
1.测试8种能够预测未检测到转录本的空间分布的整合方法，评估不同数据属性（包括 scRNA-seq 数据中的细胞数、空间数据中的spot数以及用于训练的基因数）对消耗的计算资源的影响。
2.对数据集40中的细胞数和spot数进行下采样，并对数据集6中的共享基因数进行下采样。对于可以进行spot细胞类型组成预测的10种整合方法，我们模拟了一个大型数据集（10,000 spots, 20,000 cells）来评估每种方法消耗的计算机资源。
3.我们对模拟数据集中的细胞数、spot数和细胞类型数进行下采样，然后评估这些数据属性对计算资源消耗的影响。

名称解释

基准真相Ground truth：是一个相对概念，是指相对于新的测量方式得到的测量值，作为基准的，由已有的、可靠的测量方式得到的测量值。我们往往会利用基准真相，对新的测量方式进行校准，以降低新测量方式的误差和提高新测量方式的准确性。