The purpose of this vignette is to illustrate the feasibility of reflecting the material in the online tutorial for scvi-tools 0.20.0 in Bioconductor. The authors of the tutorial describe it as producing
a joint latent representation of cells, denoised data for both protein and RNA
Additional tasks include
integrat[ing] datasets, and comput[ing] differential expression of RNA and protein
The integration concerns the simultaneous analysis of two datasets from 10x genomcs.
In this vignette we carry out the bulk of the tutorial activities using R and Bioconductor, reaching to scvi-tools python code via basilisk.
The following chunk will acquire (and cache, using BiocFileCache) a preprocessed version of the 10k and 5k combined CITE-seq experiments from the scvi-tools data repository.
## AnnData object with n_obs × n_vars = 10849 × 4000
## obs: 'n_genes', 'percent_mito', 'n_counts', 'batch', '_scvi_labels', '_scvi_batch'
## var: 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm', 'highly_variable_nbatches'
## uns: '_scvi_manager_uuid', '_scvi_uuid', 'hvg', 'log1p'
## obsm: 'protein_expression'
## layers: 'counts'
The totalVI variational autoencoder was fit with these data. A fitted version is retrieved and cached using
## [34mINFO [0m File [35m/home/biocbuild/bbs-3.20-bioc/tmpdir/RtmpcbgIwC/vae2_ov/[0m[95mmodel.pt[0m
## already downloaded
## [34mINFO [0m Computing empirical prior initialization for protein background.
This is an instance of an S3 class, python.builtin.object
,
defined in the reticulate package.
## [1] "scvi.model._totalvi.TOTALVI"
## [2] "scvi.model.base._rnamixin.RNASeqMixin"
## [3] "scvi.model.base._vaemixin.VAEMixin"
## [4] "scvi.model.base._archesmixin.ArchesMixin"
## [5] "scvi.model.base._base_model.BaseModelClass"
## [6] "scvi._types.TunableMixin"
## [7] "python.builtin.object"
Some fields of interest that are directly available from the instance include an indicator of the trained state, the general parameters used to train, and the “anndata” (annotated data) object that includes the input counts and various results of preprocessing:
## [1] TRUE
## TotalVI Model with the following params:
## n_latent: 20, gene_dispersion: gene, protein_dispersion: protein, gene_likelihood: nb, latent_distribution: normal
## AnnData object with n_obs × n_vars = 10849 × 4000
## obs: 'n_genes', 'percent_mito', 'n_counts', 'batch', '_scvi_labels', '_scvi_batch'
## var: 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm', 'highly_variable_nbatches'
## uns: '_scvi_manager_uuid', '_scvi_uuid', 'hvg', 'log1p'
## obsm: 'protein_expression'
## layers: 'counts'
The structure of the VAE is reported using
This is quite voluminous and is provided in an appendix.
The negative “evidence lower bound” (ELBO) is a criterion that is minimized in order to produce a fitted autoencoder. The scvi-tools totalVAE elgorithm creates a nonlinear projection of the inputs to a 20-dimensional latent space, and a decoder that transforms object positions in the latent space to positions in the space of observations that are close to the original input positions.
The negative ELBO values are computed for samples from the training data and for “left out” validation samples. Details on the validation assessment would seem to be part of pytorch lightning. More investigation of scvi-tools code and doc are in order.
h = vae$history
npts = nrow(h$elbo_train)
plot(seq_len(npts), as.numeric(h$elbo_train[[1]]), ylim=c(1200,1400),
type="l", col="blue", main="Negative ELBO over training epochs",
ylab="neg. ELBO", xlab="epoch")
graphics::legend(300, 1360, lty=1, col=c("blue", "orange"), legend=c("training", "validation"))
graphics::lines(seq_len(npts), as.numeric(h$elbo_validation[[1]]), type="l", col="orange")
On a CPU, the following can take a long time.
NE = vae$get_normalized_expression(n_samples=25L,
return_mean=TRUE,
transform_batch=c("PBMC10k", "PBMC5k")
)
We provide the totalVI-based denoised quantities in
## rna_nmlzd prot_nmlzd
## [1,] 10849 10849
## [2,] 4000 14
Note that these have features as columns, samples (cells) as rows.
## [1] "AL645608.8" "HES4" "ISG15" "TTLL10" "TNFRSF18"
## [6] "TNFRSF4"
## [1] "CD3_TotalSeqB" "CD4_TotalSeqB" "CD8a_TotalSeqB" "CD14_TotalSeqB"
## [5] "CD15_TotalSeqB" "CD16_TotalSeqB"
We have stored a fully loaded anndata instance for retrieval to inspect the latent space and clustering produced by the tutorial notebook procedure. The images produced here do not agree exactly with what I see in the colab pages for 0.20.0. The process was run in Jetstream2, not in colab.
full = getTotalVI5k10kAdata()
# class distribution
cllabs = full$obs$leiden_totalVI
blabs = full$obs$batch
table(cllabs)
## cllabs
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 2287 1899 1138 866 787 660 637 587 461 375 334 261 260 208 59 19
## 16
## 11
um = full$obsm$get("X_umap")
dd = data.frame(umap1=um[,1], umap2=um[,2], clust=factor(cllabs), batch=blabs)
ggplot(dd, aes(x=umap1, y=umap2, colour=clust)) + geom_point(size=.05) +
guides(color = guide_legend(override.aes = list(size = 4)))
Effectiveness at accommodating the two-batch design is suggested by the mixed representation of the batches in all the Leiden clusters.
We focus on four of the ADT. Points (cells) in the UMAP projection given above are colored by the estimated abundance of the proteins quantified via (normalized) ADT abundance. Complementary expression of CD4 and CD8a is suggested by the configurations in the middle two panels.
pro4 = denoised$prot_nmlzd[,1:4]
names(pro4) = gsub("_.*", "", names(pro4))
wprot = cbind(dd, pro4)
mm = melt(wprot, id.vars=c("clust", "batch", "umap1", "umap2"))
utils::head(mm,3)
## clust batch umap1 umap2 variable value
## 1 0 PBMC10k 1.0286884 14.58188 CD3 0.6282051
## 2 0 PBMC10k 0.8443313 15.29502 CD3 11.9159498
## 3 0 PBMC10k 0.8057814 13.96269 CD3 8.9384489
We have shown that all the results of the totalVI application in the tutorial are readily accessible with utilities in scviR. Additional work on details of differential expression are present in the tutorial and can be explored by the interested reader/user.
The structure of the VAE is reported using
## TOTALVAE(
## (encoder): EncoderTOTALVI(
## (encoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=4016, out_features=256, bias=True)
## (1): BatchNorm1d(256, eps=0.001, momentum=0.01, affine=True, track_running_stats=True)
## (2): None
## (3): ReLU()
## (4): Dropout(p=0.2, inplace=False)
## )
## (Layer 1): Sequential(
## (0): Linear(in_features=258, out_features=256, bias=True)
## (1): BatchNorm1d(256, eps=0.001, momentum=0.01, affine=True, track_running_stats=True)
## (2): None
## (3): ReLU()
## (4): Dropout(p=0.2, inplace=False)
## )
## )
## )
## (z_mean_encoder): Linear(in_features=256, out_features=20, bias=True)
## (z_var_encoder): Linear(in_features=256, out_features=20, bias=True)
## (l_gene_encoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=4016, out_features=256, bias=True)
## (1): BatchNorm1d(256, eps=0.001, momentum=0.01, affine=True, track_running_stats=True)
## (2): None
## (3): ReLU()
## (4): Dropout(p=0.2, inplace=False)
## )
## )
## )
## (l_gene_mean_encoder): Linear(in_features=256, out_features=1, bias=True)
## (l_gene_var_encoder): Linear(in_features=256, out_features=1, bias=True)
## )
## (decoder): DecoderTOTALVI(
## (px_decoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=22, out_features=256, bias=True)
## (1): BatchNorm1d(256, eps=0.001, momentum=0.01, affine=True, track_running_stats=True)
## (2): None
## (3): ReLU()
## (4): Dropout(p=0.2, inplace=False)
## )
## )
## )
## (px_scale_decoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=278, out_features=4000, bias=True)
## (1): None
## (2): None
## (3): None
## (4): None
## )
## )
## )
## (px_scale_activation): Softmax(dim=-1)
## (py_back_decoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=22, out_features=256, bias=True)
## (1): BatchNorm1d(256, eps=0.001, momentum=0.01, affine=True, track_running_stats=True)
## (2): None
## (3): ReLU()
## (4): Dropout(p=0.2, inplace=False)
## )
## )
## )
## (py_back_mean_log_alpha): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=278, out_features=14, bias=True)
## (1): None
## (2): None
## (3): None
## (4): None
## )
## )
## )
## (py_back_mean_log_beta): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=278, out_features=14, bias=True)
## (1): None
## (2): None
## (3): None
## (4): None
## )
## )
## )
## (py_fore_decoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=22, out_features=256, bias=True)
## (1): BatchNorm1d(256, eps=0.001, momentum=0.01, affine=True, track_running_stats=True)
## (2): None
## (3): ReLU()
## (4): Dropout(p=0.2, inplace=False)
## )
## )
## )
## (py_fore_scale_decoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=278, out_features=14, bias=True)
## (1): None
## (2): None
## (3): ReLU()
## (4): None
## )
## )
## )
## (sigmoid_decoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=22, out_features=256, bias=True)
## (1): BatchNorm1d(256, eps=0.001, momentum=0.01, affine=True, track_running_stats=True)
## (2): None
## (3): ReLU()
## (4): Dropout(p=0.2, inplace=False)
## )
## )
## )
## (px_dropout_decoder_gene): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=278, out_features=4000, bias=True)
## (1): None
## (2): None
## (3): None
## (4): None
## )
## )
## )
## (py_background_decoder): FCLayers(
## (fc_layers): Sequential(
## (Layer 0): Sequential(
## (0): Linear(in_features=278, out_features=14, bias=True)
## (1): None
## (2): None
## (3): None
## (4): None
## )
## )
## )
## )
## )
## signature: (*input, **kwargs)
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB 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
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] reshape2_1.4.4 ggplot2_3.5.1
## [3] scviR_1.6.0 SingleCellExperiment_1.28.0
## [5] SummarizedExperiment_1.36.0 Biobase_2.66.0
## [7] GenomicRanges_1.58.0 GenomeInfoDb_1.42.0
## [9] IRanges_2.40.0 S4Vectors_0.44.0
## [11] BiocGenerics_0.52.0 MatrixGenerics_1.18.0
## [13] matrixStats_1.4.1 shiny_1.9.1
## [15] basilisk_1.18.0 reticulate_1.39.0
## [17] BiocStyle_2.34.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 gridExtra_2.3 rlang_1.1.4
## [4] magrittr_2.0.3 scater_1.34.0 compiler_4.4.1
## [7] RSQLite_2.3.7 dir.expiry_1.14.0 png_0.1-8
## [10] vctrs_0.6.5 stringr_1.5.1 pkgconfig_2.0.3
## [13] crayon_1.5.3 fastmap_1.2.0 magick_2.8.5
## [16] dbplyr_2.5.0 XVector_0.46.0 labeling_0.4.3
## [19] scuttle_1.16.0 utf8_1.2.4 promises_1.3.0
## [22] rmarkdown_2.28 UCSC.utils_1.2.0 ggbeeswarm_0.7.2
## [25] tinytex_0.53 purrr_1.0.2 bit_4.5.0
## [28] xfun_0.48 zlibbioc_1.52.0 cachem_1.1.0
## [31] beachmat_2.22.0 jsonlite_1.8.9 blob_1.2.4
## [34] highr_0.11 later_1.3.2 DelayedArray_0.32.0
## [37] BiocParallel_1.40.0 irlba_2.3.5.1 parallel_4.4.1
## [40] R6_2.5.1 stringi_1.8.4 bslib_0.8.0
## [43] RColorBrewer_1.1-3 limma_3.62.0 jquerylib_0.1.4
## [46] Rcpp_1.0.13 bookdown_0.41 knitr_1.48
## [49] httpuv_1.6.15 Matrix_1.7-1 tidyselect_1.2.1
## [52] abind_1.4-8 yaml_2.3.10 viridis_0.6.5
## [55] codetools_0.2-20 curl_5.2.3 plyr_1.8.9
## [58] lattice_0.22-6 tibble_3.2.1 basilisk.utils_1.18.0
## [61] withr_3.0.2 evaluate_1.0.1 BiocFileCache_2.14.0
## [64] pillar_1.9.0 BiocManager_1.30.25 filelock_1.0.3
## [67] generics_0.1.3 munsell_0.5.1 scales_1.3.0
## [70] xtable_1.8-4 glue_1.8.0 pheatmap_1.0.12
## [73] tools_4.4.1 BiocNeighbors_2.0.0 ScaledMatrix_1.14.0
## [76] grid_4.4.1 colorspace_2.1-1 GenomeInfoDbData_1.2.13
## [79] beeswarm_0.4.0 BiocSingular_1.22.0 vipor_0.4.7
## [82] cli_3.6.3 rsvd_1.0.5 fansi_1.0.6
## [85] S4Arrays_1.6.0 viridisLite_0.4.2 dplyr_1.1.4
## [88] gtable_0.3.6 sass_0.4.9 digest_0.6.37
## [91] SparseArray_1.6.0 ggrepel_0.9.6 farver_2.1.2
## [94] memoise_2.0.1 htmltools_0.5.8.1 lifecycle_1.0.4
## [97] httr_1.4.7 statmod_1.5.0 mime_0.12
## [100] bit64_4.5.2