单细胞两个数据集之间的比较：scmap
1 Introduction
随着越来越多的scRNA-seq数据集的产生，在它们之间进行比较就显得非常关键。一个核心的应用是比较不同实验室收集的类似生物来源的数据集，以确保注释和分析是一致的。此外，随着大量reference的出现，例如人类细胞图谱(HCA)，一项重要的应用将是将来自新样本(例如疾病组织)的细胞投射到上，以确定定成分的差异，或检测新的细胞类型。
scmap是一种将从scRNA-seq实验中获得的细胞投射到不同实验中鉴定的细胞类型或细胞上的方法。
2 SingleCellExperiment class
scmap是建立在Bioconductor的singlecellexperexperiment类之上的。
请阅读有关如何从您自己的数据创建一个singlecellexperexperiment的相应简介。
这里我们将展示一个关于如何做到这一点的小示例，但请注意，它不是一个全面的指南。
3 scmapinput
如果你已经有了一个SingleCellExperiment对象, 那么继续下一章。
如果你有一个包含表达量数据的matrix或dataframe，那么你首先需要创建一个包含你的数据的singlecellexperexperiment对象。为了便于说明，我们将使用scmap提供的示例表达矩阵。
数据集(yan)代表了90个来自人类胚胎的细胞的FPKM基因表达量。作者(Yan et al.)在原始出版物(ann数据框架)中定义了所有细胞的发育阶段。我们将在后面的投影中使用这些数据。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager") BiocManager::install("scmap")
library(SingleCellExperiment)
library(scmap)
head(ann)
## cell_type1
## Oocyte..1.RPKM. zygote
## Oocyte..2.RPKM. zygote
## Oocyte..3.RPKM. zygote
## Zygote..1.RPKM. zygote
## Zygote..2.RPKM. zygote
## Zygote..3.RPKM. zygote
yan[1:3, 1:3]
## Oocyte..1.RPKM. Oocyte..2.RPKM. Oocyte..3.RPKM.
## C9orf152 0.0 0.0 0.0
## RPS11 1219.9 1021.1 931.6
## ELMO2 7.0 12.2 9.3
注意，细胞类型信息必须存储在singlecellexperexperiment对象的rowData slot的cell_type1列中。现在让我们创建一个yan数据集的singlecellexperexperiment对象。
sce <- SingleCellExperiment(assays = list(normcounts = as.matrix(yan)),
colData = ann)
logcounts(sce) <- log2(normcounts(sce) + 1)
# use gene names as feature symbols
rowData(sce)$feature_symbol <- rownames(sce)
# remove features with duplicated names
sce <- sce[!duplicated(rownames(sce)), ]
sce
## class: SingleCellExperiment
## dim: 20214 90 ## metadata(0):
## assays(2): normcounts logcounts
## rownames(20214): C9orf152 RPS11 ... CTSC AQP7
## rowData names(1): feature_symbol
## colnames(90): Oocyte..1.RPKM. Oocyte..2.RPKM. ...
## Late.blastocyst..3..Cell.7.RPKM. Late.blastocyst..3..Cell.8.RPKM.
## colData names(1): cell_type1
## reducedDimNames(0):
## altExpNames(0):
4 Feature selection
一旦我们有了一个singlecellexperexperiment对象，我们就可以运行scmap。首先，我们需要从输入数据集中选择最有用的特征(基因)
sce <- selectFeatures(sce, suppress_plot = FALSE)
## Warning in linearModel(object, n_features):
## Your object does not contain counts() slot. Dropouts were calculated using logcounts() slot...
image.png
用红色突出的gene将用于进一步的分析(投影)。这些features存储在输入对象的rowData slot的scmap_features列中。默认情况下，scmap选择500个特征基因(也可以通过设置n_features参数来控制)
table(rowData(sce)$scmap_features)
## FALSE TRUE
## 19714 500
5 scmap-cluster
5.1 Index
参考数据集的scmap-cluster索引是通过查找每个cluster的中值基因表达来创建的。默认情况下，scmap使用reference中colData槽的cell_type1列来标识cluster。其他列可通过调整cluster_col参数手动选择。
sce <- indexCluster(sce)
indexCluster函数自动写入reference数据集的metadata slot的scmap_cluster_index。
head(metadata(sce)$scmap_cluster_index)
## zygote 2cell 4cell 8cell 16cell blast
## ABCB4 5.788589 6.2258580 5.935134 0.6667119 0.000000 0.000000
## ABCC6P1 7.863625 7.7303559 8.322769 7.4303689 4.759867 0.000000
## ABT1 0.320773 0.1315172 0.000000 5.9787977 6.100671 4.627798
## ACCSL 7.922318 8.4274290 9.662611 4.5869260 1.768026 0.000000
## ACOT11 0.000000 0.0000000 0.000000 6.4677243 7.147798 4.057444
## ACOT9 4.877394 4.2196038 5.446969 4.0685468 3.827819 0.000000
我们也可以将指数可视化
heatmap(as.matrix(metadata(sce)$scmap_cluster_index))
image.png
5.2 Projection
一旦生成了scmap-cluster索引，我们就可以使用它来将我们的数据集投影到它自己身上(只是为了演示目的)。
这可以一次用一个索引完成，但是如果索引以列表的形式提供，scmap还允许同时投影到多个索引。
scmapCluster_results <- scmapCluster(
projection = sce,
index_list = list(
yan = metadata(sce)$scmap_cluster_index
)
)
5.3 Results
scmap-cluster将queery数据集投影到index_list中定义的所有projections上。细胞类型标签赋值的结果合并成一个矩阵。
head(scmapCluster_results$scmap_cluster_labs)
## yan
## [1,] "zygote"
## [2,] "zygote"
## [3,] "zygote"
## [4,] "2cell"
## [5,] "2cell"
## [6,] "2cell"
相应的相似性存储在scmap_cluster_siml项中
head(scmapCluster_results$scmap_cluster_siml)
## yan
## [1,] 0.9947609
## [2,] 0.9951257
## [3,] 0.9955916
## [4,] 0.9934012
## [5,] 0.9953694
## [6,] 0.9871041
scmap还提供了所有reference数据集的组合结果(选择对应于跨reference数据集的最大相似度的标签)
head(scmapCluster_results$combined_labs)
## [1] "zygote" "zygote" "zygote" "2cell" "2cell" "2cell"
5.4 Visualisation
scmap-cluster的结果可以可视化为Sankey图，以显示细胞cluster是如何匹配的。
注意，Sankey图只有在query和reference数据集都被聚类的情况下才会显示信息，但没有必要为query指定有意义的标签(cluster1、cluster2等就足够了):
plot(
getSankey(
colData(sce)$cell_type1,
scmapCluster_results$scmap_cluster_labs[,'yan'],
plot_height = 400
)
)
image.png
6 scmap-cell
与scmap-cluster相反，scmap-cell将输入数据集的细胞投影到reference的单个细胞中，而不是投影到细胞cluster中。
6.1 Stochasticity随机性
scmap-cell包含了k-means 步骤，这使得它具有随机性，也就是说，多次运行它将会提供稍微不同的结果。因此，我们将固定一个随机seed，以便用户能够准确地重现我们的结果:
set.seed(1)
6.2 Index
在scmap-cell索引中，使用的是乘积量化器算法（product quantizer algorithm），通过基于特征子集的k-means聚类找到一组亚中心点来识别reference中的每个细胞。
sce <- indexCell(sce)
## Parameter M was not provided, will use M = n_features / 10 (if n_features <= 1000), where n_features is the number of selected features, and M = 100 otherwise.
## Parameter k was not provided, will use k = sqrt(number_of_cells)
与scmap-cluster索引不同，scmap-cell索引包含关于每个cell的信息，因此不容易显示。scmap-cell索引由两项组成
names(metadata(sce)$scmap_cell_index)
## [1] "subcentroids" "subclusters"
6.2.1 Sub-centroids
subcentroids包含由所选特征、乘积量化器算法的k和M参数定义的低维子空间的subcentroids坐标(参见?indexCell)。
length(metadata(sce)$scmap_cell_index$subcentroids)
## [1] 50
dim(metadata(sce)$scmap_cell_index$subcentroids[[1]])
## [1] 10 9
metadata(sce)$scmap_cell_index$subcentroids[[1]][,1:5]
## 1 2 3 4 5
## ZAR1L 0.072987697 0.2848353 0.33713297 0.26694708 0.3051086
## SERPINF1 0.179135680 0.3784345 0.35886481 0.39453521 0.4326297
## GRB2 0.439712934 0.4246024 0.23308320 0.43238208 0.3247221
## GSTP1 0.801498298 0.1464230 0.14880665 0.19900079 0.0000000
## ABCC6P1 0.005544482 0.4358565 0.46276591 0.40280401 0.3989602
## ARGFX 0.341212258 0.4284664 0.07629512 0.47961460 0.1296112
## DCT 0.004323311 0.1943568 0.32117489 0.21259776 0.3836451
## C15orf60 0.006681366 0.1862540 0.28346531 0.01123282 0.1096438
## SVOPL 0.003004345 0.1548237 0.33551596 0.12691677 0.2525819
## NLRP9 0.101524942 0.3223963 0.40624639 0.30465156 0.4640308
在yan数据集中:
-
yan数据集包含 N=90 个细胞 - 我们挑选 f=500 个基因(
scmap默认) -
M计算为 f/10=50(scmap默认值为 f≤1000).M是低维子空间的个数 - 任意低维子空间的特征个数等于f/M=10
6.2.2 Sub-clusters
subclusters包含给定细胞所属的每一个亚中心的低维子空间索引
dim(metadata(sce)$scmap_cell_index$subclusters)
## [1] 50 90
metadata(sce)$scmap_cell_index$subclusters[1:5,1:5]
## Oocyte..1.RPKM. Oocyte..2.RPKM. Oocyte..3.RPKM. Zygote..1.RPKM.
## [1,] 6 6 6 6
## [2,] 5 5 5 5
## [3,] 5 5 5 5
## [4,] 3 3 3 3
## [5,] 6 6 6 6
## Zygote..2.RPKM.
## [1,] 6
## [2,] 5
## [3,] 5
## [4,] 3
## [5,] 6
6.3 Projection
一旦生成了scmap-cell索引，我们就可以使用它们来投影baron数据集。这可以一次用一个索引完成，但是如果索引以列表的形式提供，scmap允许同时投影到多个索引。
scmapCell_results <- scmapCell(
sce,
list(
yan = metadata(sce)$scmap_cell_index
)
)
6.4 Results
scmapCell_results包含列表中每个引用数据集的投影结果
names(scmapCell_results)
## [1] "yan"
对于每个数据集，有两个矩阵。cells matrix包含投影数据集的给定细胞最接近的reference数据集细胞的前10个(scmap默认值)细胞 id
scmapCell_results$yan$cells[,1:3]
## Oocyte..1.RPKM. Oocyte..2.RPKM. Oocyte..3.RPKM.
## [1,] 1 1 1
## [2,] 2 2 2
## [3,] 3 3 3
## [4,] 11 11 11
## [5,] 5 5 5
## [6,] 6 6 6
## [7,] 7 7 7
## [8,] 12 8 12
## [9,] 9 9 9
## [10,] 10 10 10
similaritiesmatrix 包含相关的cosine 相似度
scmapCell_results$yan$similarities[,1:3]
## Oocyte..1.RPKM. Oocyte..2.RPKM. Oocyte..3.RPKM. ## [1,] 0.9742737 0.9736593 0.9748542 ## [2,] 0.9742274 0.9737083 0.9748995 ## [3,] 0.9742274 0.9737083 0.9748995 ## [4,] 0.9693955 0.9684169 0.9697731 ## [5,] 0.9698173 0.9688538 0.9701976 ## [6,] 0.9695394 0.9685904 0.9699759 ## [7,] 0.9694336 0.9686058 0.9699198 ## [8,] 0.9694091 0.9684312 0.9697699 ## [9,] 0.9692544 0.9684312 0.9697358 ## [10,] 0.9694336 0.9686058 0.9699198
6.5 Cluster annotation
如果细胞cluster注释对于reference数据集是可用的，scmap-cell除了可以找到前10个最近的邻居外，还可以使用reference的标签来注释投影数据集的细胞。它通过查看前3个最近的邻居(scmap默认值)来做到这一点，如果它们都属于reference中的同一个cluster，并且它们的最大相似度高于阈值(scmap默认值为0.5)，则将一个投影细胞分配给相应的reference cluster。
scmapCell_clusters <- scmapCell2Cluster(
scmapCell_results,
list(
as.character(colData(sce)$cell_type1)
)
)
scmap-cell的结果格式与scmap-cluster提供的结果格式相同(见上文)
head(scmapCell_clusters$scmap_cluster_labs)
## yan
## [1,] "zygote"
## [2,] "zygote"
## [3,] "zygote"
## [4,] "unassigned"
## [5,] "unassigned"
## [6,] "unassigned"
相应的相似性存储在scmap_cluster_siml项中
head(scmapCell_clusters$scmap_cluster_siml)
## yan
## [1,] 0.9742737
## [2,] 0.9737083
## [3,] 0.9748995
## [4,] NA
## [5,] NA
## [6,] NA
head(scmapCell_clusters$combined_labs)
## [1] "zygote" "zygote" "zygote" "unassigned" "unassigned" ## [6] "unassigned"
6.6 Visualisation
plot(
getSankey(
colData(sce)$cell_type1,
scmapCell_clusters$scmap_cluster_labs[,"yan"],
plot_height = 400
)
)
参考：http://bioconductor.org/packages/release/bioc/vignettes/scmap/inst/doc/scmap.html
总结：
总的来说，scmap采用了一种类似于将reads mapping到参考基因组上的策略。首先需要人为选择一个数据集作为所谓的reference，然后利用scmap将其他的数据集逐一映射到这个参考数据集上，这样从横向平铺变成纵向叠加就可以直接避免了因为样本间异质性和数据集之间的相互干扰而导致的分群混乱。
scmap中主要有两个方法：一个叫scmap-cluster，另一个叫scmap-cell。就是将细胞分到cluster中和将细胞分到细胞中的区别。
对于scmap-cluster，scmap需要选择一定数量的gene合集（500），将它们的表达中值作为该cluster的质心（centroid），这样在映射的时候，只需要计算待分群细胞和每个cluster的质心之间的相似性即可。
scmap-cluster一共使用三中相似性度量方法：皮尔斯相关系数（pearson），斯皮尔曼相关系数（Spearman）和余弦相似性（cosine）。算法默认要求这三种相关性度量中，至少依据其中的两种将其分配到群是相同的，否则认定该细胞为“unassigned”，并且，这三种方法中得到的最高值也要大于设定阈值（默认为0.7），否则依旧被判定为“unassigned”。
而scmap-cell则只用了一种相似性计算方法，即余弦相似性，用来找细胞C在reference数据集中的最邻近（the nearest neighbor）。在找寻最近邻居的过程中，可以改变的参数是选择将细胞C周围的多少个相邻细胞进行计算。
实战
读取两个10x数据，构建一个SingleCellExperiment对象。进行分析。由于没有colData，没有cellType信息，这里使用seurat_clusters信息代替celltype,后续可以用自动注释工具注释cluster。
如下图所示，左边样本query有8个cluster，右边样本reference有7个cluster，那么到底左边的cluster中的细胞是否跟右边的cluster细胞相关性如何呢？
plot_zoom_png (3).png
plot_zoom_png.png
plot_zoom_png (1).png
plot_zoom_png (2).png
我们可以看到：左边样本的cluster0,1，以及34567都跟右边样本cluster6有非常高的相似性，在85%左右。这也是非常期待的结果。
2022-01-25_162804.png
完整代码见公众号：生信诊断所
