Seurat 4.0 || 单细胞BMNC多模态参考数据集
男，
一个长大了才会遇到的帅哥，
稳健，潇洒，大方，靠谱。
一段生信缘，一棵技能树。
生信技能树核心成员，单细胞天地特约撰稿人，简书创作者，单细胞数据科学家。
前情回顾
Seurat 4.0 || 您的单细胞数据分析工具箱上新啦
Seurat 4.0 || 单细胞多模态数据整合算法WNN
Seurat 4.0 || 分析scRNA和膜蛋白数据
Seurat 4.0 || WNN整合scRNA和scATAC数据
Seurat 4.0 || 单细胞PBMC多模态参考数据集
正文
多模态数据越来越多地用来分析单细胞的状态，在之前的文章中我们介绍了PBMC的多模态数据集，它使我们可以快速方便地定义PBMC细胞类型。这里我们绘制了一个人类骨髓单核细胞(BMNC)数据集，这些细胞来自8个捐赠者，由人类细胞图谱（HCA）制作。我们使用人类BMNC的CITE-seq参考数据集，并使用加权最近邻分析(WNN)进行分析。
这里展示了与的PBMC示例相同的参考数据映射功能。此外，我们还将演示:
如何构造一个监督的PCA (sPCA)转换 如何映射多个数据集到同一个参考数据集上 优化步骤，进一步提高映射速度
library(Seurat)
# Both datasets are available through SeuratData
library(SeuratData)
#load reference data
# InstallData("bmcite")
bm <- LoadData(ds = "bmcite")
bm
An object of class Seurat
17034 features across 30672 samples within 2 assays
Active assay: RNA (17009 features, 2000 variable features)
1 other assay present: ADT
1 dimensional reduction calculated: spca
#load query data
InstallData('hcabm40k')
library(hcabm40k.SeuratData)
hcabm40k
# hcabm40k <- LoadData(ds = "hcabm40k")
hcabm40k
An object of class Seurat
17369 features across 40000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
head(hcabm40k@meta.data)
orig.ident nCount_RNA nFeature_RNA
MantonBM1_HiSeq_8-CCCAATCGTATGCTTG-1 MantonBM1 3235 856
MantonBM1_HiSeq_8-CTGCTGTAGGACACCA-1 MantonBM1 1566 653
MantonBM1_HiSeq_3-CAGCTGGGTACATGTC-1 MantonBM1 2103 716
MantonBM1_HiSeq_7-GCTTCCAAGATGTAAC-1 MantonBM1 1798 750
MantonBM1_HiSeq_6-AGTGAGGCAGGGTTAG-1 MantonBM1 11540 2760
MantonBM1_HiSeq_5-ATGTGTGCACCGAATT-1 MantonBM1 25616 1966
参考数据集包含一个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)
head(bm@meta.data)
orig.ident nCount_RNA nFeature_RNA nCount_ADT nFeature_ADT lane donor celltype.l1 celltype.l2 RNA.weight
a_AAACCTGAGCTTATCG-1 bmcite 7546 2136 1350 25 HumanHTO4 batch1 Progenitor cells Prog_RBC 0.4827011
a_AAACCTGAGGTGGGTT-1 bmcite 1029 437 2970 25 HumanHTO1 batch1 T cell gdT 0.2417890
a_AAACCTGAGTACATGA-1 bmcite 1111 429 2474 23 HumanHTO5 batch1 T cell CD4 Naive 0.5077136
a_AAACCTGCAAACCTAC-1 bmcite 2741 851 4799 25 HumanHTO3 batch1 T cell CD4 Memory 0.4313079
a_AAACCTGCAAGGTGTG-1 bmcite 2099 843 5434 25 HumanHTO2 batch1 Mono/DC CD14 Mono 0.5685085
a_AAACCTGCACGGTAGA-1 bmcite 2291 783 4658 25 HumanHTO6 batch1 B cell Naive B 0.4255879
DimPlot(bm, group.by = "celltype.l2", reduction = "wnn.umap")
正如我们在手稿中所描述的，我们首先计算一个“监督的”PCA。这标识了WNN图结构的transcriprome转换。使得蛋白质和RNA测量的加权组合能够“监督”PCA，并突出最相关的变异来源。在计算此转换之后，我们可以将其投影到查询数据集中。我们也可以计算和投射一个PCA投影，但是建议在处理由WNN分析构建的多模态引用时使用sPCA。
sPCA计算执行一次，然后可以快速地投影到每个查询数据集中。
bm <- ScaleData(bm, assay = 'RNA')
bm <- RunSPCA(bm, assay = 'RNA', graph = 'wsnn')
因为我们将把多个查询数据映射到同一个参考数据集上，所以我们可以缓存只涉及该参考数据的特定步骤。这个步骤是可选的，但是可以提高映射多个样本的速度。我们计算参考数据集在sPCA空间中的前50个邻居，并将这些信息存储在spca中（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
)
如果你想保存和加载一个用method = "annoy"生成的Neighbor对象的缓存索引，可以使用SaveAnnoyIndex()/LoadAnnoyIndex()函数。重要的是，这个索引通常不能保存到RDS或RDA文件中，因此对于包含它的Seurat对象，它在R会话重启或saveRDS/readRDS期间将不能正确持久。相反，每次R重启时，使用LoadAnnoyIndex()将Annoy 索引添加到邻居对象，或者从RDS加载引用Seurat对象。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")
在这里，我们将演示如何将多个供体骨髓样本映射到多模式骨髓参考数据上。这些查询数据集来自人类细胞图谱(HCA)免疫细胞图谱骨髓数据集，并可通过SeuratData获得。这个数据集是提供8位捐赠者的合并对象。我们首先将数据分割回8个单独的Seurat对象，每个原始捐赠者分别映射。
hcabm40k.batches <- SplitObject(hcabm40k, split.by = "orig.ident")
$MantonBM1
An object of class Seurat
17369 features across 5000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
$MantonBM2
An object of class Seurat
17369 features across 5000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
$MantonBM3
An object of class Seurat
17369 features across 5000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
$MantonBM4
An object of class Seurat
17369 features across 5000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
$MantonBM5
An object of class Seurat
17369 features across 5000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
$MantonBM6
An object of class Seurat
17369 features across 5000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
$MantonBM7
An object of class Seurat
17369 features across 5000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
$MantonBM8
An object of class Seurat
17369 features across 5000 samples within 1 assay
Active assay: RNA (17369 features, 0 variable features)
然后，我们以与参考数据相同的方式均一化查询。在这里，通过NormalizeData()使用log均一化对参考数据进行均一化。如果参考数据已经使用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()
head(hcabm40k@meta.data)
orig.ident nCount_RNA nFeature_RNA predicted.celltype.score predicted.celltype
MantonBM1_HiSeq_8-CCCAATCGTATGCTTG-1 MantonBM1 3235 856 0.9230445 Memory B
MantonBM1_HiSeq_8-CTGCTGTAGGACACCA-1 MantonBM1 1566 653 0.7258886 NK
MantonBM1_HiSeq_3-CAGCTGGGTACATGTC-1 MantonBM1 2103 716 0.8573307 CD4 Memory
MantonBM1_HiSeq_7-GCTTCCAAGATGTAAC-1 MantonBM1 1798 750 1.0000000 CD14 Mono
MantonBM1_HiSeq_6-AGTGAGGCAGGGTTAG-1 MantonBM1 11540 2760 0.8110870 Prog_Mk
MantonBM1_HiSeq_5-ATGTGTGCACCGAATT-1 MantonBM1 25616 1966 1.0000000 Plasmablast
我们可以可视化查询细胞中的基因表达、聚类预测分数和(估算)表面蛋白表达水平:
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
References
[1] 手稿: https://satijalab.org/seurat/v4.0/www.satijalab.org/v4preprint
[2] https://satijalab.org/seurat/v4.0/reference_mapping.html