--- title: "Introduction to microSTASIS" author: - name: Pedro Sánchez Sánchez email: bio.pedro.technology@gmail.com output: BiocStyle::html_document: self_contained: yes toc: true toc_float: true toc_depth: 2 code_folding: show date: "`r doc_date()`" package: "`r pkg_ver('microSTASIS')`" vignette: > %\VignetteIndexEntry{Introduction to microSTASIS} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", crop = NULL ## Related to https://stat.ethz.ch/pipermail/bioc-devel/2020-April/016656.html ) ``` ```{r vignetteSetup, echo=FALSE, warning = FALSE} ## Track time spent on making the vignette startTime <- Sys.time() ## Bib setup suppressPackageStartupMessages(library("RefManageR")) ## Write bibliography information bib <- c( R = citation(), BiocStyle = citation("BiocStyle")[1], knitr = citation("knitr")[1], RefManageR = citation("RefManageR")[1], rmarkdown = citation("rmarkdown")[1], sessioninfo = citation("sessioninfo")[1], testthat = citation("testthat")[1], microSTASIS = citation("microSTASIS")[1] ) ``` # Introduction **µSTASIS**, or microSTASIS, was firstly develop for assessing the stability of microbiota over time. The initial idea behind it was to warrant the proper processing of high-dimensional, sparse, undetermined, right skewed, overdispersed and compositional data. Specifically, it aims to fill the gap of temporal stability metrics in the compositional analysis of microbiome data. However, it can also work if the data are not transformed for belonging to the [Aitchison simplex](https://en.wikipedia.org/wiki/Compositional_data). On one side, the assumptions in which µSTASIS is based are two: the individual-specific composition and the within-individual variability over time. On the other side, the output, which we called mS score, is easy to interpret and provides a contextualized and intuitive metric to estimate temporal microbiota stability. Also, the package incorporates leave-p-out (or k) cross-validation routines, that compute the mean absolute error, and multiple functions to visualize the result. Therefore, sort of the questions that could be answered are related with robust time-resolved definitions of stability to assess microbiota behaviour against perturbations or to search for specific compositions that present a dynamic equilibrium around a central attractor state rather than overcoming a threshold of no return. ## Algorithm Firstly, two samples from the same individual has to be paired. For that, one can merge the sequential paired times (t1 with t2, t2 with t3...) or use a non-sequential order (t1 with t3). Therefore, from a single large data matrix, we would generate more that are smaller, that is, with fewer observations (samples). Once we have the paired samples, it is time for the main algorithm: iterative clustering. Concretely, Hartigan-Wong k-means algorithm is used as many times as possible for stressing out paired samples from the same individuals to test if they remain together for multiple numbers of clusters over a whole data set of individuals. Also, this is corrected in those cases where the distance between samples from the same individuals is lower than the distances to their respective centroids. ![](../man/figures/microSTASIS_algo.png) ## Why Bioconductor? Although the package was initially released at CRAN, we feel that releasing it at Bioconductor could add better experience to the user. Some functions were internally modified leveraging on other Bioconductor packages as well as we added interoperability of our workflow with `TreeSummarizedExperiment`. * `r Biocpkg("microSTASIS")` `r Citep(bib[["microSTASIS"]])` # Basics ## Install `microSTASIS` `r Biocpkg("microSTASIS")` is a `R` package available via the [Bioconductor](http://bioconductor.org) repository (formerly at [CRAN](https://cran.r-project.org/package=microSTASIS)) which can be installed by using the following commands in your `R` session: ```{r "install", eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install("microSTASIS") ## Check that you have a valid Bioconductor installation BiocManager::valid() ``` ## Citing `microSTASIS` We hope that `r Biocpkg("microSTASIS")` will be useful for your research. Please use the following information to cite the package and the overall approach. Thank you! ```{r "citation"} ## Citation info citation("microSTASIS") ``` # Quick start to using `microSTASIS` ```{r "start"} ## Load the package library("microSTASIS") ``` The input data can be a matrix where the rownames include ID, common pattern and sampling time and columns are the microbial features i.e. amplicon sequence variants (ASVs) or operational taxonomic units (OTUs). Alternatively, it can be a TreeSummarizedExperiment object containing the matrix in an `assay` slot or within an `altExps` (alternative experiment) slot. Then, the ID and sampling time info should be in the `colData` slot. ```{r "input"} ## Show the input data data(clr) clr[1:8, 1:6] ``` The first step is to subset an initial matrix of microbiome data with multiple sampling points. The`pairedTimes()` function already do it for every possible paired times in a `sequential = TRUE` way or for specific times points, for example: `sequential = FALSE, specifiedTimePoints = c("1", "3")`. The output is a list with length equal to the number of paired times. ```{r "pairedTimes"} ## Subseting the initial data to a list of multiple paired times times <- pairedTimes(data = clr, sequential = TRUE, common = "_0_") ``` Alternatively, the input can also be a `TreeSummarizedExperiment`, like below. ```{r "TSE"} ## Loading the TSE package suppressPackageStartupMessages(library(TreeSummarizedExperiment)) ``` ```{r "makeObject"} ## Creating a TreeSummarizedExperiment object ID_common_time <- rownames(clr) samples <- 1:dim(clr)[1] metadata <- data.frame(samples, stringr::str_split(ID_common_time, "_", simplify = TRUE)) colnames(metadata)[2:4] <- c("ID", "common", "time_point") both_tse <- TreeSummarizedExperiment(assays = list(counts = t(clr)), colData = metadata) altExp(both_tse) <- SummarizedExperiment(assays = SimpleList(data = t(clr)), colData = metadata) assayNames(altExp(both_tse)) <- "data" ``` ```{r "pairedTimesTSE"} ## Subseting the data as two different TSE objects to two list of multiple paired times TSE_times <- pairedTimes(data = both_tse, sequential = TRUE, assay = "counts", alternativeExp = NULL, ID = "ID", timePoints = "time_point") TSE_altExp_times <- pairedTimes(data = both_tse, sequential = TRUE, assay = "data", alternativeExp = "unnamed1", ID = "ID", timePoints = "time_point") ``` It is pretty important not to forget `sequential = FALSE` when the user wants to get `specifiedTimePoints`. ## Main algorithm Then, it is time for `iterativeClustering()`, which can be done in parallel and runs for every item in the list. ```{r "iterativeClustering"} ## Main algorithm applied to the matrix-derived object mS <- iterativeClustering(pairedTimes = times, common = "_") ``` Note that the downstream functions work equally once we have the `pairedTimes()`output. Also, `BPPARAM` can be specified as following `BiocParallel` guides. ```{r "iterativeClusteringTSE"} ## Main algorithm applied to the TSE-derived object mS_TSE <- iterativeClustering(pairedTimes = TSE_times, common = "_") ``` ## Visualization There are two main functions for visualizing the results: a scatter plot with boxplots on the side and a heatmap. ```{r "visualization"} ## Prepare the result for the later visualization functions results <- mSpreviz(results = mS, times = times) ## Scatter plot of the stability scores per patient and time plotmSscatter(results = results, order = "median", times = c("t1_t25", "t25_t26"), gridLines = TRUE, sideScale = 0.2) ## Heatmap of the stability scores per patient and time plotmSheatmap(results = results, order = "mean", times = c("t1_t25", "t25_t26"), label = TRUE) ``` Both functions can be sorted by the stability score by the argument `order`. The interpretation of both visualizations is similar. One the one side we can see how the stability scores differ a little bit between the sampling points, although the boxplots show that the values range from 0 to 1. Equally, the heatmap allows to see what are the most stable patients over sampling points. ## Cross-validation The stability results can be validated by cross validation (CV) in two ways: leave-one-out or leave-k-out. `k` is the number of individuals removed in each time that `iterativeClustering()` is run internally; it is equal to 1 for the LOOCV. ```{r "CV"} ## Cross validation in a parallelized way with different k values cv_klist_k2 <- BiocParallel::bpmapply(iterativeClusteringCV, name = names(times), k = rep(2L, 3), MoreArgs = list(pairedTimes = times, results = mS, common = "_0_")) ``` Then, the result can be displayed in two ways: as the mean absolute error or as a visual representation of how the mS score change. ```{r "error"} ## Calculate the mean absolute error after computing the cross validation MAE_t1_t25 <- mSerrorCV(pairedTime = times$t1_t25, CVklist = cv_klist_k2[[1]], k = 2L) ``` Remember that for both `mSerrorCV()` and `plotmSlinesCV()`, `k` argument has to be the same previously used for extracting the input for `CVklist`. For that, we recommend to save it in another local variable and then use it for the arguments of the different functions. ```{r "error viz"} ## Prepare the result for the later visualization functions MAE <- mSpreviz(results = list(MAE_t1_t25), times = list(t1_t25 = times$t1_t25)) ## Heatmap of the mean absolute error values per patient and time plotmSheatmap(results = MAE, times = c("t1_t25", "t25_t26"), label = TRUE, high = 'red2', low = 'forestgreen', midpoint = 5) ``` The visualizations allows to look at the reliability of the stability metric for our data set and to detect some cases where the score should not be taken doubtlessly. ```{r "lines"} ## Line plot of the stability scores per patient at a given paired times ## Both the proper score after iterative clustering and the cross validation ones plotmSlinesCV(pairedTime = times$t1_t25, CVklist = cv_klist_k2[[1]], k = 2L) ``` The lines visualizations should be seen as a way to track how the stability score for an individual in a given paired time changes in dependence to the data set i.e. when removing other individuals. Variability is expected, but heavy drops or increases not. ## Metadata integration The analysis can also be integrated with metadata of the individuals in the following way. ```{r "metadata"} ## Create a metadata data frame metadata <- data.frame(Sample = rownames(clr), age = c(rep("youth", 65), rep("old", 131-65))) ## Modify the metadata to the proper format that match with results group <- mSmetadataGroups(metadata = metadata, samples = metadata$Sample, common = "_0_", individuals = results$individual, variable = "age") ## Boxplot with individual stability scores split by groups plotmSdynamics(results, groups = group, points = TRUE, linetype = 0) ``` This plot allow to track general trends of stability in different groups across time. # Reproducibility The `r Biocpkg("microSTASIS")` package `r Citep(bib[["microSTASIS"]])` was made possible thanks to: * R `r Citep(bib[["R"]])` * `r Biocpkg("BiocStyle")` `r Citep(bib[["BiocStyle"]])` * `r CRANpkg("knitr")` `r Citep(bib[["knitr"]])` * `r CRANpkg("RefManageR")` `r Citep(bib[["RefManageR"]])` * `r CRANpkg("rmarkdown")` `r Citep(bib[["rmarkdown"]])` * `r CRANpkg("sessioninfo")` `r Citep(bib[["sessioninfo"]])` * `r CRANpkg("testthat")` `r Citep(bib[["testthat"]])` This package was developed using `r BiocStyle::Biocpkg("biocthis")`. Code for creating the vignette. ```{r createVignette, eval=FALSE} ## Create the vignette library("rmarkdown") system.time(render("microSTASIS.Rmd", "BiocStyle::html_document")) ## Extract the R code library("knitr") knit("microSTASIS.Rmd", tangle = TRUE) ``` Date the vignette was generated. ```{r reproduce1, echo=FALSE} ## Date the vignette was generated Sys.time() ``` Wallclock time spent generating the vignette. ```{r reproduce2, echo=FALSE} ## Processing time in seconds totalTime <- diff(c(startTime, Sys.time())) round(totalTime, digits = 3) ``` # R session information ```{r reproduce3, echo=FALSE} ## Session info library("sessioninfo") options(width = 120) session_info() ``` # Bibliography This vignette was generated using `r Biocpkg("BiocStyle")` `r Citep(bib[["BiocStyle"]])` with `r CRANpkg("knitr")` `r Citep(bib[["knitr"]])` and `r CRANpkg("rmarkdown")` `r Citep(bib[["rmarkdown"]])` running behind the scenes. Citations made with `r CRANpkg("RefManageR")` `r Citep(bib[["RefManageR"]])`. ```{r vignetteBiblio, results = "asis", echo = FALSE, warning = FALSE, message = FALSE} ## Print bibliography PrintBibliography(bib, .opts = list(hyperlink = "to.doc", style = "html")) ```