引言

本系列开启R中单细胞RNA-seq数据分析教程[1]，持续更新，欢迎关注，转发！

2.3. 使用 LIGER 进行数据整合

除了 Harmony 和 Seurat，Evan Macosko 实验室开发的 LIGAR 也是被基准论文重点介绍的另一个数据整合工具。LIGAR 通过集成非负矩阵分解来识别共享和数据集特定的因素，以进行联合分析。该方法的详细数学原理可以在 论文 中找到。它作为 R 语言中的 liger 包实现，并提供了一个适用于 Seurat 对象的包装器，这也依赖于 R 语言中的额外包 SeuratWrappers。

library(liger)
library(SeuratWrappers)

seurat <- merge(seurat_DS1, seurat_DS2) %>%
    FindVariableFeatures(nfeatures = 3000)
seurat <- ScaleData(seurat, split.by = "orig.ident", do.center = FALSE)
seurat <- RunOptimizeALS(seurat, k = 20, lambda = 5, split.by = "orig.ident")
seurat <- RunQuantileAlignSNF(seurat, split.by = "orig.ident")
seurat <- RunUMAP(seurat, dims = 1:ncol(seurat[["iNMF"]]), reduction = "iNMF")
seurat <- FindNeighbors(seurat, reduction = "iNMF", dims = 1:ncol(Embeddings(seurat, "iNMF"))) %>%
    FindClusters(resolution = 0.6)

# You may also want to save the object
saveRDS(seurat, file="integrated_liger.rds")

和之前一样，接下来将使用 UMAP 来展示整合后的结果，包括数据集、聚类情况以及一些特征图。

plot1 <- UMAPPlot(seurat, group.by="orig.ident")
plot2 <- UMAPPlot(seurat, label = T)
plot3 <- FeaturePlot(seurat, c("FOXG1","EMX1","DLX2","LHX9"), ncol=2, pt.size = 0.1)
((plot1 / plot2) | plot3) + plot_layout(width = c(1,2))

如果你想提升 LIGER 整合的效果，除了在 FindVariableFeatures 函数中设置 nfeatures 参数之外，RunOptimizeALS 函数中的参数也很重要，比如 k 和 lambda。LIGAR 提供了 suggestK 和 suggestLambda 两个函数来帮助设定这两个参数。遗憾的是，这两个参数并没有对应的 Seurat 包装函数，或者需要使用独立的 liger 包及其 LIGER 数据类型来使用这两个函数，而且这两个函数运行起来相当慢。人们也可以根据一些原则进行猜测和调整，比如数据中如果有更多的子结构，可能需要更大的 k 值；更大的 lambda 值会更强地抑制数据集特定的效应，这可能有助于更好地混合不同数据集中的细胞，但也可能以过度整合为代价（例如混合具有不同表达特征的细胞）。

2.4. 使用 MNN 进行数据整合

MNN，由 EMBL-EBI 的 John Marioni 实验室开发，是最早为单细胞 RNA 测序数据整合或批次校正开发的算法之一。它通过两个不同样本/批次之间的细胞的相互最近邻来估计细胞特定的校正向量，以引入对查询细胞降维（例如 PCA）的校正。MNN 还引入了一个排序机制，因此也支持两个以上样本/批次的整合。

RunFastMNN 函数接受一个包含多个 Seurat 对象的列表作为输入，其中每个 Seurat 对象代表一个单独的样本或批次。如果需要根据某个元数据列将一个 Seurat 对象拆分成多个子对象，可以使用 Seurat 包提供的 SplitObject 函数来实现。

library(SeuratWrappers)

seurat_samples <- SplitObject(seurat, "orig.ident")
seurat_mnn <- RunFastMNN(seurat_samples)
seurat[['mnn']] <- CreateDimReducObject(Embeddings(seurat_mnn, "mnn")[colnames(seurat),], key="MNN_")
seurat <- RunUMAP(seurat, dims = 1:20, reduction = "mnn")
seurat <- FindNeighbors(seurat, reduction = "mnn", dims = 1:20) %>%
    FindClusters(resolution = 0.6)

# You may also want to save the object
saveRDS(seurat, file="integrated_mnn.rds")

接下来，可以通过 UMAP 嵌入图来评估整合方法的效果。

plot1 <- UMAPPlot(seurat, group.by="orig.ident")
plot2 <- UMAPPlot(seurat, label = T)
plot3 <- FeaturePlot(seurat, c("FOXG1","EMX1","DLX2","LHX9"), ncol=2, pt.size = 0.1)
((plot1 / plot2) | plot3) + plot_layout(width = c(1,2))

整合的结果看起来很有前景。在大多数情况下，MNN 使用默认参数就能取得很好的效果。当然，用户也可以通过调整一些参数来优化整合，比如改变特征数量或提供一套完全定制的特征集给整合过程。这可以通过在 RunFastMNN 包装函数中设置 features 参数来实现。此外，用户还可以向原始函数传递更多参数（例如，在 batchelor 包中的 fastMNN 函数，可以指定计算的 PCs 数量）。

