此教程介绍了将查询数据集映射到参考数据集的过程。在此示例中，我们映射了 10X Genomics的（2，700 个PBMC）的scRNA-seq 数据集，到我们最近发表的 CITE-seq 参考的 162，000 PBMC，使用了 228 种抗体。我们选择这个例子来演示由参考数据集指导的监督分析如何有助于列举在无监督的分析中难以找到的细胞状态。在第二个例子中，我们演示了如何从不同个体中人类的BMNC数据集进行连续映射到一致的参考集中。
我们以前已经演示了如何使用参考映射方法在查询数据集中注释细胞。在 Seurat v4 中，我们大大提高了整合任务（包括参考映射）的速度和内存要求，还包括将查询细胞投影到以前计算的 UMAP 图上的新功能。
在此教程中，我们演示了如何使用先前已建立的参考集来注释 scRNA-seq 查询集：
- 根据参考集的细胞状态注释每个查询集细胞
- 将每个查询细胞投影到以前计算的 UMAP 图中
- 估算CITE-seq 参考集中测量的表面蛋白质的水平
要运行此教程，请安装Seurat v4。此外，您还需要安装SeuratDisk包。
install.packages("Seurat")
remotes::install_github("mojaveazure/seurat-disk")
library(Seurat)
library(SeuratDisk)
library(ggplot2)
library(patchwork)
示例1：绘制人类外周血细胞图
一个多模式PBMC参考数据集
我们加载参考集（下载这里），并可视化预先计算的UMAP。此参考集以 h5Seurat 文件的形式存储。
reference <- LoadH5Seurat("../data/pbmc_multimodal.h5seurat")
DimPlot(object = reference, reduction = "wnn.umap", group.by = "celltype.l2", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()
映射
为了演示此多模式参考集的映射，我们将使用10x Genomics生成的 2，700 个 PBMC 数据集。
library(SeuratData)
InstallData('pbmc3k')
参考集使用SCTransform()标准化，使用相同的方法标准化查询集。
pbmc3k <- SCTransform(pbmc3k, verbose = FALSE)
然后，我们在参考集和查询集之间找到锚点。我们使用预计算监督 PCA （spca） 转换来展示这一示例。我们建议将受监督的 PCA 用于 CITE-seq 数据集，并在此教程的下一部分演示如何计算此转换。但是，您也可以使用标准的 PCA 转换。
anchors <- FindTransferAnchors(
reference = reference,
query = pbmc3k,
normalization.method = "SCT",
reference.reduction = "spca",
dims = 1:50
)
然后，我们将细胞类型标签和蛋白质数据从参考集转移到查询集。此外，我们将查询数据投影到参考集的 UMAP 结构上。
pbmc3k <- MapQuery(
anchorset = anchors,
query = pbmc3k,
reference = reference,
refdata = list(
celltype.l1 = "celltype.l1",
celltype.l2 = "celltype.l2",
predicted_ADT = "ADT"
),
reference.reduction = "spca",
reduction.model = "wnn.umap"
)
探索映射结果
我们现在可以可视化 2，700 个查询细胞。它们已投影到参考集的 UMAP 图中，每个亚群都收到了两个级别（第 1 级和第 2 级）的注释。
p1 = DimPlot(pbmc3k, reduction = "ref.umap", group.by = "predicted.celltype.l1", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()
p2 = DimPlot(pbmc3k, reduction = "ref.umap", group.by = "predicted.celltype.l2", label = TRUE, label.size = 3 ,repel = TRUE) + NoLegend()
p1 + p2
参考映射数据集帮助我们识别以前混淆在无监督的查询数据集分析中的细胞类型。举比如血浆树突状细胞（pDC）、造血干细胞和祖细胞（HSPC）、调节T细胞（Treg）、CD8幼稚T细胞、CD56+NK细胞、记忆和幼稚的B细胞以及血浆细胞。
每个预测的得分在 0 到 1 之间。
FeaturePlot(pbmc3k, features = c("pDC", "CD16 Mono", "Treg"), reduction = "ref.umap", cols = c("lightgrey", "darkred"), ncol = 3) & theme(plot.title = element_text(size = 10))
我们可以通过探索传统标记基因的表达来验证我们的预测。例如，CLEC4C 和 LIRA4 已被报道为 pDC 的标记，这与我们的预测一致。同样，如果我们执行差异表达来识别 Tregs 的标记，我们会识别一组传统marker，包括 RTKN2、CTLA4、FOXP3 和 IL2RA。
Idents(pbmc3k) <- 'predicted.celltype.l2'
VlnPlot(pbmc3k, features = c("CLEC4C", "LILRA4"), sort = TRUE) + NoLegend()
treg_markers <- FindMarkers(pbmc3k, ident.1 = "Treg", only.pos = TRUE, logfc.threshold = 0.1)
print(head(treg_markers))
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## FOXP3 8.428389e-41 0.2002914 0.152 0.002 1.059617e-36
## RP11-1399P15.1 1.702575e-21 0.4749271 0.273 0.021 2.140478e-17
## RTKN2 3.890305e-17 0.2339929 0.121 0.005 4.890892e-13
## TIGIT 3.773433e-16 0.5934260 0.424 0.065 4.743960e-12
## F5 4.035005e-16 0.2318413 0.121 0.005 5.072808e-12
## FRS2 2.983908e-15 0.1906088 0.152 0.009 3.751369e-11
最后，我们可以可视化根据CITE-seq参考推断出的表面蛋白质的大致水平。
DefaultAssay(pbmc3k) <- 'predicted_ADT'
# see a list of proteins: rownames(pbmc3k)
FeaturePlot(pbmc3k, features = c("CD3-1", "CD45RA", "IgD"), reduction = "ref.umap", cols = c("lightgrey", "darkgreen"), ncol = 3)
计算新的 UMAP 可视化图
在前面的示例中，我们可以在映射到参考组的 UMAP 图后可视化查询组细胞。保持一致的可视化有助于解释新的数据集。但是，如果细胞在查询数据集中存在而不在参考集存在，则它们将投影到参考集中最相似的细胞中。这是 UMAP 包所建立的预设行为和功能，但可能会掩盖查询集中感兴趣的新细胞类型。
在我们的手稿中，我们映射了一个包含发育和分化的中性粒细胞的查询数据集，这个细胞类型不包括在我们的参考集中。我们发现，在合并参考和查询集后计算新的 UMAP（"从头可视化"）有助于识别这些类型，如文中补充图 8 所示。在"从头"可视化中，查询集中特异的细胞状态仍然被分离。在此示例中，2，700 PBMC 不包含特异的细胞状态，但我们演示了如何计算可视化他们。
我们强调，如果用户试图映射基础样本不是 PBMC 的数据集，或包含参考集中不存在的细胞类型，计算一个从头可视化的图是解释其数据集的重要一步。
#merge reference and query
reference$id <- 'reference'
pbmc3k$id <- 'query'
refquery <- merge(reference, pbmc3k)
refquery[["spca"]] <- merge(reference[["spca"]], pbmc3k[["ref.spca"]])
refquery <- RunUMAP(refquery, reduction = 'spca', dims = 1:50)
DimPlot(refquery, group.by = 'id', shuffle = TRUE)
示例2：绘制人体骨髓细胞图
一个多模式骨髓单核细胞BMNC 参考数据集
作为第二个例子，我们绘制了由人类细胞图谱产生的8个单个捐献者骨髓单核细胞（BMNC）的数据集。我们使用加权邻近分析（WNN）中的人类BMNC的CITE-seq作为参考。
此教程显示的参考映射功能与前面的 PBMC 示例相同。此外，我们还演示：
- 如何构建受监督的 PCA （sPCA） 转换
- 如何连续的将多个数据集映射到相同的参考集
- 优化步骤，以进一步提高映射速度
# Both datasets are available through SeuratData
library(SeuratData)
#load reference data
InstallData("bmcite")
bm <- LoadData(ds = "bmcite")
#load query data
InstallData('hcabm40k')
hcabm40k <- LoadData(ds = "hcabm40k")
参考数据集包含一个WNN图，反映了RNA和蛋白质数据在这个CITE-seq实验中的加权组合。
我们可以根据此图进行 UMAP 可视化。设置return.model = TRUE，使我们能够投影查询数据集到此UMAP图上。
bm <- RunUMAP(bm, nn.name = "weighted.nn", reduction.name = "wnn.umap",
reduction.key = "wnnUMAP_", return.model = TRUE)
DimPlot(bm, group.by = "celltype.l2", reduction = "wnn.umap")
计算 一个sPCA 转换
正如我们的手稿所描述的，我们首先计算一个"监督"的PCA。它可识别能封装 WNN 图形结构的转录组数据的转换，允许蛋白质和RNA测量的加权组合到"监督"的PCA，并突出最相关的变异来源。计算此转换后，我们可以将其投影到查询数据集上。我们还可以计算和投影 PCA ，但在使用 WNN 分析构建的多模式参考集时建议使用 sPCA。
sPCA 计算一旦完成，可以快速投影到每个查询数据集上。
bm <- ScaleData(bm, assay = 'RNA')
bm <- RunSPCA(bm, assay = 'RNA', graph = 'wsnn')
计算缓存的邻近索引
由于我们将多个查询样本映射到相同的参考集，我们可以缓存仅涉及参考集的特定步骤。此步骤是可选的，但在映射多个样本时将提高速度。
我们计算参考sPCA空间中的前50个邻近值。我们将此信息存储在参考集的Seurat 对象中的spca.annoy.neighbors，并通过cache.index = TRUE缓存烦人的索引数据结构。
bm <- FindNeighbors(
object = bm,
reduction = "spca",
dims = 1:50,
graph.name = "spca.annoy.neighbors",
k.param = 50,
cache.index = TRUE,
return.neighbor = TRUE,
l2.norm = TRUE
)
查询数据集预处理
在这里，我们将演示映射多个捐赠者骨髓样本查询集到多模式骨髓参考集。这些查询数据集来自人体细胞图谱 （HCA） 中的免疫细胞骨髓数据，可通过 SeuratData 获得。此数据集以单个合并对象的身份提供，有 8 个捐赠者。我们首先将数据拆分为 8 个单独的 Seurat 对象，每个捐赠者可单独映射一个对象。
library(dplyr)
library(SeuratData)
InstallData('hcabm40k')
hcabm40k.batches <- SplitObject(hcabm40k, split.by = "orig.ident")
然后，我们用和参考集相同的方式将查询集标准化。在这里，参考是标准化使用NormalizeData()。如果引用已使用SCTransform()标准化，则查询集也必须使用SCTransform()标准化。
hcabm40k.batches <- lapply(X = hcabm40k.batches, FUN = NormalizeData, verbose = FALSE)
映射
然后，我们在每个捐赠者查询数据集和多模式参考集之间找到锚点。此命令经过优化，通过传递预先计算的参考邻近集，并关闭锚点过滤来最大限度地缩短映射时间。
anchors <- list()
for (i in 1:length(hcabm40k.batches)) {
anchors[[i]] <- FindTransferAnchors(
reference = bm,
query = hcabm40k.batches[[i]],
k.filter = NA,
reference.reduction = "spca",
reference.neighbors = "spca.annoy.neighbors",
dims = 1:50
)
}
然后，我们单独映射每个数据集。
for (i in 1:length(hcabm40k.batches)) {
hcabm40k.batches[[i]] <- MapQuery(
anchorset = anchors[[i]],
query = hcabm40k.batches[[i]],
reference = bm,
refdata = list(
celltype = "celltype.l2",
predicted_ADT = "ADT"),
reference.reduction = "spca",
reduction.model = "wnn.umap"
)
}
探索映射结果
现在映射已完成，我们可以可视化单个对象的结果
p1 <- DimPlot(hcabm40k.batches[[1]], reduction = 'ref.umap', group.by = 'predicted.celltype', label.size = 3)
p2 <- DimPlot(hcabm40k.batches[[2]], reduction = 'ref.umap', group.by = 'predicted.celltype', label.size = 3)
p1 + p2 + plot_layout(guides = "collect")
我们还可以将所有对象合并到一个数据集中。注意，它们都已整合到参考集所定义的公共空间中。然后，我们可以一起可视化结果。
# Merge the batches
hcabm40k <- merge(hcabm40k.batches[[1]], hcabm40k.batches[2:length(hcabm40k.batches)], merge.dr = "ref.umap")
DimPlot(hcabm40k, reduction = "ref.umap", group.by = "predicted.celltype", label = TRUE, repel = TRUE, label.size = 3) + NoLegend()
我们可以可视化查询集细胞中的基因表达、cluster预测分数、（计算）表面蛋白质表达水平：
p3 <- FeaturePlot(hcabm40k, features = c("rna_TRDC", "rna_MPO", "rna_AVP"), reduction = 'ref.umap',
max.cutoff = 3, ncol = 3)
# cell type prediction scores
DefaultAssay(hcabm40k) <- 'prediction.score.celltype'
p4 <- FeaturePlot(hcabm40k, features = c("CD16 Mono", "HSC", "Prog-RBC"), ncol = 3,
cols = c("lightgrey", "darkred"))
# imputed protein levels
DefaultAssay(hcabm40k) <- 'predicted_ADT'
p5 <- FeaturePlot(hcabm40k, features = c("CD45RA", "CD16", "CD161"), reduction = 'ref.umap',
min.cutoff = 'q10', max.cutoff = 'q99', cols = c("lightgrey", "darkgreen") ,
ncol = 3)
p3 / p4 / p5
