1 Introduction
隨著越來(lái)越多的scRNA-seq數(shù)據(jù)集的產(chǎn)生迈勋，在它們之間進(jìn)行比較就顯得非常關(guān)鍵炬灭。一個(gè)核心的應(yīng)用是比較不同實(shí)驗(yàn)室收集的類似生物來(lái)源的數(shù)據(jù)集，以確保注釋和分析是一致的靡菇。此外重归，隨著大量reference的出現(xiàn)，例如人類細(xì)胞圖譜(HCA)厦凤，一項(xiàng)重要的應(yīng)用將是將來(lái)自新樣本(例如疾病組織)的細(xì)胞投射到上鼻吮，以確定定成分的差異，或檢測(cè)新的細(xì)胞類型较鼓。
scmap是一種將從scRNA-seq實(shí)驗(yàn)中獲得的細(xì)胞投射到不同實(shí)驗(yàn)中鑒定的細(xì)胞類型或細(xì)胞上的方法椎木。
2 SingleCellExperiment class
scmap是建立在Bioconductor的singlecellexperexperiment類之上的。
請(qǐng)閱讀有關(guān)如何從您自己的數(shù)據(jù)創(chuàng)建一個(gè)singlecellexperexperiment的相應(yīng)簡(jiǎn)介博烂。
這里我們將展示一個(gè)關(guān)于如何做到這一點(diǎn)的小示例香椎，但請(qǐng)注意，它不是一個(gè)全面的指南禽篱。
3 scmapinput
如果你已經(jīng)有了一個(gè)SingleCellExperiment對(duì)象, 那么繼續(xù)下一章畜伐。
如果你有一個(gè)包含表達(dá)量數(shù)據(jù)的matrix或dataframe，那么你首先需要?jiǎng)?chuàng)建一個(gè)包含你的數(shù)據(jù)的singlecellexperexperiment對(duì)象躺率。為了便于說(shuō)明玛界，我們將使用scmap提供的示例表達(dá)矩陣万矾。
數(shù)據(jù)集(yan)代表了90個(gè)來(lái)自人類胚胎的細(xì)胞的FPKM基因表達(dá)量。作者(Yan et al.)在原始出版物(ann數(shù)據(jù)框架)中定義了所有細(xì)胞的發(fā)育階段慎框。我們將在后面的投影中使用這些數(shù)據(jù)勤众。
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
注意，細(xì)胞類型信息必須存儲(chǔ)在singlecellexperexperiment對(duì)象的rowData slot的cell_type1列中±鹪啵現(xiàn)在讓我們創(chuàng)建一個(gè)yan數(shù)據(jù)集的singlecellexperexperiment對(duì)象们颜。
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
一旦我們有了一個(gè)singlecellexperexperiment對(duì)象，我們就可以運(yùn)行scmap猎醇。首先窥突，我們需要從輸入數(shù)據(jù)集中選擇最有用的特征(基因)
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...
用紅色突出的gene將用于進(jìn)一步的分析(投影)。這些features存儲(chǔ)在輸入對(duì)象的rowData slot的scmap_features列中硫嘶。默認(rèn)情況下阻问，scmap選擇500個(gè)特征基因(也可以通過(guò)設(shè)置n_features參數(shù)來(lái)控制)
table(rowData(sce)$scmap_features)
## FALSE TRUE
## 19714 500
5 scmap-cluster
5.1 Index
參考數(shù)據(jù)集的scmap-cluster索引是通過(guò)查找每個(gè)cluster的中值基因表達(dá)來(lái)創(chuàng)建的。默認(rèn)情況下沦疾，scmap使用reference中colData槽的cell_type1列來(lái)標(biāo)識(shí)cluster称近。其他列可通過(guò)調(diào)整cluster_col參數(shù)手動(dòng)選擇。
sce <- indexCluster(sce)
indexCluster函數(shù)自動(dòng)寫入reference數(shù)據(jù)集的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
我們也可以將指數(shù)可視化
heatmap(as.matrix(metadata(sce)$scmap_cluster_index))
5.2 Projection
一旦生成了scmap-cluster索引刨秆，我們就可以使用它來(lái)將我們的數(shù)據(jù)集投影到它自己身上(只是為了演示目的)。
這可以一次用一個(gè)索引完成忆畅，但是如果索引以列表的形式提供衡未，scmap還允許同時(shí)投影到多個(gè)索引。
scmapCluster_results <- scmapCluster(
projection = sce,
index_list = list(
yan = metadata(sce)$scmap_cluster_index
)
)
5.3 Results
scmap-cluster將queery數(shù)據(jù)集投影到index_list中定義的所有projections上家凯。細(xì)胞類型標(biāo)簽賦值的結(jié)果合并成一個(gè)矩陣缓醋。
head(scmapCluster_results$scmap_cluster_labs)
## yan
## [1,] "zygote"
## [2,] "zygote"
## [3,] "zygote"
## [4,] "2cell"
## [5,] "2cell"
## [6,] "2cell"
相應(yīng)的相似性存儲(chǔ)在scmap_cluster_siml項(xiàng)中
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數(shù)據(jù)集的組合結(jié)果(選擇對(duì)應(yīng)于跨reference數(shù)據(jù)集的最大相似度的標(biāo)簽)
head(scmapCluster_results$combined_labs)
## [1] "zygote" "zygote" "zygote" "2cell" "2cell" "2cell"
5.4 Visualisation
scmap-cluster的結(jié)果可以可視化為Sankey圖，以顯示細(xì)胞cluster是如何匹配的绊诲。
注意送粱，Sankey圖只有在query和reference數(shù)據(jù)集都被聚類的情況下才會(huì)顯示信息，但沒(méi)有必要為query指定有意義的標(biāo)簽(cluster1掂之、cluster2等就足夠了):
plot(
getSankey(
colData(sce)$cell_type1,
scmapCluster_results$scmap_cluster_labs[,'yan'],
plot_height = 400
)
)
6 scmap-cell
與scmap-cluster相反抗俄，scmap-cell將輸入數(shù)據(jù)集的細(xì)胞投影到reference的單個(gè)細(xì)胞中，而不是投影到細(xì)胞cluster中板惑。
6.1 Stochasticity隨機(jī)性
scmap-cell包含了k-means 步驟橄镜，這使得它具有隨機(jī)性，也就是說(shuō)冯乘，多次運(yùn)行它將會(huì)提供稍微不同的結(jié)果。因此晒夹，我們將固定一個(gè)隨機(jī)seed裆馒，以便用戶能夠準(zhǔn)確地重現(xiàn)我們的結(jié)果:
set.seed(1)
6.2 Index
在scmap-cell索引中姊氓，使用的是乘積量化器算法（product quantizer algorithm），通過(guò)基于特征子集的k-means聚類找到一組亞中心點(diǎn)來(lái)識(shí)別reference中的每個(gè)細(xì)胞喷好。
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索引包含關(guān)于每個(gè)cell的信息，因此不容易顯示梗搅。scmap-cell索引由兩項(xiàng)組成
names(metadata(sce)$scmap_cell_index)
## [1] "subcentroids" "subclusters"
6.2.1 Sub-centroids
subcentroids包含由所選特征禾唁、乘積量化器算法的k和M參數(shù)定義的低維子空間的subcentroids坐標(biāo)(參見?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數(shù)據(jù)集中:
-
yan數(shù)據(jù)集包含 N=90 個(gè)細(xì)胞 - 我們挑選 f=500 個(gè)基因(
scmap默認(rèn)) -
M計(jì)算為 f/10=50(scmap默認(rèn)值為 f≤1000).M是低維子空間的個(gè)數(shù) - 任意低維子空間的特征個(gè)數(shù)等于f/M=10
6.2.2 Sub-clusters
subclusters包含給定細(xì)胞所屬的每一個(gè)亞中心的低維子空間索引
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索引无切，我們就可以使用它們來(lái)投影baron數(shù)據(jù)集荡短。這可以一次用一個(gè)索引完成，但是如果索引以列表的形式提供哆键，scmap允許同時(shí)投影到多個(gè)索引掘托。
scmapCell_results <- scmapCell(
sce,
list(
yan = metadata(sce)$scmap_cell_index
)
)
6.4 Results
scmapCell_results包含列表中每個(gè)引用數(shù)據(jù)集的投影結(jié)果
names(scmapCell_results)
## [1] "yan"
對(duì)于每個(gè)數(shù)據(jù)集，有兩個(gè)矩陣籍嘹。cells matrix包含投影數(shù)據(jù)集的給定細(xì)胞最接近的reference數(shù)據(jù)集細(xì)胞的前10個(gè)(scmap默認(rèn)值)細(xì)胞 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 包含相關(guān)的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
如果細(xì)胞cluster注釋對(duì)于reference數(shù)據(jù)集是可用的闪盔，scmap-cell除了可以找到前10個(gè)最近的鄰居外，還可以使用reference的標(biāo)簽來(lái)注釋投影數(shù)據(jù)集的細(xì)胞辱士。它通過(guò)查看前3個(gè)最近的鄰居(scmap默認(rèn)值)來(lái)做到這一點(diǎn)泪掀，如果它們都屬于reference中的同一個(gè)cluster，并且它們的最大相似度高于閾值(scmap默認(rèn)值為0.5)颂碘，則將一個(gè)投影細(xì)胞分配給相應(yīng)的reference cluster族淮。
scmapCell_clusters <- scmapCell2Cluster(
scmapCell_results,
list(
as.character(colData(sce)$cell_type1)
)
)
scmap-cell的結(jié)果格式與scmap-cluster提供的結(jié)果格式相同(見上文)
head(scmapCell_clusters$scmap_cluster_labs)
## yan
## [1,] "zygote"
## [2,] "zygote"
## [3,] "zygote"
## [4,] "unassigned"
## [5,] "unassigned"
## [6,] "unassigned"
相應(yīng)的相似性存儲(chǔ)在scmap_cluster_siml項(xiàng)中
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
總結(jié)：
總的來(lái)說(shuō)，scmap采用了一種類似于將reads mapping到參考基因組上的策略凭涂。首先需要人為選擇一個(gè)數(shù)據(jù)集作為所謂的reference祝辣，然后利用scmap將其他的數(shù)據(jù)集逐一映射到這個(gè)參考數(shù)據(jù)集上，這樣從橫向平鋪?zhàn)兂煽v向疊加就可以直接避免了因?yàn)闃颖鹃g異質(zhì)性和數(shù)據(jù)集之間的相互干擾而導(dǎo)致的分群混亂切油。
scmap中主要有兩個(gè)方法：一個(gè)叫scmap-cluster蝙斜，另一個(gè)叫scmap-cell。就是將細(xì)胞分到cluster中和將細(xì)胞分到細(xì)胞中的區(qū)別澎胡。
對(duì)于scmap-cluster孕荠，scmap需要選擇一定數(shù)量的gene合集（500），將它們的表達(dá)中值作為該cluster的質(zhì)心（centroid）攻谁，這樣在映射的時(shí)候稚伍，只需要計(jì)算待分群細(xì)胞和每個(gè)cluster的質(zhì)心之間的相似性即可。
scmap-cluster一共使用三中相似性度量方法：皮爾斯相關(guān)系數(shù)（pearson）戚宦，斯皮爾曼相關(guān)系數(shù)（Spearman）和余弦相似性（cosine）个曙。算法默認(rèn)要求這三種相關(guān)性度量中，至少依據(jù)其中的兩種將其分配到群是相同的受楼，否則認(rèn)定該細(xì)胞為“unassigned”垦搬，并且呼寸，這三種方法中得到的最高值也要大于設(shè)定閾值（默認(rèn)為0.7），否則依舊被判定為“unassigned”猴贰。
而scmap-cell則只用了一種相似性計(jì)算方法对雪，即余弦相似性，用來(lái)找細(xì)胞C在reference數(shù)據(jù)集中的最鄰近（the nearest neighbor）米绕。在找尋最近鄰居的過(guò)程中瑟捣，可以改變的參數(shù)是選擇將細(xì)胞C周圍的多少個(gè)相鄰細(xì)胞進(jìn)行計(jì)算。
實(shí)戰(zhàn)
讀取兩個(gè)10x數(shù)據(jù)栅干，構(gòu)建一個(gè)SingleCellExperiment對(duì)象迈套。進(jìn)行分析。由于沒(méi)有colData非驮，沒(méi)有cellType信息交汤，這里使用seurat_clusters信息代替celltype,后續(xù)可以用自動(dòng)注釋工具注釋cluster。
如下圖所示劫笙，左邊樣本query有8個(gè)cluster芙扎，右邊樣本reference有7個(gè)cluster，那么到底左邊的cluster中的細(xì)胞是否跟右邊的cluster細(xì)胞相關(guān)性如何呢填大？
我們可以看到：左邊樣本的cluster0,1戒洼，以及34567都跟右邊樣本cluster6有非常高的相似性，在85%左右允华。這也是非常期待的結(jié)果圈浇。
完整代碼見公眾號(hào)：生信診斷所
