Computing RNA velocity in a Bioconductor framework
28 October 2020
Package
velociraptor 1.0.0
Contents
1 Overview
This package provides a lightweight interface between the Bioconductor SingleCellExperiment data structure
and the scvelo Python package for RNA velocity calculations.
The interface is comparable to that of many other SingleCellExperiment-compatible functions,
allowing users to plug in RNA velocity calculations into the existing Bioconductor analysis framework.
To demonstrate, we will use a data set from Hermann et al. (2018), provided via the scRNAseq package.
This data set contains gene-wise estimates of spliced and unspliced UMI counts for 2,325 mouse spermatogenic cells.
library(scRNAseq)
sce <- HermannSpermatogenesisData()
sce## class: SingleCellExperiment
## dim: 54448 2325
## metadata(0):
## assays(2): spliced unspliced
## rownames(54448): ENSMUSG00000102693.1 ENSMUSG00000064842.1 ...
## ENSMUSG00000064369.1 ENSMUSG00000064372.1
## rowData names(0):
## colnames(2325): CCCATACTCCGAAGAG AATCCAGTCATCTGCC ... ATCCACCCACCACCAG
## ATTGGTGGTTACCGAT
## colData names(1): celltype
## reducedDimNames(0):
## altExpNames(0):2 Downsampling for demonstration
The full data set requires up to 12 GB of memory for the example usage presented in this vignette. For demonstration purposes, we downsample the data set to the first 500 cells. Feel free to skip this downsampling step if you have access to sufficient memory.
sce <- sce[, 1:500]3 Basic workflow
We assume that feature selection has already been performed by the user using any method (see here for some suggestions). We will use the log-expression method from the scran method here.
library(scuttle)
sce <- logNormCounts(sce, assay.type=1)
library(scran)
dec <- modelGeneVar(sce)
top.hvgs <- getTopHVGs(dec, n=2000)We can plug these choices into the scvelo() function with our SingleCellExperiment object.
Here, we will use the "spliced" count matrix as a proxy for the typical exonic count matrix.
While this count matrix is not required for the velocity estimation in itself, it is used internally to retrieve the principal component space that is used by scvelo to find nearest neighbors.
There are some subtle differences between the spliced count matrix and the typical exonic count matrix - see ?scvelo for some commentary about this -
but the spliced counts are generally a satisfactory replacement.
By default, scvelo() uses the steady-state approach to estimate velocities.
The stochastic and dynamical models implemented in scvelo are accessible as well, by modifying the mode argument.
library(velociraptor)
velo.out <- scvelo(sce, subset.row=top.hvgs, assay.X="spliced")
velo.out## class: SingleCellExperiment
## dim: 2000 500
## metadata(4): neighbors velocity_params velocity_graph
## velocity_graph_neg
## assays(6): X spliced ... Mu velocity
## rownames(2000): ENSMUSG00000117819.1 ENSMUSG00000081984.3 ...
## ENSMUSG00000022965.8 ENSMUSG00000094660.2
## rowData names(3): velocity_gamma velocity_r2 velocity_genes
## colnames(500): CCCATACTCCGAAGAG AATCCAGTCATCTGCC ... CACCTTGTCGTAGGAG
## TTCCCAGAGACTAAGT
## colData names(7): velocity_self_transition root_cells ...
## velocity_confidence velocity_confidence_transition
## reducedDimNames(1): X_pca
## altExpNames(0):The scvelo() function produces a SingleCellExperiment containing all of the outputs of the calculation in Python.
Of particular interest is the velocity_pseudotime vector that captures the relative progression of each cell
along the biological process driving the velocity vectors.
We can visualize this effect below in a \(t\)-SNE plot generated by scater on the top HVGs.
library(scater)
set.seed(100)
sce <- runPCA(sce, subset_row=top.hvgs)
sce <- runTSNE(sce, dimred="PCA")
sce$velocity_pseudotime <- velo.out$velocity_pseudotime
plotTSNE(sce, colour_by="velocity_pseudotime")
It is also straightforward to embed the velocity vectors into our desired low-dimensional space, as shown below for the \(t\)-SNE coordinates. This uses a grid-based approach to summarize the per-cell vectors into local representatives for effective visualization.
embedded <- embedVelocity(reducedDim(sce, "TSNE"), velo.out)
grid.df <- gridVectors(reducedDim(sce, "TSNE"), embedded)
library(ggplot2)
plotTSNE(sce, colour_by="velocity_pseudotime") +
geom_segment(data=grid.df, mapping=aes(x=start.1, y=start.2,
xend=end.1, yend=end.2), arrow=arrow(length=unit(0.05, "inches")))
And that’s it, really.
4 Advanced options
scvelo() interally performs a PCA step that we can bypass by supplying our own PC coordinates.
Indeed, it is often the case that we have already performed PCA in the earlier analysis steps,
so we can just re-use those results to (i) save time and (ii) improve consistency with the other steps.
Here, we computed the PCA coordinates in runPCA() above, so let’s just recycle that:
# Only setting assay.X= for the initial AnnData creation,
# it is not actually used in any further steps.
velo.out2 <- scvelo(sce, assay.X=1, subset.row=top.hvgs, use.dimred="PCA")
velo.out2## class: SingleCellExperiment
## dim: 2000 500
## metadata(4): neighbors velocity_params velocity_graph
## velocity_graph_neg
## assays(6): X spliced ... Mu velocity
## rownames(2000): ENSMUSG00000117819.1 ENSMUSG00000081984.3 ...
## ENSMUSG00000022965.8 ENSMUSG00000094660.2
## rowData names(3): velocity_gamma velocity_r2 velocity_genes
## colnames(500): CCCATACTCCGAAGAG AATCCAGTCATCTGCC ... CACCTTGTCGTAGGAG
## TTCCCAGAGACTAAGT
## colData names(7): velocity_self_transition root_cells ...
## velocity_confidence velocity_confidence_transition
## reducedDimNames(1): X_pca
## altExpNames(0):We also provide an option to use the scvelo pipeline without modification,
i.e., relying on their normalization and feature selection.
This sacrifices consistency with other Bioconductor workflows but enables perfect mimicry of a pure Python-based analysis.
In this case, arguments like subset.row= are simply ignored.
velo.out3 <- scvelo(sce, assay.X=1, use.theirs=TRUE)
velo.out3## class: SingleCellExperiment
## dim: 54448 500
## metadata(5): pca neighbors velocity_params velocity_graph
## velocity_graph_neg
## assays(6): X spliced ... Mu velocity
## rownames(54448): ENSMUSG00000102693.1 ENSMUSG00000064842.1 ...
## ENSMUSG00000064369.1 ENSMUSG00000064372.1
## rowData names(4): velocity_gamma velocity_r2 velocity_genes varm
## colnames(500): CCCATACTCCGAAGAG AATCCAGTCATCTGCC ... CACCTTGTCGTAGGAG
## TTCCCAGAGACTAAGT
## colData names(11): initial_size_spliced initial_size_unspliced ...
## velocity_confidence velocity_confidence_transition
## reducedDimNames(1): X_pca
## altExpNames(0):Advanced users can tinker with the settings of individual scvelo steps
by setting named lists of arguments in the scvelo.params= argument.
For example, to tinker with the behavior of the recover_dynamics step, we could do:
velo.out4 <- scvelo(sce, assay.X=1, subset.row=top.hvgs,
scvelo.params=list(recover_dynamics=list(max_iter=20)))
velo.out4## class: SingleCellExperiment
## dim: 2000 500
## metadata(4): neighbors velocity_params velocity_graph
## velocity_graph_neg
## assays(6): X spliced ... Mu velocity
## rownames(2000): ENSMUSG00000117819.1 ENSMUSG00000081984.3 ...
## ENSMUSG00000022965.8 ENSMUSG00000094660.2
## rowData names(3): velocity_gamma velocity_r2 velocity_genes
## colnames(500): CCCATACTCCGAAGAG AATCCAGTCATCTGCC ... CACCTTGTCGTAGGAG
## TTCCCAGAGACTAAGT
## colData names(7): velocity_self_transition root_cells ...
## velocity_confidence velocity_confidence_transition
## reducedDimNames(1): X_pca
## altExpNames(0):5 Session information
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-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] scater_1.18.0 ggplot2_3.3.2
## [3] velociraptor_1.0.0 scran_1.18.0
## [5] scuttle_1.0.0 scRNAseq_2.3.17
## [7] SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0
## [9] Biobase_2.50.0 GenomicRanges_1.42.0
## [11] GenomeInfoDb_1.26.0 IRanges_2.24.0
## [13] S4Vectors_0.28.0 BiocGenerics_0.36.0
## [15] MatrixGenerics_1.2.0 matrixStats_0.57.0
## [17] knitr_1.30 BiocStyle_2.18.0
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 ggbeeswarm_0.6.0
## [3] colorspace_1.4-1 ellipsis_0.3.1
## [5] bluster_1.0.0 XVector_0.30.0
## [7] BiocNeighbors_1.8.0 farver_2.0.3
## [9] bit64_4.0.5 interactiveDisplayBase_1.28.0
## [11] AnnotationDbi_1.52.0 xml2_1.3.2
## [13] sparseMatrixStats_1.2.0 jsonlite_1.7.1
## [15] Rsamtools_2.6.0 dbplyr_1.4.4
## [17] shiny_1.5.0 BiocManager_1.30.10
## [19] compiler_4.0.3 httr_1.4.2
## [21] dqrng_0.2.1 basilisk_1.2.0
## [23] assertthat_0.2.1 Matrix_1.2-18
## [25] fastmap_1.0.1 lazyeval_0.2.2
## [27] limma_3.46.0 later_1.1.0.1
## [29] BiocSingular_1.6.0 htmltools_0.5.0
## [31] prettyunits_1.1.1 tools_4.0.3
## [33] rsvd_1.0.3 igraph_1.2.6
## [35] gtable_0.3.0 glue_1.4.2
## [37] GenomeInfoDbData_1.2.4 dplyr_1.0.2
## [39] rappdirs_0.3.1 Rcpp_1.0.5
## [41] vctrs_0.3.4 Biostrings_2.58.0
## [43] ExperimentHub_1.16.0 rtracklayer_1.50.0
## [45] zellkonverter_1.0.0 DelayedMatrixStats_1.12.0
## [47] xfun_0.18 stringr_1.4.0
## [49] beachmat_2.6.0 mime_0.9
## [51] lifecycle_0.2.0 irlba_2.3.3
## [53] ensembldb_2.14.0 statmod_1.4.35
## [55] XML_3.99-0.5 AnnotationHub_2.22.0
## [57] edgeR_3.32.0 scales_1.1.1
## [59] zlibbioc_1.36.0 basilisk.utils_1.2.0
## [61] hms_0.5.3 promises_1.1.1
## [63] ProtGenerics_1.22.0 AnnotationFilter_1.14.0
## [65] yaml_2.2.1 curl_4.3
## [67] gridExtra_2.3 memoise_1.1.0
## [69] reticulate_1.18 biomaRt_2.46.0
## [71] stringi_1.5.3 RSQLite_2.2.1
## [73] BiocVersion_3.12.0 GenomicFeatures_1.42.0
## [75] filelock_1.0.2 BiocParallel_1.24.0
## [77] rlang_0.4.8 pkgconfig_2.0.3
## [79] bitops_1.0-6 evaluate_0.14
## [81] lattice_0.20-41 purrr_0.3.4
## [83] labeling_0.4.2 GenomicAlignments_1.26.0
## [85] cowplot_1.1.0 bit_4.0.4
## [87] tidyselect_1.1.0 magrittr_1.5
## [89] bookdown_0.21 R6_2.4.1
## [91] magick_2.5.0 generics_0.0.2
## [93] DelayedArray_0.16.0 DBI_1.1.0
## [95] withr_2.3.0 pillar_1.4.6
## [97] RCurl_1.98-1.2 tibble_3.0.4
## [99] crayon_1.3.4 BiocFileCache_1.14.0
## [101] rmarkdown_2.5 viridis_0.5.1
## [103] progress_1.2.2 locfit_1.5-9.4
## [105] grid_4.0.3 blob_1.2.1
## [107] digest_0.6.27 xtable_1.8-4
## [109] httpuv_1.5.4 munsell_0.5.0
## [111] openssl_1.4.3 viridisLite_0.3.0
## [113] beeswarm_0.2.3 vipor_0.4.5
## [115] askpass_1.1References
Hermann, Brian P, Keren Cheng, Anukriti Singh, Lorena Roa-De La Cruz, Kazadi N Mutoji, I-Chung Chen, Heidi Gildersleeve, et al. 2018. “The Mammalian Spermatogenesis Single-Cell Transcriptome, from Spermatogonial Stem Cells to Spermatids.” Cell Rep. 25 (6):1650–1667.e8.