Seurat软件学习6-多模型参考映射的方法

原创

小胡子刺猬的生信学习123

发布于 2022-10-29 17:47:56

5880

发布于 2022-10-29 17:47:56

文章被收录于专栏：文献分享及代码学习文献分享及代码学习

Seurat软件学习1-多个模型得数据进行整合：/developer/article/2130078

Seurat软件学习2-scrna数据整合分析：/developer/article/2131431

Seurat软件学习3-scrna数据整合分析注释数据集：/developer/article/2133583

Seurat软件学习4-使用RPCA进行快速整合数据集：/developer/article/2134684

Seurat软件学习5-scRNA-Seq和scATAC-Seq数据整合：/developer/article/2136814

Intro: Seurat v4 Reference Mapping

本节介绍了将查询数据集映射到Seurat中带注释的query数据集的过程。在这个例子中，我们将10Xgenomic发表的2700个PBMC的第一个scRNA-seq数据集之一映射到我们最近描述的用228个抗体测量的162,000个PBMC的Cite-seq参考。我们选择这个例子来演示由参考数据集指导的监督分析如何帮助枚举使用非监督分析难以找到的细胞状态。在第二个例子中，我们演示了如何将来自不同个体的人类BMNC的人类细胞图谱数据集序列映射到一致的参考数据集上。

我们以前已经演示了如何使用引用映射方法来注释查询数据集中的单元格标签。在Seurat v4中，我们大幅提高了集成任务的速度和内存存贮，包括引用映射，还包括将查询单元格投影到先前计算的UMAP可视化上的新功能。

在这个小节中，我们演示了如何使用先前建立的参考来解释scRNA-seq查询：

1.基于一组参考定义的细胞状态来注释每个查询细胞

2.将每个拟要做分析的细胞投影到先前计算的UMAP可视化中

3.在CITE-SEQ参考中根据蛋白的预测水平进行计算

要运行此节，请安装cran上提供的seurat v4。此外，您还需要安装SeuratDisk包。

install.packages("Seurat")
remotes::install_github("mojaveazure/seurat-disk")
library(Seurat)
library(SeuratDisk)
library(ggplot2)
library(patchwork)

Example 1: Mapping human peripheral blood cells

多组学的PBMC参考数据集我们从最近的论文中加载参考数据集(参见软件教程原链接)，并可视化预计算的UMAP。此引用存储为h5Seurat文件，这是一种支持在磁盘上存储多模式Seurat对象的格式。

先加载数据集。

reference <- LoadH5Seurat("../data/pbmc_multimodal.h5seurat")
DimPlot(object = reference, reduction = "wnn.umap", group.by = "celltype.l2", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()

Mapping

为了对到这种多model的数据集进行分析，我们将使用由10xgenomic公司的SeuratData数据集，是有2700个PBMC的数据集。

library(SeuratData)
InstallData('pbmc3k')

引用数据中是使用SCTransform()标准化的，因此我们在这里使用相同的方法来标准化数据集。

pbmc3k <- SCTransform(pbmc3k, verbose = FALSE)

然后，我们在引用和query数据集之间找到锚点。正如流程简介中所描述的，我们在这个例子中使用了预计算监督主元分析(SPCA)变换。我们建议对CITE-SEQ数据集使用有监督的主成分分析，并在本教程的的下一个数据集上演示如何计算这种转换。但是，您也可以使用标准的PCA变换。

anchors <- FindTransferAnchors(
  reference = reference,
  query = pbmc3k,
  normalization.method = "SCT",
  reference.reduction = "spca",
  dims = 1:50
)

然后，我们将细胞类型标签和蛋白质数据从引用传输到query数据集上。此外，我们将查询数据映射到引用的数据集的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"
)

What is MapQuery doing?

MapQuery()是三个函数的合成函数：TransferData()、IntegrateEmbedding()和ProjectUMAP()。TransferData()用于传输单元类型标签和分配ADT值。IntegrateEmbedding()和ProjectUMAP()用于将查询数据映射到引用的UMAP结构。使用中间函数执行此操作的等价代码如下：

pbmc3k <- TransferData(
  anchorset = anchors, 
  reference = reference,
  query = pbmc3k,
  refdata = list(
    celltype.l1 = "celltype.l1",
    celltype.l2 = "celltype.l2",
    predicted_ADT = "ADT")
)
pbmc3k <- IntegrateEmbeddings(
  anchorset = anchors,
  reference = reference,
  query = pbmc3k, 
  new.reduction.name = "ref.spca"
)
pbmc3k <- ProjectUMAP(
  query = pbmc3k, 
  query.reduction = "ref.spca", 
  reference = reference, 
  reference.reduction = "spca", 
  reduction.model = "wnn.umap"
)

分析映射数据结果。

