Libraries

We need of course destiny, scran for preprocessing, and some tidyverse niceties.

library(conflicted)
library(destiny)
suppressPackageStartupMessages(library(scran))
library(purrr)
library(ggplot2)
library(SingleCellExperiment)

Data

Let’s use data from the scRNAseq[1] package. If necessary, install it via BiocManager::install('scRNAseq').

[1] Risso D, Cole M (2019). scRNAseq: A Collection of Public Single-Cell RNA-Seq Datasets.

# The parts of the help we’re interested in
help('scRNAseq-package', package = 'scRNAseq') %>% repr::repr_html() %>%
    stringr::str_extract_all(stringr::regex('<p>The dataset.*?</p>', dotall = TRUE)) %>% unlist() %>%
    paste(collapse = '') %>% knitr::raw_html()

The dataset fluidigm contains 65 cells from Pollen et al. (2014), each sequenced at high and low coverage.

The dataset th2 contains 96 T helper cells from Mahata et al. (2014).

The dataset allen contains 379 cells from the mouse visual cortex. This is a subset of the data published in Tasic et al. (2016).

379 cells seems sufficient to see something!

allen <- scRNAseq::ReprocessedAllenData()

## see ?scRNAseq and browseVignettes('scRNAseq') for documentation

## loading from cache

## see ?scRNAseq and browseVignettes('scRNAseq') for documentation

## loading from cache

## see ?scRNAseq and browseVignettes('scRNAseq') for documentation

## loading from cache

## see ?scRNAseq and browseVignettes('scRNAseq') for documentation

## loading from cache

## see ?scRNAseq and browseVignettes('scRNAseq') for documentation

## loading from cache

## see ?scRNAseq and browseVignettes('scRNAseq') for documentation

## loading from cache

Preprocessing

We’ll mostly stick to the scran vignette here. Let’s add basic information to the data and choose what to work with.

As scran expects the raw counts in the counts assay, we rename the more accurate RSEM counts to counts. Our data has ERCC spike-ins in an altExp slot:

rowData(allen)$Symbol <- rownames(allen)
rowData(allen)$EntrezID <- AnnotationDbi::mapIds(org.Mm.eg.db::org.Mm.eg.db, rownames(allen), 'ENTREZID', 'ALIAS')

## 

## 'select()' returned 1:many mapping between keys and columns

rowData(allen)$Uniprot <- AnnotationDbi::mapIds(org.Mm.eg.db::org.Mm.eg.db, rownames(allen), 'UNIPROT', 'ALIAS', multiVals = 'list')

## 'select()' returned 1:many mapping between keys and columns

assayNames(allen)[assayNames(allen) == 'rsem_counts'] <- 'counts'
assayNames(altExp(allen, 'ERCC'))[assayNames(altExp(allen, 'ERCC')) == 'rsem_counts'] <- 'counts'
allen

## class: SingleCellExperiment 
## dim: 20816 379 
## metadata(2): SuppInfo which_qc
## assays(4): tophat_counts cufflinks_fpkm counts rsem_tpm
## rownames(20816): 0610007P14Rik 0610009B22Rik ... Zzef1 Zzz3
## rowData names(3): Symbol EntrezID Uniprot
## colnames(379): SRR2140028 SRR2140022 ... SRR2139341 SRR2139336
## colData names(22): NREADS NALIGNED ... Animal.ID passes_qc_checks_s
## reducedDimNames(0):
## mainExpName: endogenous
## altExpNames(1): ERCC

Now we can use it to renormalize the data. We normalize the counts using the spike-in size factors and logarithmize them into logcounts.

allen <- computeSpikeFactors(allen, 'ERCC')
allen <- logNormCounts(allen)
allen