2.5. 利用 RSS 整合到 BrainSpan 数据

Seurat、Harmony、LIGER 和 MNN 可能是目前最广泛使用的通用单细胞 RNA 测序数据整合方法，但也存在其他方法和概念可以应用于数据整合。其中一个概念是，如果存在一个包含多个样本的参考数据集，并且这些样本间的差异能够反映出样本中细胞类型的异质性，那么通过将每个细胞与这些参考样本的转录组相似性作为表示，而不是直接使用其转录组轮廓，可以有效地过滤掉技术噪声，同时保留关键信息。

进行这种分析时，首先需要一个优质的参考数据集。对于脑类器官样本而言，艾伦脑图谱（Allen Brain Atlas）提供的 BrainSpan 人类大脑批量 RNA-seq 数据集，涵盖了从早期胎儿发育到成人的阶段，是一个非常好的选择。

ref_brainspan <- readRDS("data/ext/brainspan_fetal.rds")

下一步，将计算每个细胞与参考样本之间的相似度，或者进行标准化的皮尔逊相关性分析。simspec 包提供了一个方便的函数来完成这一步骤。分析结果作为一种降维表示，会保存在 Seurat 对象中（默认名称为 "rss"）。之后，可以利用这种降维表示进行 tSNE/UMAP 可视化和聚类分析等。

library(simspec)
seurat <- merge(seurat_DS1, seurat_DS2)
seurat <- ref_sim_spectrum(seurat, ref)
seurat <- RunUMAP(seurat, reduction="rss", dims = 1:ncol(Embeddings(seurat, "rss")))
seurat <- FindNeighbors(seurat, reduction = "rss", dims = 1:ncol(Embeddings(seurat, "rss"))) %>%
    FindClusters(resolution = 0.6)

plot1 <- UMAPPlot(seurat, group.by="orig.ident")
plot2 <- UMAPPlot(seurat, label = T)
plot3 <- FeaturePlot(seurat, c("FOXG1","EMX1","DLX2","LHX9"), ncol=2, pt.size = 0.1)
((plot1 / plot2) | plot3) + plot_layout(width = c(1,2))

观察到了令人满意的细胞轨迹，两个样本的细胞似乎以一种合理的模式混合。然而，你可能已经注意到，在比较聚类结果和特定基因（例如 LHX9）表达时存在问题。

即使你对这一结果感到满意，RCA/RSS 有一个明显的局限：必须有一个优质的参考数据集可供使用，以便能够在不损失太多信息的情况下计算相似性。如果你的数据集中包含一些在参考数据中完全未被覆盖的有趣信号，这些信号很可能会被忽略。由于 RSS 完全基于与参考数据的相似性来表示数据，如果参考数据没有变化，那么改善结果的空间就非常有限。在 ref_sim_spectrum 函数中，唯一可能通过改变而有益的参数是 method 参数，它定义了要计算的相关性类型。默认情况下，使用的是皮尔逊相关（method = "pearson"），但也可以选择使用斯皮尔曼相关（method = "spearman"）。

2.6. 使用 CSS 进行数据整合

最终，将尝试本教程中的最后一种数据整合方法，即 RCA/RSS 的扩展版本。与使用外部参考数据集通过相似性来表示数据中的细胞不同，CSS 首先对每个待整合的单细胞 RNA 测序样本进行细胞聚类，然后使用这些聚类得到的平均表达谱作为参考，来计算相似性。

library(simspec)
seurat <- merge(seurat_DS1, seurat_DS2) %>%
    FindVariableFeatures(nfeatures = 3000) %>%
    ScaleData() %>%
    RunPCA(npcs = 50)
seurat <- cluster_sim_spectrum(seurat, label_tag = "orig.ident", cluster_resolution = 0.3)
seurat <- RunUMAP(seurat, reduction="css", dims = 1:ncol(Embeddings(seurat, "css")))
seurat <- FindNeighbors(seurat, reduction = "css", dims = 1:ncol(Embeddings(seurat, "css"))) %>%
    FindClusters(resolution = 0.6)

plot1 <- UMAPPlot(seurat, group.by="orig.ident")
plot2 <- UMAPPlot(seurat, label = T)
plot3 <- FeaturePlot(seurat, c("FOXG1","EMX1","DLX2","LHX9"), ncol=2, pt.size = 0.1)
((plot1 / plot2) | plot3) + plot_layout(width = c(1,2))

CSS 在合并数据集时基于定义的 PCA 对每个数据集进行聚类，因此在 FindVariableFeatures 函数中的 nfeatures 参数以及 cluster_sim_spectrum 函数中的 dims 参数都会影响所使用的主成分。此外，CSS 对每个数据集独立进行聚类，其中 cluster_sim_spectrum 函数的 cluster_resolution 参数（默认值为 cluster_resolution = 0.6）定义了聚类的分辨率。更高的分辨率能够捕捉数据中更细微的结构，这可能有助于更好地保留数据的结构特征，但同时也可能保留了更多的数据集特有的差异。

[1]Source: https://github.com/quadbio/scRNAseq_analysis_vignette