--- title: 'RNA-seq vignette: dimensionality reduction with omicsGMF' author: - name: Alexandre Segers bibliography: sgdGMF.bib date: "14/01/2025" output: BiocStyle::html_document: toc: true toc_depth: 3 BiocStyle::pdf_document: default package: omicsGMF abstract: | RNA-seq vignette for the omicsGMF package. This vignette aims to provide a detailed description of a matrix factorization on RNA-seq, which can be used for dimensionality reduction and visualization of RNA-seq data. vignette: > %\VignetteIndexEntry{RNASeq-vignette: omicsGMF} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, echo = FALSE} library(knitr) ``` # Introduction `omicsGMF` is an R package for generalized matrix factorization and missing value imputation in omics data. It is designed for dimensionality reduction and visualization, specifically handling count data and missing values efficiently. Unlike conventional PCA, `omicsGMF` does not require log-transformation of RNA-seq data or prior imputation of proteomics data. A key advantage of `omicsGMF` is its ability to control for known sample- and feature-level covariates, such as batch effects. This improves downstream analyses like clustering. Additionally, omicsGMF includes model selection to optimize the number of latent confounders, ensuring an optimal dimensionality for analysis. Its stochastic optimization algorithms allow it to remain fast while handling these complex data structures. `omicsGMF` builds on the `sgdGMF` framework provided in the `sgdGMF` CRAN package, and provides easy integration with `SingleCellExperiment`, `SummarizedExperiment`, and `QFeature` classes, with adapted default values for the optimization arguments when dealing with omics data. All details about the `sgdGMF` framework, such as the adaptive learning rates, exponential gradient averaging and subsampling of the data are described in our preprint [@Castiglione2024]. There, we show the use of the `sgdGMF-framework` on single-cell RNA-seq data. In our newest preprint [@Segers2025], we show how `omicsGMF` can be used to visualize (single-cell) proteomics data and impute missing values. This vignette provides a step-by-step workflow for using `omicsGMF` for dimensionality reduction of omics data. The main function are: 1. `runCVGMF` or `calculateCVGMF` performs cross-validation to determine the optimal number of latent confounders. These results can be visualized using `plotCV`. This cross-validation avoids arbitrarily choosing `ncomponents`, but requires some computational time. An alternative is `calculateRankGMF`, which performs an eigenvalue decomposition on the deviance residuals. This allows for model selection based on a scree plot using `plotRank`, for example using the elbow method. 2. `runGMF` or `calculateGMF` estimates the latent confounders and the rotation matrix, and estimates the respective parameters of the sample-level and feature-level covariates. 3. `plotGMF` plots the samples using its decomposition. 4. `imputeGMF` creates a new assay with missing values imputed using the estimates of `runGMF`. We here apply `omicsGMF` on RNA-seq data. # Package installation `sgdGMF` can be installed through CRAN. `omicsGMF` can be installed from github, and will be soon available through Bioconductor. ```{r, eval=FALSE} if(!requireNamespace("sgdGMF", quietly = TRUE)) install.packages("sgdGMF") BiocManager::install("omicsGMF") ``` ```{r, echo = TRUE, warning=FALSE, message=FALSE} library(sgdGMF) library(omicsGMF) library(dplyr) library(scuttle) set.seed(100) ``` # RNA-seq analysis To perform dimensionality reduction on RNA-seq data, one can use the original count matrices, without normalizing or log-transforming the sequencing counts to the Gaussian scale. By using `family = poisson()`, `omicsGMF` optimizes the dimensionality reduction with respect to the likelihood of the Poisson family. Further, by including a known covariate matrix, X, `omicsGMF` corrects for known confounders in the dimensionality reduction. First, we simulate a small dataset using the `scuttle` package. For sake of exposition, we will further account for the `Treatment` covariate from the `colData`. `omicsGMF` can internally correct for these treatment effects, and therefore does not require prior correction with other tools. ```{r} example_sce <- mockSCE(ncells = 20, ngenes = 50) X <- model.matrix(~Treatment, colData(example_sce)) ``` A recommended step is to estimate the optimal dimensionality in the model by using cross-validation. This cross-validation masks a proportion of the values as missing, and tries to reconstruct these. Using the out-of-sample deviances, one can estimate the optimal dimensionality of the latent space. This cross-validation can be done with the `runCVGMF` or `calculateCVGMF` function, which builds on the sgdgmf.cv function from the `sgdGMF` package. Although the `sgdGMF` framework allows great flexibility regarding the optimization algorithm, sensible default values are here introduced for omics data. One should choose the correct distribution family (`family`), the number of components in the dimensionality reduction for which the cross-validation is run (`ncomponents`), and the known covariate matrices to account for (`X` and `Z`). Also, one should select the right assay that is used for dimensionality reduction (`exprs_values` or `assay.type`). Visualization of the cross-validation results can be done using `plotCV`. In case that multiple cross-validation results are available in the `metadata`, it is possible to visualize these by giving all names of the metadata slots. The optimal dimensionality is the one that has the lowest out-of-sample deviances. ```{r} example_sce <- runCVGMF( example_sce, X = X, # Covariate matrix exprs_values = "counts", # Use raw counts (no normalization) family = poisson(), # Poisson model for RNA-seq count data ncomponents = 1:5, # Test components from 1 to 5 ntop = 50 # Use top 50 most variable genes ) metadata(example_sce)$cv_GMF %>% group_by(ncomp) %>% summarise(mean_dev = mean(dev), mean_aic = mean(aic), mean_bic = mean(bic), mean_mae = mean(mae), mean_mse = mean(mse)) plotCV(example_sce, name = "cv_GMF") ``` If the dataset is large or you are unsure about the optimal range of components to test, an alternative is the scree plot approach. This method uses PCA on deviance residuals to estimate eigenvalues, providing a fast approximation of the optimal dimensionality. This can be done using `runRankGMF` or `calculateRankGMF` followed by `plotRank` or `screeplot_rank` respectively. Note that now, the `maxcomp` argument can be defined, which is the number of eigenvalues computed. ```{r} example_sce <- runRankGMF( example_sce, X = X, exprs_values="counts", family = poisson(), maxcomp = 10, ntop = 50) plotRank(example_sce, maxcomp = 10) ``` After choosing the number of components to use in the final dimensionality reduction, `runGMF` or `calculateGMF` can be used. Again, one should select the distribution family (`family`), the dimensionality (`ncomponents`), the known covariate matrices to account for (`X` and `Z`) and the assay used (`exprs_values` or `assay.type`). Unlike `runPCA`, `runGMF` uses all features by default. If you want to select the most variable genes instead, set `ntop`. The results are stored in the `reducedDim` slot of the `SingleCellExperiment` object. Additional information such as the `rotation` matrix, parameter estimates, the optimization history of `sgdGMF` framework and many more are available in the `attributes`. See `runGMF` for all outputs. ```{r} example_sce <- runGMF( example_sce, X = X, exprs_values="counts", family = poisson(), ncomponents = 3, # Use optimal dimensionality, here arbitrarily chosen as 3 ntop = 50, name = "GMF") ``` ```{r, results = 'hide'} reducedDimNames(example_sce) head(reducedDim(example_sce)) names(attributes(reducedDim(example_sce, type = "GMF"))) head(attr(reducedDim(example_sce, type = "GMF"), "rotation")) tail(attr(reducedDim(example_sce, type = "GMF"), "trace")) ``` To visualize the reduced dimensions, you can use `plotReducedDim` from the `scater` package, specifying "GMF" as the dimension reduction method. Alternatively, the `plotGMF` function provides a direct wrapper for this. ```{r} plotReducedDim(example_sce, dimred = "GMF", colour_by = "Mutation_Status") ``` ```{r} sessionInfo() ```