Generating Samples of Gene Expression Data with Variational Autoencoders
Dongmin Jung
December 15, 2021
1 Introduction
A fundamental problem in biomedical research is the low number of observations, mostly due to a lack of available biosamples, prohibitive costs, or ethical reasons. By augmenting a few real observations with artificially generated samples, their analysis could lead to more robust and higher reproducibility. One possible solution to the problem is the use of generative models, which are statistical models of data that attempt to capture the entire probability distribution from the observations. Using the variational autoencoder (VAE), a well-known deep generative model, this package is aimed to generate samples with gene expression data, especially for single-cell RNA-seq data. Furthermore, the VAE can use conditioning to produce specific cell types or subpopulations. The conditional VAE (CVAE) allows us to create targeted samples rather than completely random ones.
Autoencoders are an unsupervised neural networks that perform data compression from multidimensional to a preferred dimensionality. They reconstruct input data using the hidden layer weights calculated by encoding. The basic idea of an autoencoder is to have an output layer with the same dimensionality as the inputs. The idea is to try to reconstruct each dimension exactly by passing it through the network. It is common but not necessary for an autoencoder to have a symmetric architecture between the input and output. The number of units in the middle layer is typically fewer than that in the input or output. After training an autoencoder, it is not necessary to use both the encoder and decoder portions. For example, when using the approach for dimensionality reduction, one can use the encoder portion in order to create the reduced representations of the data. The reconstructions of the decoder might not be required at all. As a result, an autoencoder is capable of performing dimension reduction. The objective function of this neural network encompasses reconstruction loss. The loss function uses the sum of squared differences between the input and the output in order to force the output to be as similar as possible to the input. Also, the cross-entropy can used as a loss function for quantifying the difference between two probability distributions.
Another interesting application of the autoencoder is one in which we use only the decoder portion of the network. Variational autoencoders are based on Bayesian inference in which the compressed representation follows probability distribution. This constraint differentiates the VAE from standard autoencoder. The VAE can generate new data while conventional autoencoders fail. For example, one might add a term to the loss function to enforce the fact that the hidden variables are drawn from a Gaussian distribution. Then, one might repeatedly draw samples from this Gaussian distribution and use only the decoder portion of the network in order to generate samples of the original data. In this autoencoder, bottleneck vector (latent vector) is replaced by two vectors, namely, mean vector and standard deviation vector. The overall loss function \(J = L + \lambda R\) of the VAE is expressed as a weighted sum of the reconstruction loss \(L\) and the regularization loss \(R\), where \(\lambda > 0\) is the regularization parameter. The term “variational” comes from the close relationship between the regularization and the variational inference method in statistics. One can use a variety of choices for the reconstruction error, and we will use the binary cross-entropy loss between the input and output. The regularization loss is simply the Kullback-Leibler divergence measure of the conditional distributions of the hidden representations of particular points with respect to the standard multivariate Gaussian distribution. Small values of \(\lambda\) will favor exact reconstruction, and the approach will behave like a traditional autoencoder.
One can apply conditioning to variational autoencoders in order to obtain some interesting results. The basic idea in the conditional variational autoencoder is to add an additional conditional input. From an implementation perspective, we can encode category information as a one-hot representation, indicating to the model which class is at the input. One can use an autoencoder for embedding multimodal data in a joint latent space. Multimodal data is essentially data in which the input features are heterogeneous. In addition, by separating the samples into different classes, the data points within the same category become more similar, enhancing the modeling capacity and sample quality of the CVAE.
2 Example
2.1 VAE
Consider artificial data. The data consist of 1000 genes and three groups of 100 samples. Each group has 100 differentially expressed genes.
if (keras::is_keras_available() & reticulate::py_available()) {
library(VAExprs)
### simulate differentially expressed genes
set.seed(1)
g <- 3
n <- 100
m <- 1000
mu <- 5
sigma <- 5
mat <- matrix(rnorm(n*m*g, mu, sigma), m, n*g)
rownames(mat) <- paste0("gene", seq_len(m))
colnames(mat) <- paste0("cell", seq_len(n*g))
group <- factor(sapply(seq_len(g), function(x) {
rep(paste0("group", x), n)
}))
names(group) <- colnames(mat)
mu_upreg <- 6
sigma_upreg <- 10
deg <- 100
for (i in seq_len(g)) {
mat[(deg*(i-1) + 1):(deg*i), group == paste0("group", i)] <-
mat[1:deg, group==paste0("group", i)] + rnorm(deg, mu_upreg, sigma_upreg)
}
# positive expression only
mat[mat < 0] <- 0
x_train <- as.matrix(t(mat))
# heatmap
heatmap(mat, Rowv = NA, Colv = NA,
col = colorRampPalette(c('green', 'red'))(100),
scale = "none")
}## Loading required package: keras## Loading required package: mclust## Package 'mclust' version 5.4.8
## Type 'citation("mclust")' for citing this R package in publications.
The VAE model can be built by using the function “fit_vae” with gene expression data and the cell annotation from the object “sce”. The overall loss function of the VAE is expressed as a weighted sum of the reconstruction loss and the regularization loss. The reconstruction loss is the binary cross-entropy loss between the input and output and the regularization loss is simply the Kullback-Leibler divergence measure. Note that the same dataset is used for training and validation.
if (keras::is_keras_available() & reticulate::py_available()) {
# model parameters
batch_size <- 32
original_dim <- 1000
intermediate_dim <- 512
epochs <- 100
# VAE
vae_result <- fit_vae(x_train = x_train, x_val = x_train,
encoder_layers = list(layer_input(shape = c(original_dim)),
layer_dense(units = intermediate_dim,
activation = "relu")),
decoder_layers = list(layer_dense(units = intermediate_dim,
activation = "relu"),
layer_dense(units = original_dim,
activation = "sigmoid")),
epochs = epochs, batch_size = batch_size,
use_generator = FALSE,
callbacks = keras::callback_early_stopping(
monitor = "val_loss",
patience = 10,
restore_best_weights = TRUE))
}## Loaded Tensorflow version 2.4.1## normalizing...## training...The function “plot_vae” draws the plot for model architecture.
if (keras::is_keras_available() & reticulate::py_available()) {
# model architecture
plot_vae(vae_result$model)
}The function “gen_exprs” can generate samples with expression data by using the trained model.
if (keras::is_keras_available() & reticulate::py_available()) {
# sample generation
set.seed(1)
gen_sample_result <- gen_exprs(vae_result, num_samples = 100)
# heatmap
heatmap(cbind(t(x_train), t(gen_sample_result$x_gen)),
col = colorRampPalette(c('green', 'red'))(100),
Rowv=NA)
}## generating...## post-processing...
The function “plot_aug” uses reduced dimension plots for augmented data visualization.
if (keras::is_keras_available() & reticulate::py_available()) {
# plot for augmented data
plot_aug(gen_sample_result, "PCA")
}
2.2 CVAE
The “yan” data set is single-cell RNA sequencing data with 20214 genes and 90 cells from human preimplantation embryos and embryonic stem cells at different passages. The rows in the dataset correspond to genes and columns correspond to cells. The “SingleCellExperiment” class can be used to store and manipulate single-cell genomics data. It extends the “RangedSummarizedExperiment” class and follows similar conventions. The object “sce” can be created by the data “yan” with cell type annotation “ann”.
if (keras::is_keras_available() & reticulate::py_available()) {
library(VAExprs)
library(SC3)
library(SingleCellExperiment)
# create a SingleCellExperiment object
sce <- SingleCellExperiment::SingleCellExperiment(
assays = list(counts = as.matrix(yan)),
colData = ann
)
# define feature names in feature_symbol column
rowData(sce)$feature_symbol <- rownames(sce)
# remove features with duplicated names
sce <- sce[!duplicated(rowData(sce)$feature_symbol), ]
# remove genes that are not expressed in any samples
sce <- sce[which(rowMeans(assay(sce)) > 0),]
dim(assay(sce))
# model parameters
batch_size <- 32
original_dim <- 19595
intermediate_dim <- 256
epochs <- 100
# model
cvae_result <- fit_vae(object = sce,
encoder_layers = list(layer_input(shape = c(original_dim)),
layer_dense(units = intermediate_dim,
activation = "relu")),
decoder_layers = list(layer_dense(units = intermediate_dim,
activation = "relu"),
layer_dense(units = original_dim,
activation = "sigmoid")),
epochs = epochs, batch_size = batch_size,
use_generator = TRUE,
callbacks = keras::callback_early_stopping(
monitor = "loss",
patience = 20,
restore_best_weights = TRUE))
# model architecture
plot_vae(cvae_result$model)
}## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: GenomicRanges## Loading required package: stats4## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following object is masked from 'package:keras':
##
## normalize## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## Filter, Find, Map, Position, Reduce, anyDuplicated, append,
## as.data.frame, basename, cbind, colnames, dirname, do.call,
## duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
## lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
## pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which.max, which.min## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:base':
##
## I, expand.grid, unname## Loading required package: IRanges## Loading required package: GenomeInfoDb## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians## pre-processing...## normalizing...## training...if (keras::is_keras_available() & reticulate::py_available()) {
# sample generation
set.seed(1)
gen_sample_result <- gen_exprs(cvae_result, 100,
batch_size, use_generator = TRUE)
# plot for augmented data
plot_aug(gen_sample_result, "PCA")
}## generating...## post-processing...
3 Session information
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
##
## 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
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0
## [3] Biobase_2.54.0 GenomicRanges_1.46.1
## [5] GenomeInfoDb_1.30.0 IRanges_2.28.0
## [7] S4Vectors_0.32.3 BiocGenerics_0.40.0
## [9] MatrixGenerics_1.6.0 matrixStats_0.61.0
## [11] SC3_1.22.0 VAExprs_1.0.1
## [13] mclust_5.4.8 keras_2.7.0
##
## loaded via a namespace (and not attached):
## [1] DeepPINCS_1.2.1 BiocParallel_1.28.3
## [3] tfruns_1.5.0 SnowballC_0.7.0
## [5] ggplot2_3.3.5 scater_1.22.0
## [7] digest_0.6.29 foreach_1.5.1
## [9] htmltools_0.5.2 viridis_0.6.2
## [11] float_0.2-6 fansi_0.5.0
## [13] magrittr_2.0.1 ScaledMatrix_1.2.0
## [15] tm_0.7-8 cluster_2.1.2
## [17] doParallel_1.0.16 ROCR_1.0-11
## [19] colorspace_2.0-2 rrcov_1.6-0
## [21] rvest_1.0.2 rappdirs_0.3.3
## [23] ggrepel_0.9.1 sylly_0.1-6
## [25] WriteXLS_6.3.0 xfun_0.29
## [27] dplyr_1.0.7 crayon_1.4.2
## [29] RCurl_1.98-1.5 CatEncoders_0.1.1
## [31] jsonlite_1.7.2 zeallot_0.1.0
## [33] iterators_1.0.13 glue_1.5.1
## [35] stopwords_2.3 gtable_0.3.0
## [37] zlibbioc_1.40.0 XVector_0.34.0
## [39] DelayedArray_0.20.0 webchem_1.1.2
## [41] BiocSingular_1.10.0 DEoptimR_1.0-9
## [43] mlapi_0.1.0 scales_1.1.1
## [45] mvtnorm_1.1-3 pheatmap_1.0.12
## [47] DBI_1.1.1 rngtools_1.5.2
## [49] data.tree_1.0.0 Rcpp_1.0.7
## [51] xtable_1.8-4 viridisLite_0.4.0
## [53] gradDescent_3.0 reticulate_1.22
## [55] proxy_0.4-26 text2vec_0.6
## [57] rsvd_1.0.5 matlab_1.0.2
## [59] htmlwidgets_1.5.4 httr_1.4.2
## [61] DiagrammeR_1.0.6.1 RColorBrewer_1.1-2
## [63] ellipsis_0.3.2 pkgconfig_2.0.3
## [65] farver_2.1.0 scuttle_1.4.0
## [67] sass_0.4.0 utf8_1.2.2
## [69] here_1.0.1 later_1.3.0
## [71] tidyselect_1.1.1 labeling_0.4.2
## [73] rlang_0.4.12 PRROC_1.3.1
## [75] munsell_0.5.0 tools_4.1.2
## [77] visNetwork_2.1.0 generics_0.1.1
## [79] rsparse_0.5.0 evaluate_0.14
## [81] stringr_1.4.0 fastmap_1.1.0
## [83] yaml_2.2.1 textstem_0.1.4
## [85] RhpcBLASctl_0.21-247.1 knitr_1.36
## [87] robustbase_0.93-9 lgr_0.4.3
## [89] purrr_0.3.4 doRNG_1.8.2
## [91] sparseMatrixStats_1.6.0 mime_0.12
## [93] koRpus_0.13-8 whisker_0.4
## [95] slam_0.1-49 xml2_1.3.3
## [97] tokenizers_0.2.1 compiler_4.1.2
## [99] rstudioapi_0.13 beeswarm_0.4.0
## [101] png_0.1-7 e1071_1.7-9
## [103] koRpus.lang.en_0.1-4 tibble_3.1.6
## [105] pcaPP_1.9-74 bslib_0.3.1
## [107] stringi_1.7.6 highr_0.9
## [109] lattice_0.20-45 Matrix_1.4-0
## [111] tensorflow_2.7.0 vctrs_0.3.8
## [113] stringdist_0.9.8 pillar_1.6.4
## [115] lifecycle_1.0.1 jquerylib_0.1.4
## [117] BiocNeighbors_1.12.0 data.table_1.14.2
## [119] cowplot_1.1.1 bitops_1.0-7
## [121] irlba_2.3.5 httpuv_1.6.4
## [123] sylly.en_0.1-3 ttgsea_1.2.1
## [125] R6_2.5.1 promises_1.2.0.1
## [127] gridExtra_2.3 vipor_0.4.5
## [129] codetools_0.2-18 assertthat_0.2.1
## [131] rprojroot_2.0.2 GenomeInfoDbData_1.2.7
## [133] parallel_4.1.2 grid_4.1.2
## [135] beachmat_2.10.0 class_7.3-19
## [137] rmarkdown_2.11 DelayedMatrixStats_1.16.0
## [139] shiny_1.7.1 NLP_0.2-1
## [141] base64enc_0.1-3 ggbeeswarm_0.6.04 References
Aggarwal, C. C. (2018). Neural Networks and Deep Learning. Springer.
Al-Jabery, K., Obafemi-Ajayi, T., Olbricht, G., & Wunsch, D. (2019). Computational Learning Approaches to Data Analytics in Biomedical Applications. Academic Press.
Cinelli, L. P., Marins, M. A., da Silva, E. A. B., & Netto, S. L. (2021). Variational Methods for Machine Learning with Applications to Deep Networks. Springer.
Das, H., Pradhan, C., & Dey, N. (2020). Deep Learning for Data Analytics: Foundations, Biomedical Applications, and Challenges. Academic Press.
Marouf, M., Machart, P., Bansal, V., Kilian, C., Magruder, D. S., Krebs, C. F., & Bonn, S. (2020). Realistic in silico generation and augmentation of single-cell RNA-seq data using generative adversarial networks. Nature communications, 11(1), 1-12.
Pedrycz, W., & Chen, S. M. (Eds.). (2020). Deep Learning: Concepts and Architectures. Springer.
Yan, W. Q. (2020). Computational Methods for Deep Learning: Theoretic, Practice and Applications. Springer.