# detecting relevant genes with destiny 3

# Single Cell RNA-Sequencing data and gene relevance

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

    ## Warning: replacing previous import 'utils::findMatches' by
    ## 'S4Vectors::findMatches' when loading 'AnnotationDbi'

    ## 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-04-26T01:44:35.512+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-04-26T01:44:35.512+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>