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.
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.
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.
- microSTASIS (Sanchez-Sanchez, Santonja, and Benitez-Paez, 2022)
Quick start to using microSTASIS
## 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.
## Show the input data
data(clr)
clr[1:8, 1:6]
#> 1 2 3 4 6 7
#> 003_0_1 2.891583 2.996944 -1.1114255 -1.1114255 -1.8606574 1.3875060
#> 003_0_25 7.430753 7.365052 4.4745478 4.7213331 0.9239730 -0.8677865
#> 003_0_26 6.359208 6.314002 1.3350641 1.7537744 0.1956298 -0.2743738
#> 003_0_27 7.364809 7.372548 4.5086757 4.6234512 -2.7098118 -2.7098118
#> 004_0_1 4.033265 4.028176 -2.3723362 -2.3723362 1.8855559 -2.3723362
#> 004_0_25 4.666157 4.539863 -1.2180974 -1.2180974 -2.0048403 0.5886198
#> 004_0_26 5.005840 5.082213 -0.2471546 -0.2471546 7.2118445 -1.5034582
#> 004_0_27 4.492273 4.449713 3.5114435 -0.9009753 8.4833123 -0.9009753
The first step is to subset an initial matrix of microbiome data with multiple sampling points. ThepairedTimes() 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.
## 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.
## Loading the TSE package
suppressPackageStartupMessages(library(TreeSummarizedExperiment))
## 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"
## 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.
## 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.
## 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.
## 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)
#> Registered S3 method overwritten by 'ggside':
#> method from
#> +.gg ggplot2
#> Warning: Removed 16 rows containing non-finite values (`stat_boxplot()`).
#> Warning: Removed 16 rows containing missing values (`geom_point()`).
## Heatmap of the stability scores per patient and time
plotmSheatmap(results = results, order = "mean", times = c("t1_t25", "t25_t26"),
label = TRUE)
#> Warning: Removed 16 rows containing missing values (`geom_text()`).
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.
## 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.
## 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.
## 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.
## 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)
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
#> returning NA
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.
Bibliography
This vignette was generated using BiocStyle (Oleś, 2023)
with knitr (Xie, 2023) and rmarkdown (Allaire, Xie, Dervieux et al., 2023) running behind the scenes.
Citations made with RefManageR (McLean, 2017).
[1]
J. Allaire, Y. Xie, C. Dervieux, et al.
rmarkdown: Dynamic Documents for R.
R package version 2.25.
2023.
URL: https://github.com/rstudio/rmarkdown.
[2]
M. W. McLean.
“RefManageR: Import and Manage BibTeX and BibLaTeX References in R”.
In: The Journal of Open Source Software (2017).
DOI: 10.21105/joss.00338.
[3]
A. Oleś.
BiocStyle: Standard styles for vignettes and other Bioconductor documents.
R package version 2.30.0.
2023.
DOI: 10.18129/B9.bioc.BiocStyle.
URL: https://bioconductor.org/packages/BiocStyle.
[4]
R Core Team.
R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing.
Vienna, Austria, 2023.
URL: https://www.R-project.org/.
[5]
P. Sanchez-Sanchez, F. J. Santonja, and A. Benitez-Paez.
“Assessment of human microbiota stability across longitudinal
samples using iteratively growing-partitioned clustering”.
In: Briefings in Bioinformatics 23.2 (Feb. 2022). bbac055.
ISSN: -2577.
DOI: 10.1093/bib/bbac055.
URL: https://doi.org/10.1093/bib/bbac055.
[6]
H. Wickham.
“testthat: Get Started with Testing”.
In: The R Journal 3 (2011), pp. 5–10.
URL: https://journal.r-project.org/archive/2011-1/RJournal_2011-1_Wickham.pdf.
[7]
H. Wickham, W. Chang, R. Flight, et al.
sessioninfo: R Session Information.
R package version 1.2.2.
2021.
URL: https://CRAN.R-project.org/package=sessioninfo.
[8]
Y. Xie.
knitr: A General-Purpose Package for Dynamic Report Generation in R.
R package version 1.44.
2023.
URL: https://yihui.org/knitr/.
