Seurat 4.0学习笔记(测试数据与参考数据集的比对)

Multimodal 参考数据的比对

Intro: Seurat v4 Reference Mapping
这篇文章介绍了将查询数据集比对到Seurat中并且将其注释的过程。在此示例中,我们将10X Genomics of 2,700 PBMC发布的第一个scRNA-seq数据集映射到我们最近描述的用228种抗体测量的162,000 PBMC的CITE-seq参考序列。我们选择了此示例来说明参考数据集指导下的监督分析如何帮助枚举在无监督分析中很难找到的细胞状态。在第二个示例中,我们演示了如何将从不同个体剖析的人BMNC的人细胞图集数据连续映射到一致的参考上。

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

在此文章中,我们演示了如何使用先前建立的参考数据来解释scRNA-seq查询:

  • 根据一组参考定义的单元格状态注释每个查询单元格
  • 将每个查询单元格投影到先前计算的UMAP可视化文件中
  • 估算在CITE-seq参考序列中测量的表面蛋白的水平

请安装CRAN上可用的Seurat v4。此外,您将需要安装该SeuratDisk软件包。

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

实施例1:比对人外周血细胞

一个Multimodal PBMC参考数据集

我们从最近的 preprint中加载参考序列(在此处下载),并可视化预先计算的UMAP。此参考存储为h5Seurat文件,该格式允许在磁盘上存储多模式Seurat对象(有关h5Seurat的更多详细信息,请参见此处)。SeuratDisk

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

在这里插入图片描述

映射
为了演示到该multimodal参考的映射,我们将使用由10x Genomics生成并可通过访问的2,700个PBMC的数据集SeuratData。

library(SeuratData)
InstallData('pbmc3k')

参考序列是使用规范化的SCTransform()处理的,因此我们在此处使用相同的方法来规范化处理查询数据。

pbmc3k <- SCTransform(pbmc3k, verbose = FALSE)

然后,我们在引用和查询之间找到锚点。如手稿中所述,在此示例中,我们使用了预先计算的 *supervised PCA(spca)*转换。我们建议对CITE-seq数据集使用 supervised PCA,并在此文章的后续演示如何计算此转换。但是,也可以使用标准的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"
)

什么是MapQuery?
MapQuery()封装了三个函数:TransferData()、IntegrateEmbeddings()和ProjectUMAP()。TransferData()用于传递细胞类型标签,并输入ADT值。IntegrateEmbeddings()和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

在这里插入图片描述

参考映射数据集可帮助我们识别以前在查询数据集的无监督分析中混合的单元格类型。仅举几个例子,包括浆细胞样树突状细胞(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的标记物,我们将鉴定出一组标准标记物,包括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:映射人骨髓细胞

一个Multimodal BMNC参考数据集
作为第二个例子,我们绘制了由人类细胞图集产生的来自八个个体供体的人类骨髓单核(BMNC)细胞的数据集。作为参考,我们使用人类BMNC的CITE-seq参考,我们使用加权最近邻分析(WNN)对其进行了分析。

此小插图具有与上一个选项卡上的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图,反映了此CITE-seq实验中RNA和蛋白质数据的加权组合。

我们可以基于该图计算UMAP可视化。我们设置return.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") 

在这里插入图片描述

计算sPCA转换
如我们的手稿中所述,我们首先计算一个“监督”的PCA。这确定了转录组数据的转换,该转换最好地封装了WNN图的结构。这样可以对蛋白质和RNA的测量结果进行加权组合,以“监督” PCA,并突出显示最相关的变异来源。计算完此转换后,我们可以将其投影到查询数据集上。我们还可以计算和投影PCA投影,但是在处理已通过WNN分析构建的多峰参考时,建议使用sPCA。

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

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

计算 a cached neighbor index
由于我们将多个查询样本映射到同一参考,因此我们可以缓存仅涉及参考的特定步骤。此步骤是可选的,但在映射多个样本时将提高速度。

我们计算参考的sPCA空间中的前50个邻居。我们将此信息存储spca.annoy.neighbors在参考Seurat对象内的对象中,并且还缓存了烦人的索引数据结构(通过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()

在这里插入图片描述

我们可以可视化查询单元格中的基因表达,簇预测分数和(估算的)表面蛋白水平:

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

在这里插入图片描述

  • 2
    点赞
  • 14
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值