## class: SingleCellExperiment 
## dim: 20816 379 
## metadata(2): SuppInfo which_qc
## assays(5): tophat_counts cufflinks_fpkm counts rsem_tpm logcounts
## rownames(20816): 0610007P14Rik 0610009B22Rik ... Zzef1 Zzz3
## rowData names(3): Symbol EntrezID Uniprot
## colnames(379): SRR2140028 SRR2140022 ... SRR2139341 SRR2139336
## colData names(23): NREADS NALIGNED ... passes_qc_checks_s sizeFactor
## reducedDimNames(0):
## mainExpName: endogenous
## altExpNames(1): ERCC

We also use the spike-ins to detect highly variable genes more accurately:

decomp <- modelGeneVarWithSpikes(allen, 'ERCC')
rowData(allen)$hvg_order <- order(decomp$bio, decreasing = TRUE)

We create a subset of the data containing only rasonably highly variable genes:

allen_hvg <- subset(allen, hvg_order <= 5000L)

Let’s create a Diffusion map. For rapid results, people often create a PCA first, which can be stored in your SingleCellExperiment before creating the Diffusion map or simply created implicitly using DiffusionMap(..., n_pcs = <number>).

However, even with many more principal components than necessary to get a nicely resolved Diffusion Map, the close spatial correspondence between diffusion components and genes are lost.

#reducedDim(allen_hvg, 'pca') <- irlba::prcomp_irlba(t(assay(allen, 'logcounts')), 50)$x

The chosen distance metric has big implications on your results, you should try at least cosine and rankcor.

set.seed(1)
dms <- c('euclidean', 'cosine', 'rankcor') %>% #, 'l2'
    set_names() %>%
    map(~ DiffusionMap(allen_hvg, distance = ., knn_params = list(method = 'covertree')))

## Warning in DiffusionMap(allen_hvg, distance = ., knn_params = list(method =
## "covertree")): You have 5000 genes. Consider passing e.g. n_pcs = 50 to speed
## up computation.

## Warning in DiffusionMap(allen_hvg, distance = ., knn_params = list(method =
## "covertree")): You have 5000 genes. Consider passing e.g. n_pcs = 50 to speed
## up computation.

## Warning in DiffusionMap(allen_hvg, distance = ., knn_params = list(method =
## "covertree")): You have 5000 genes. Consider passing e.g. n_pcs = 50 to speed
## up computation.

TODO: wide plot

dms %>%
    imap(function(dm, dist) plot(dm, 1:2, col_by = 'driver_1_s') + ggtitle(dist)) %>%
    cowplot::plot_grid(plotlist = ., nrow = 1)

grs <- map(dms, gene_relevance)

TODO: wide plot

gms <- imap(grs, function(gr, dist) plot(gr, iter_smooth = 0) + ggtitle(dist))
cowplot::plot_grid(plotlist = gms, nrow = 1)

As you can see, despite the quite different embedding, the rankcor and Cosine diffusion Maps display a number of the same driving genes.

gms[-1] %>% map(~ .$ids[1:10]) %>% purrr::reduce(intersect) %>% cat(sep = ' ')

httr::GET('https://www.uniprot.org/uniprot/', query = list(
    columns = 'id,genes,comment(TISSUE SPECIFICITY)',
    format = 'tab',
    query = rowData(allen)$Uniprot[gms$cosine$ids[1:6]] %>% unlist() %>% paste(collapse = ' or ')
)) %>% httr::content(type = 'text/tab-separated-values', encoding = 'utf-8', )

## Rows: 0 Columns: 1
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): {"timestamp":"2023-10-24T21:04:11.123+00:00","status":500,"error":"...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## # A tibble: 0 × 1
## # ℹ 1 variable:
## #   {"timestamp":"2023-10-24T21:04:11.123+00:00","status":500,"error":"Internal Server Error","message":"","path":"/uniprotkb/columns=id,genes,comment(TISSUE%20SPECIFICITY)&format=tab&query=Q8K2G4%20or%20Q922E6%20or%20Q8QZR0%20or%20Q505Q3%20or%20Q99KR8%20or%20Q923G2%20or%20P62254%20or%20Q5F239"} <chr>