scmap package vignette
2017-12-11
Contents
1 Introduction
As more and more scRNA-seq datasets become available, carrying out comparisons between them is key. A central application is to compare datasets of similar biological origin collected by different labs to ensure that the annotation and the analysis is consistent. Moreover, as very large references, e.g. the Human Cell Atlas (HCA), become available, an important application will be to project cells from a new sample (e.g. from a disease tissue) onto the reference to characterize differences in composition, or to detect new cell-types.
scmap is a method for projecting cells from a scRNA-seq experiment on to the cell-types or cells identified in a different experiment. A copy of the scmap manuscript is available on bioRxiv.
2 SingleCellExperiment class
scmap is built on top of the Bioconductor’s SingleCellExperiment class. Please read corresponding vignettes on how to create a SingleCellExperiment from your own data. Here we will show a small example on how to do that but note that it is not a comprehensive guide.
3 scmap input
If you already have a SingleCellExperiment object, then proceed to the next chapter.
If you have a matrix or a data frame containing expression data then you first need to create an SingleCellExperiment object containing your data. For illustrative purposes we will use an example expression matrix provided with scmap. The dataset (yan) represents FPKM gene expression of 90 cells derived from human embryo. The authors (Yan et al.) have defined developmental stages of all cells in the original publication (ann data frame). We will use these stages in projection later.
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. zygoteyan[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.3Note that the cell type information has to be stored in the cell_type1 column of the rowData slot of the SingleCellExperiment object.
Now let’s create a SingleCellExperiment object of the yan dataset:
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)
isSpike(sce, "ERCC") <- grepl("^ERCC-", 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):
## spikeNames(1): ERCC4 Feature selection
Once we have a SingleCellExperiment object we can run scmap. Firstly, we need to select the most informative features (genes) from our input dataset:
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...
Features highlighted with the red colour will be used in the futher analysis (projection).
Features are stored in the scmap_features column of the rowData slot of the input object. By default scmap selects \(500\) features (it can also be controlled by setting n_features parameter):
table(rowData(sce)$scmap_features)##
## FALSE TRUE
## 19714 5005 scmap-cluster
5.1 Index
The scmap-cluster index of a reference dataset is created by finding the median gene expression for each cluster. By default scmap uses the cell_type1 column of the colData slot in the reference to identify clusters. Other columns can be manually selected by adjusting cluster_col parameter:
sce <- indexCluster(sce)The function indexCluster automatically writes the scmap_cluster_index item of the metadata slot of the reference dataset.
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.000000One can also visualise the index:
heatmap(as.matrix(metadata(sce)$scmap_cluster_index))
5.2 Projection
Once the scmap-cluster index has been generated we can use it to project our dataset to itself (just for illustrative purposes). This can be done with one index at a time, but scmap also allows for simultaneous projection to multiple indexes if they are provided as a list:
scmapCluster_results <- scmapCluster(
projection = sce,
index_list = list(
yan = metadata(sce)$scmap_cluster_index
)
)5.3 Results
scmap-cluster projects the query dataset to all projections defined in the index_list. The results of cell label assignements are merged into one matrix:
head(scmapCluster_results$scmap_cluster_labs)## yan
## [1,] "zygote"
## [2,] "zygote"
## [3,] "zygote"
## [4,] "2cell"
## [5,] "2cell"
## [6,] "2cell"Corresponding similarities are stored in the scmap_cluster_siml item:
head(scmapCluster_results$scmap_cluster_siml)## yan
## [1,] 0.9947609
## [2,] 0.9951257
## [3,] 0.9955916
## [4,] 0.9934012
## [5,] 0.9953694
## [6,] 0.9871041scmap also provides combined results of all reference dataset (choose labels corresponding to the largest similarity across reference datasets):
head(scmapCluster_results$combined_labs)## [1] "zygote" "zygote" "zygote" "2cell" "2cell" "2cell"5.4 Visualisation
The results of scmap-cluster can be visualized as a Sankey diagram to show how cell-clusters are matched (getSankey() function). Note that the Sankey diagram will only be informative if both the query and the reference datasets have been clustered, but it is not necessary to have meaningful labels assigned to the query (cluster1, cluster2 etc. is sufficient):
plot(
getSankey(
colData(sce)$cell_type1,
scmapCluster_results$scmap_cluster_labs[,'yan'],
plot_height = 400
)
)6 scmap-cell
In contrast to scmap-cluster, scmap-cell projects cells of the input dataset to the individual cells of the reference and not to the cell clusters.
6.1 Stochasticity
scmap-cell contains k-means step which makes it stochastic, i.e. running it multiple times will provide slightly different results. Therefore, we will fix a random seed, so that a user will be able to exactly reproduce our results:
set.seed(1)6.2 Index
In the scmap-cell index is created by a product quantiser algorithm in a way that every cell in the reference is identified with a set of sub-centroids found via k-means clustering based on a subset of the features.
sce <- indexCell(sce)Unlike scmap-cluster index scmap-cell index contains information about each cell and therefore can not be easily visualised. scmap-cell index consists of two items:
names(metadata(sce)$scmap_cell_index)## [1] "subcentroids" "subclusters"6.2.1 Sub-centroids
subcentroids contains coordinates of subcentroids of low dimensional subspaces defined by selected features, k and M parameters of the product quantiser algorithm (see ?indexCell).
length(metadata(sce)$scmap_cell_index$subcentroids)## [1] 50dim(metadata(sce)$scmap_cell_index$subcentroids[[1]])## [1] 10 9metadata(sce)$scmap_cell_index$subcentroids[[1]][,1:5]## 1 2 3 4 5
## ZAR1L 0.4552243410 0.27559246 0.06636228 0.004572068 0.412884318
## SERPINF1 0.2874948804 0.38978317 0.40006850 0.012769272 0.337788854
## GRB2 0.0466392569 0.43123962 0.48846120 0.512813918 0.034120317
## GSTP1 0.0901068638 0.17828401 0.04338333 0.856549383 0.047546101
## ABCC6P1 0.3888941269 0.41834667 0.40325323 0.013910631 0.412169693
## ARGFX 0.0005263088 0.46082606 0.59345935 0.041671318 0.008650784
## DCT 0.3467325173 0.20609549 0.04589694 0.005001592 0.349504834
## C15orf60 0.3561004784 0.08340984 0.00000000 0.000000000 0.348892484
## SVOPL 0.3036048672 0.13902245 0.03672849 0.011347071 0.353105729
## NLRP9 0.4575775369 0.31341775 0.27713192 0.032822330 0.416507454In the case of our yan dataset:
yandataset contains \(N = 90\) cells- We selected \(f = 500\) features (
scmapdefault) Mwas calculated as \(f / 10 = 50\) (scmapdefault for \(f \le 1000\)).Mis the number of low dimensional subspaces- Number of features in any low dimensional subspace equals to \(f / M = 10\)
kwas calculated as \(k = \sqrt{N} \approx 9\) (scmapdefault).
6.2.2 Sub-clusters
subclusters contains for every low dimensial subspace indexies of subcentroids which a given cell belongs to:
dim(metadata(sce)$scmap_cell_index$subclusters)## [1] 50 90metadata(sce)$scmap_cell_index$subclusters[1:5,1:5]## Oocyte..1.RPKM. Oocyte..2.RPKM. Oocyte..3.RPKM. Zygote..1.RPKM.
## [1,] 1 1 1 1
## [2,] 1 1 1 1
## [3,] 6 6 6 6
## [4,] 7 7 7 7
## [5,] 2 2 2 2
## Zygote..2.RPKM.
## [1,] 1
## [2,] 1
## [3,] 6
## [4,] 7
## [5,] 26.3 Projection
Once the scmap-cell indexes have been generated we can use them to project the baron dataset. This can be done with one index at a time, but scmap allows for simultaneous projection to multiple indexes if they are provided as a list:
scmapCell_results <- scmapCell(
sce,
list(
yan = metadata(sce)$scmap_cell_index
)
)6.4 Results
scmapCell_results contains results of projection for each reference dataset in a list:
names(scmapCell_results)## [1] "yan"For each dataset there are two matricies. cells matrix contains the top 10 (scmap default) cell IDs of the cells of the reference dataset that a given cell of the projection dataset is closest to:
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,] 4 4 4
## [5,] 5 5 5
## [6,] 6 6 6
## [7,] 7 7 7
## [8,] 8 8 12
## [9,] 9 9 9
## [10,] 10 10 10similarities matrix contains corresponding cosine similarities:
scmapCell_results$yan$similarities[,1:3]## Oocyte..1.RPKM. Oocyte..2.RPKM. Oocyte..3.RPKM.
## [1,] 0.9737049 0.9725332 0.9741216
## [2,] 0.9737049 0.9725332 0.9741216
## [3,] 0.9737049 0.9725332 0.9741216
## [4,] 0.9709274 0.9700484 0.9714928
## [5,] 0.9712800 0.9703344 0.9717317
## [6,] 0.9709274 0.9700484 0.9714928
## [7,] 0.9709521 0.9699978 0.9713831
## [8,] 0.9707755 0.9698334 0.9713324
## [9,] 0.9708824 0.9699931 0.9713373
## [10,] 0.9710590 0.9701575 0.97154036.5 Cluster annotation
If cell cluster annotation is available for the reference datasets, in addition to finding top 10 nearest neighbours scmap-cell also allows to annotate cells of the projection dataset using labels of the reference. It does so by looking at the top 3 nearest neighbours (scmap default) and if they all belong to the same cluster in the reference and their maximum similarity is higher than a threshold (\(0.5\) is the scmap default) a projection cell is assigned to a corresponding reference cluster:
scmapCell_clusters <- scmapCell2Cluster(
scmapCell_results,
list(
as.character(colData(sce)$cell_type1)
)
)scmap-cell results are in the same format as the ones provided by scmap-cluster (see above):
head(scmapCell_clusters$scmap_cluster_labs)## yan
## [1,] "zygote"
## [2,] "zygote"
## [3,] "zygote"
## [4,] "unassigned"
## [5,] "unassigned"
## [6,] "unassigned"Corresponding similarities are stored in the scmap_cluster_siml item:
head(scmapCell_clusters$scmap_cluster_siml)## yan
## [1,] 0.9737049
## [2,] 0.9725332
## [3,] 0.9741216
## [4,] NA
## [5,] NA
## [6,] NAhead(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
)
)7 sessionInfo()
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.6-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] bindrcpp_0.2 scmap_1.1.5
## [3] SingleCellExperiment_1.0.0 SummarizedExperiment_1.8.0
## [5] DelayedArray_0.4.1 matrixStats_0.52.2
## [7] Biobase_2.38.0 GenomicRanges_1.30.0
## [9] GenomeInfoDb_1.14.0 IRanges_2.12.0
## [11] S4Vectors_0.16.0 BiocGenerics_0.24.0
## [13] googleVis_0.6.2 BiocStyle_2.6.1
##
## loaded via a namespace (and not attached):
## [1] reshape2_1.4.3 lattice_0.20-35 colorspace_1.3-2
## [4] htmltools_0.3.6 yaml_2.1.15 rlang_0.1.4
## [7] e1071_1.6-8 glue_1.2.0 GenomeInfoDbData_0.99.1
## [10] bindr_0.1 plyr_1.8.4 stringr_1.2.0
## [13] zlibbioc_1.24.0 munsell_0.4.3 gtable_0.2.0
## [16] codetools_0.2-15 evaluate_0.10.1 labeling_0.3
## [19] knitr_1.17 class_7.3-14 Rcpp_0.12.14
## [22] backports_1.1.1 scales_0.5.0 jsonlite_1.5
## [25] XVector_0.18.0 ggplot2_2.2.1 digest_0.6.12
## [28] stringi_1.1.6 bookdown_0.5 dplyr_0.7.4
## [31] grid_3.4.3 rprojroot_1.2 tools_3.4.3
## [34] bitops_1.0-6 magrittr_1.5 RCurl_1.95-4.8
## [37] lazyeval_0.2.1 proxy_0.4-19 tibble_1.3.4
## [40] randomForest_4.6-12 pkgconfig_2.0.1 Matrix_1.2-12
## [43] assertthat_0.2.0 rmarkdown_1.8 R6_2.2.2
## [46] compiler_3.4.3