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.

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.

1.1 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.

1.2 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)

2 Basics

2.1 Install microSTASIS

microSTASIS is a R package available via the Bioconductor repository (formerly at CRAN) which can be installed by using the following commands in your R session:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
## Check that you have a valid Bioconductor installation

2.2 Citing microSTASIS

We hope that microSTASIS will be useful for your research. Please use the following information to cite the package and the overall approach. Thank you!

## Citation info
#> To cite microSTASIS in publications, please use:
#>   Sanchez-Sanchez P, Santonja FJ, Benitez-Paez A (2022). "Assessment of
#>   human microbiota stability across longitudinal samples using
#>   iteratively growing-partitioned clustering." _Briefings in
#>   Bioinformatics_, *23*(2). ISSN -2577, doi:10.1093/bib/bbac055
#>   <https://doi.org/10.1093/bib/bbac055>, bbac055,
#>   <https://doi.org/10.1093/bib/bbac055>.
#> A BibTeX entry for LaTeX users is
#>   @Article{,
#>     title = {Assessment of human microbiota stability across longitudinal 
#>           samples using iteratively growing-partitioned clustering},
#>     author = {Pedro Sanchez-Sanchez and Francisco J. Santonja and Alfonso Benitez-Paez},
#>     journal = {Briefings in Bioinformatics},
#>     volume = {23},
#>     number = {2},
#>     year = {2022},
#>     month = {02},
#>     issn = {-2577},
#>     doi = {10.1093/bib/bbac055},
#>     url = {https://doi.org/10.1093/bib/bbac055},
#>     note = {bbac055},
#>   }

3 Quick start to using microSTASIS

## Load the package

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
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
## 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.

3.1 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 = "_")

3.2 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 outside the scale range
#> (`stat_xsideboxplot()`).
#> Warning: Removed 16 rows containing missing values or values outside the scale range
#> (`geom_point()`).