---
title: "`scmap` package vignette"
author: "Vladimir Kiselev and Martin Hemberg"
date: "`r Sys.Date()`"
output:
    BiocStyle::html_document:
        toc: true
vignette: >
  %\VignetteIndexEntry{`scmap` package vignette}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r knitr-options, echo=FALSE, message=FALSE, warning=FALSE, cache=TRUE}
library(googleVis)
op <- options(gvis.plot.tag='chart')
```

# 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](http://doi.org/10.1101/150292).

# `SingleCellExperiment` class

`scmap` is built on top of the Bioconductor’s [SingleCellExperiment](https://bioconductor.org/packages/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. 

# `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.](http://dx.doi.org/10.1038/nsmb.2660)) have defined developmental stages of all cells in the original publication (`ann` data frame). We will use these stages in projection later.

```{r , warning=FALSE, message=FALSE}
library(SingleCellExperiment)
library(scmap)

head(ann)
yan[1:3, 1:3]
```

Note 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:
```{r}
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
```

# 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:
```{r}
sce <- selectFeatures(sce, suppress_plot = FALSE)
```

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):
```{r}
table(rowData(sce)$scmap_features)
```

# scmap-cluster

## 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:
```{r}
sce <- indexCluster(sce)
```

The function `indexCluster` automatically writes the `scmap_cluster_index` item of the metadata slot of the reference dataset.
```{r}
head(metadata(sce)$scmap_cluster_index)
```

One can also visualise the index:
```{r , fig.height=7}
heatmap(as.matrix(metadata(sce)$scmap_cluster_index))
```

## 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:
```{r}
scmapCluster_results <- scmapCluster(
  projection = sce, 
  index_list = list(
    yan = metadata(sce)$scmap_cluster_index
  )
)
```

## 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:
```{r}
head(scmapCluster_results$scmap_cluster_labs)
```

Corresponding similarities are stored in the scmap_cluster_siml item:
```{r}
head(scmapCluster_results$scmap_cluster_siml)
```

`scmap` also provides combined results of all reference dataset (choose labels corresponding to the largest similarity across reference datasets):
```{r}
head(scmapCluster_results$combined_labs)
```

## 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):

```{r results='asis', tidy=FALSE, eval=FALSE}
plot(
  getSankey(
    colData(sce)$cell_type1, 
    scmapCluster_results$scmap_cluster_labs[,'yan'],
    plot_height = 400
  )
)
```

# 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.

## 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:
```{r}
set.seed(1)
```


## 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.

```{r message=FALSE, warning=FALSE}
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:

```{r}
names(metadata(sce)$scmap_cell_index)
```

### 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`). 

```{r}
length(metadata(sce)$scmap_cell_index$subcentroids)
dim(metadata(sce)$scmap_cell_index$subcentroids[[1]])
metadata(sce)$scmap_cell_index$subcentroids[[1]][,1:5]
```

In the case of our `yan` dataset:

* `yan` dataset contains $N = 90$ cells
* We selected $f = 500$ features (`scmap` default)
* `M` was calculated as $f / 10 = 50$ (`scmap` default for $f \le 1000$). `M` is the number of low dimensional subspaces
* Number of features in any low dimensional subspace equals to $f / M = 10$
* `k` was calculated as $k = \sqrt{N} \approx 9$ (`scmap` default). 

### Sub-clusters

`subclusters` contains for every low dimensial subspace indexies of `subcentroids` which a given cell belongs to:

```{r}
dim(metadata(sce)$scmap_cell_index$subclusters)
metadata(sce)$scmap_cell_index$subclusters[1:5,1:5]
```

## 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:

```{r}
scmapCell_results <- scmapCell(
  sce, 
  list(
    yan = metadata(sce)$scmap_cell_index
  )
)
```

## Results

`scmapCell_results` contains results of projection for each reference dataset in a list:
```{r}
names(scmapCell_results)
```

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:

```{r}
scmapCell_results$yan$cells[,1:3]
```

`similarities` matrix contains corresponding cosine similarities:

```{r}
scmapCell_results$yan$similarities[,1:3]
```

## 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:

```{r}
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):

```{r}
head(scmapCell_clusters$scmap_cluster_labs)
```

Corresponding similarities are stored in the `scmap_cluster_siml` item:
```{r}
head(scmapCell_clusters$scmap_cluster_siml)
```

```{r}
head(scmapCell_clusters$combined_labs)
```

## Visualisation

```{r results='asis', tidy=FALSE, eval=FALSE}
plot(
  getSankey(
    colData(sce)$cell_type1, 
    scmapCell_clusters$scmap_cluster_labs[,"yan"],
    plot_height = 400
  )
)
```

# sessionInfo()

```{r echo=FALSE}
sessionInfo()
```