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.

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 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!

## snapshotDate(): 2022-04-19
## 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:

## 
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 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.

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

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

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.

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

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

TODO: wide plot

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

## Rows: 8 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (3): Entry, Gene names, Tissue specificity
## 
## ℹ 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: 8 × 3
##   Entry  `Gene names` `Tissue specificity`                                      
##   <chr>  <chr>        <chr>                                                     
## 1 Q8QZR0 Tram1l1      <NA>                                                      
## 2 Q922E6 Fastkd2      TISSUE SPECIFICITY: Ubiquitously expressed (PubMed:187717…
## 3 Q99KR8 Fuca2        <NA>                                                      
## 4 Q923G2 Polr2h       <NA>                                                      
## 5 P62254 Ube2g1 Ube2g <NA>                                                      
## 6 Q8K2G4 Bbs7 Bbs2l1  <NA>                                                      
## 7 Q505Q3 Fuca2        <NA>                                                      
## 8 Q5F239 Ube2g1       <NA>