现在，我们可以可视化2700个查询的细胞集。它们被投影到引用定义的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

每个观测的值在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的标记，我们就鉴定了到一组经典的标记基因，包括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包建立的预期行为和功能，但可能会掩盖查询中可能感兴趣的新单元类型的存在。

在我们的流程分析中，我们绘制了一个包含发育和分化的中性粒细胞的query数据集，这不包括在我们的参考文献中。我们发现，合并参考和查询后计算新的UMAP(“从头可视化”)有助于识别这些簇，如补充图8所示。在“从头”可视化中，查询中的独特细胞状态保持分离。在这个例子中，2700 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)

Example 2: Mapping human bone marrow cells

多模式BMNC参考数据集作为第二个例子，我们映射了由人类细胞图谱产生的来自八个个体捐赠者的人骨髓单个核细胞(BMNC)的数据集。作为参考，我们使用了我们使用加权最近邻分析(WNN)分析的人类BMNC的Cite-Seq参考

此节展示了与上一个节上的PBMC示例相同的引用两个数据集的功能。此外，我们还分析了：

1.如何进行PCA的转换

2.如何连续映射多个数据集到相同的引用集

3.优化的步骤又提升了相关的运行速度

# 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图，反映了此CITE-SEQ实验中RNA和蛋白质数据的加权组合。

我们可以基于该图计算UMAP可视化。我们设置reur.model=true，这将使我们能够将查询数据集投影到这个可视化上。

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")

Computing an sPCA transformation

正如我们流程中所描述的，我们首先计算一个“有监督的”主成分分析。这确定了最佳封装WNN图结构的转录组数据的转换。这使得蛋白质和RNA测量的加权组合可以‘监督’PCA，并突出最相关的变异来源。在计算此转换之后，我们可以将其投影到查询数据集。我们还可以计算和投影PCA投影，但建议在处理已通过WNN分析构建的多模式参考时使用SPCA。

SPCA计算只执行一次，然后可以快速投影到每个查询数据集上。

bm <- ScaleData(bm, assay = 'RNA')
bm <- RunSPCA(bm, assay = 'RNA', graph = 'wsnn')

Computing a cached neighbor index

因为我们将把多个查询样本映射到同一个引用，所以我们可以缓存只涉及引用的特定步骤。此步骤是可选的，但在映射多个样品时将提高速度。

我们在参考的SPCA空间中计算前50个neighbor。我们将此信息存储在Reference 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
)

How can I save and load a cached annoy index?

如果想要保存和加载一个neighbor对象的缓存索引，可以使用SaveAnnoyIndex()/LoadAnnoyIndex()函数。重要的是，该索引不能正常保存到RDS或RDA文件中，因此它不会在R会话重新启动或为包含它的Seurat对象保存RDS/ReadRDS时正确保存。相反，每次R重新启动或从RDS加载Reference Seurat对象时，使用LoadAnnoyIndex()将Anady索引添加到邻居对象。由SaveAnnoyIndex()创建的文件可以与引用Seurat对象一起分发，并添加到引用中的邻居对象。

bm[["spca.annoy.neighbors"]]
## A Neighbor object containing the 50 nearest neighbors for 30672 cells
SaveAnnoyIndex(object = bm[["spca.annoy.neighbors"]], file = "../data/reftmp.idx")
bm[["spca.annoy.neighbors"]] <- LoadAnnoyIndex(object = bm[["spca.annoy.neighbors"]], file = "../data/reftmp.idx")

Query dataset preprocessing

在这里，我们将演示将多个供体骨髓样本映射到多模式骨髓参考。这些查询数据集来自人类细胞图谱(HCA)免疫细胞图谱骨髓数据集，可通过SeuratData获得。该数据集作为具有8个捐赠者的单个合并对象提供。我们首先将数据分成8个独立的Seurat对象，每个原始捐赠者一个对象单独映射。

library(dplyr)
library(SeuratData)
InstallData('hcabm40k')
hcabm40k.batches <- SplitObject(hcabm40k, split.by = "orig.ident")

然后，我们以与引用相同的方式规范化query。在这里，引用通过NorMalizeData()使用对数标准化进行了标准化。如果引用已使用SCTransform()标准化，则query也必须使用SCTransform()标准化。

hcabm40k.batches <- lapply(X = hcabm40k.batches, FUN = NormalizeData, verbose = FALSE)

Mapping

然后，我们在每个donor查询数据集和多模式引用之间找到锚点。此命令经过优化，通过传入一组预先计算的参考邻居并关闭锚点过滤，最大限度地减少了映射时间。

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"
  )
}

Explore the mapping results

现在已经降维完了，可以可视化一下相关的结果。

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()

我们可以可视化查询细胞中的基因表达、集群预测分数和(归因于)表面蛋白水平：

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