LegATo Tutorial
Aubrey R. Odom
Program in Bioinformatics, Boston University, Boston, MAaodom@bu.edu
January 30, 2024
LegATo_vignette.RmdIntroduction to LegATo
Streamlining longitudinal microbiome profiling in Bioconductor
Microbiome time-series data presents several distinct challenges, including complex covariate dependencies and a variety of longitudinal study designs. Furthermore, while several individual tools have been created to aid in longitudinal microbiome data analysis, an all-inclusive R-based toolkit has remained elusive.
To address these concerns, LegATo (Longitudinal mEtaGenomic Analysis TOolkit) has risen from the ashes. LegATo is a suite of open-source software tools for longitudinal microbiome analysis that is extendable to several different study forms, all while promising optimal ease-of-use for researchers. This toolkit will allow researchers to determine which microbial taxa are affected over time by perturbations such as onset of disease or lifestyle choices, and to predict the effects of these perturbations over time, including changes in composition or stability of commensal bacteria.
Currently, LegATo entertains a number of data cleaning, aggregation, modeling and testing procedures. As we develop or learn of new methods, we will add to the toolkit accordingly. We will soon add hierarchical clustering tools and multivariate generalized estimating equations (JGEEs) to adjust for the compositional nature of microbiome data.
Getting started with LegATo
Compatibility with MultiAssayExperiment and
SummarizedExperiment objects
The convenience of LegATo is in part credited to the integration of
MultiAssayExperiment and SummarizedExperiment
objects. These are reliable and clean data structures developed by the
as part of the MultiAssayExperiment
and SummarizedExperiment
packages.
The MultiAssayExperiment container concisely stores a
SummarizedExperiment object, which aggregates data matrices
along with annotation information, metadata, and reduced dimensionality
data (PCA, t-SNE, etc.). To learn more about proper usage and context of
these objects, you may want to take a look at the MultiAssayExperiment
package vignette and SummarizedExperiment
package vignette.
To install these two packages, run the following code:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SummarizedExperiment")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SummarizedExperiment")Installation and setup
In order to install the development version of LegATo from Github, run the following code:
if (!requireNamespace("devtools", quietly=TRUE))
install.packages("devtools")
devtools::install_github("wejlab/LegATo")Now we’ll load the LegATo package into our library:
library(LegATo)Formatting input data for use with LegATo
LegATo functions expect data that is properly formatted in a
MultiAssayExperiment container that is compatible with the
animalcules
package. Obtaining properly formatted data can be achieved easily
via the create_formatted_MAE() function from the LegATo
package, as detailed below.
A note on compatibility with other packages
If you’ve analyzed your metagenomic or 16S data via the MetaScope package,
you can aggregate your sample outputs with user-input metadata using the
convert_animalcules() or
convert_animalcules_silva() functions. These functions
create a taxonomy table and collate all data into a LegATo- and
animalcules-compatible MultiAssayExperiment object.
You can also format your data analyzed elsewhere (like PathoScope
2.0) with the animalcules
package in the upload step with the R Shiny app, and then select
“Download Animalcules File” to obtain a MAE object that can
be used with LegATo.
Requirements for using create_formatted_MAE()
For the purposes of this example, we’ll look at some example inputs unrelated to the main HIV-E infant data:
counts <- system.file("extdata", "counts.csv", package = "LegATo") |>
read.csv(row.names = 1) |>
dplyr::rename_with(function(x) stringr::str_replace(x, "\\.", "-"))
tax <- system.file("extdata", "tax.csv", package = "LegATo") |> read.csv(row.names = 1)
sample <- system.file("extdata", "sample.csv", package = "LegATo") |> read.csv(row.names = 1)Now we’ll look at the formatting of each:
ndim <- 5
counts[seq_len(ndim), seq_len(ndim)] |>
knitr::kable(caption = "Counts Table Preview",
label = NA)| X-1 | X-2 | X-3 | X-4 | X-5 | |
|---|---|---|---|---|---|
| Acinetobacter_beijerinckii | 430 | 18103 | 810 | 103 | 74 |
| Acinetobacter_bouvetii | 0 | 0 | 0 | 0 | 0 |
| Acinetobacter_guillouiae | 0 | 5 | 0 | 0 | 0 |
| Acinetobacter_gyllenbergii | 6 | 1 | 1 | 1 | 11 |
| Acinetobacter_indicus | 0 | 0 | 9 | 0 | 0 |
| superkingdom | phylum | class | order | family | genus | species | |
|---|---|---|---|---|---|---|---|
| Acinetobacter_beijerinckii | Bacteria | Proteobacteria | Gammaproteobacteria | Moraxellales | Moraxellaceae | Acinetobacter | Acinetobacter_beijerinckii |
| Acinetobacter_bouvetii | Bacteria | Proteobacteria | Gammaproteobacteria | Moraxellales | Moraxellaceae | Acinetobacter | Acinetobacter_bouvetii |
| Acinetobacter_guillouiae | Bacteria | Proteobacteria | Gammaproteobacteria | Moraxellales | Moraxellaceae | Acinetobacter | Acinetobacter_guillouiae |
| Acinetobacter_gyllenbergii | Bacteria | Proteobacteria | Gammaproteobacteria | Moraxellales | Moraxellaceae | Acinetobacter | Acinetobacter_gyllenbergii |
| Acinetobacter_indicus | Bacteria | Proteobacteria | Gammaproteobacteria | Moraxellales | Moraxellaceae | Acinetobacter | Acinetobacter_indicus |
| Sample | Subject | Sex | Month | Group | Pairing | HairLength | Age | |
|---|---|---|---|---|---|---|---|---|
| X-1 | X-1 | S1 | Male | 1 | A | 1 | 48.62753 | 21.9 |
| X-2 | X-2 | S2 | Male | 1 | B | 1 | 53.53195 | 67.7 |
| X-3 | X-3 | S3 | Female | 1 | A | 2 | 47.61090 | 78.9 |
| X-4 | X-4 | S4 | Female | 1 | B | 2 | 49.72870 | 48.6 |
| X-5 | X-5 | S5 | Female | 1 | A | 3 | 59.04509 | 48.4 |
You should note the following:
- The
countstable is composed of raw microbial counts. There is one row for each taxon, and one column for each sample. - The
taxtable is composed of taxonomy lineages. There is one row for each taxon, and any number of columns each taxonomic level (e.g., superkingdom, genus, species). The lowest level of taxonomy serves as the rownames for the table, and is also stored in its own column. Taxonomic levels increase in granularity from left to right. - The
sampletable contains entries for each sample, with unique sample IDs on the rows. The columns are metadata for the samples, such as subjects or units on which repeated measures were taken, groupings, pairings, and covariates. rownames(counts) == rownames(tax)colnames(counts) == colnames(sample)
Use create_formatted_MAE()
So once your data are formatted correctly, you can pretty easily use
create_formatted_MAE() like so:
output <- create_formatted_MAE(counts_dat = counts,
tax_dat = tax,
metadata_dat = sample)
class(output)## [1] "MultiAssayExperiment"
## attr(,"package")
## [1] "MultiAssayExperiment"
MultiAssayExperiment::assays(output)## List of length 1
## names(1): MicrobeGenetics
SummarizedExperiment::assays(output[["MicrobeGenetics"]])## List of length 1
## names(1): MGXAdding information to your metadata later in your analysis
If information needs to be added to your data object at some point in
the analysis, it is easiest to manipulate the raw data objects
(potentially via parse_MAE_SE()) and then recreate the
MAE object with create_formatted_MAE().
Example data
To illustrate the capabilities of LegATo, we will turn to a published dataset from the following paper:
Odom-Mabey AR, Gill CJ, Pieciak R et al. Characterization of longitudinal nasopharyngeal microbiome patterns in maternally HIV-exposed Zambian infants [version 1; peer review: 1 approved with reservations]. Gates Open Res 2022, 6:143 (https://doi.org/10.12688/gatesopenres.14041.1)
The raw dataset is archived in Zenodo:
Zenodo: Underlying data for ‘Characterization of longitudinal nasopharyngeal microbiome patterns in maternally HIV-exposed Zambian infants’. https://doi.org/10.5281/zenodo.725531324
Further details on how the dataset was altered for inclusion in this package are provided here.
Example data context and structure
The example data consists of 167 NP swabs of healthy HIV-exposed, uninfected (HEU; n=10) infants and their HIV(+) mothers and HIV-unexposed, uninfected (HUU; n=10) infants and their HIV(-) mothers. A total of 7 samples were collected per infant, with some missingness in the data. These swabs were identified from a sample library collected in Lusaka, Zambia between 2015 and 2016.
The analysis objective is to parse the association between the NP resident bacteria and infant HIV exposure during the first 3.5 months (14 weeks) of life, a critical time in microbiome maturation.
All data was processed with PathoScope 2.0, and the
sample outputs were aggregated with the animalcules R
package.
For this analysis, we will work with the infant data only and subset the data accordingly.
dat <- system.file("extdata", "MAE.RDS", package = "LegATo") |>
readRDS()
dat_subsetted <- MultiAssayExperiment::subsetByColData(dat,
dat$MothChild == "Infant")This leaves us with 129 samples in our analysis.
Initial data exploration
MAE manipulation
LegATo has several functions for manipulating
MAE objects, delineated below:
Clean up animalcules-formatted MAEs
If your data were aggregated via the animalcules
package, as was done for the example data, the samples were aggregated
at the strain taxonomic level and have residual taxonomy IDs present.
The clean_animalcules_MAE() function cleans up data
suffering from these issues.
dat_cleaned <- clean_animalcules_MAE(dat_subsetted)## Registered S3 method overwritten by 'httr':
## method from
## print.response rmutilFilter MAE
Many metagenomic pipelines identify taxon abundances at extremely
small levels, which can be noisy to deal with in an analysis. The
filter_animalcules_MAE function smoothly transforms reads
belonging to taxa with an overall genera threshold under the
filter_prop (filter proportion) argument, which we will set
as 0.005.
dat_filt <- filter_animalcules_MAE(dat_cleaned, 0.05)Parse MAE to extract data
If you want to take a closer look at your data, you can easily
extract it into the counts, taxonomy, and sample metadata tables using
the parse_MAE_SE() function.
parsed <- parse_MAE_SE(dat_filt, which_assay = "MicrobeGenetics", type = "MAE")
parsed$counts[seq_len(5), seq_len(5)] |>
knitr::kable(caption = "Counts Table")| X1755 | X2216 | X2431 | X2561 | X2699 | |
|---|---|---|---|---|---|
| Corynebacterium accolens | 0 | 9643 | 0 | 0 | 43435 |
| Corynebacterium ammoniagenes | 0 | 0 | 0 | 0 | 0 |
| Corynebacterium appendicis | 0 | 0 | 0 | 0 | 0 |
| Corynebacterium argentoratense | 0 | 0 | 0 | 4 | 0 |
| Corynebacterium aurimucosum | 0 | 0 | 0 | 58 | 0 |
| Sample | Subject | HIVStatus | MothChild | timepoint | Age | pairing | |
|---|---|---|---|---|---|---|---|
| X1755 | X1755 | 0469-1 | Control | Infant | 0 | 7 | 1 |
| X2216 | X2216 | 0507-1 | Control | Infant | 0 | 6 | 2 |
| X2431 | X2431 | 0539-1 | Control | Infant | 0 | 9 | 3 |
| X2561 | X2561 | 0554-1 | HIV | Infant | 0 | 6 | 1 |
| X2699 | X2699 | 0620-1 | HIV | Infant | 0 | 4 | 2 |
| superkingdom | phylum | class | order | family | genus | species | |
|---|---|---|---|---|---|---|---|
| Corynebacterium accolens | Bacteria | Actinobacteria | Actinomycetia | Corynebacteriales | Corynebacteriaceae | Corynebacterium | Corynebacterium accolens |
| Corynebacterium ammoniagenes | Bacteria | Actinobacteria | Actinomycetia | Corynebacteriales | Corynebacteriaceae | Corynebacterium | Corynebacterium ammoniagenes |
| Corynebacterium appendicis | Bacteria | Actinobacteria | Actinomycetia | Corynebacteriales | Corynebacteriaceae | Corynebacterium | Corynebacterium appendicis |
| Corynebacterium argentoratense | Bacteria | Actinobacteria | Actinomycetia | Corynebacteriales | Corynebacteriaceae | Corynebacterium | Corynebacterium argentoratense |
| Corynebacterium aurimucosum | Bacteria | Actinobacteria | Actinomycetia | Corynebacteriales | Corynebacteriaceae | Corynebacterium | Corynebacterium aurimucosum |
Summarizing microbial counts data
It is fairly convenient to summarize the average number of reads for
a report with get_summary_table(). The table groups by
user-provided discrete covariates.
group_vars <- c("HIVStatus", "MothChild")
get_summary_table(dat_filt, group_vars) |>
knitr::kable(caption = "Summary Table", label = NA)| HIVStatus | MothChild | mean_reads | sd_reads | min_reads | max_reads | num_total |
|---|---|---|---|---|---|---|
| Control | Infant | 156209.2 | 408000.52 | 24176 | 3016276 | 68 |
| HIV | Infant | 90008.2 | 51395.08 | 13218 | 202079 | 61 |
The get_top_taxa function outputs a
data.frame that lists taxa in order of relative
abundance.
best_genus <- get_top_taxa(dat_filt, "genus")
best_genus |> knitr::kable(caption = "Table of genera, ranked by abundance")| taxon | allmeans |
|---|---|
| Dolosigranulum | 0.1892754 |
| Streptococcus | 0.1730879 |
| Moraxella | 0.1713469 |
| Staphylococcus | 0.1554891 |
| Corynebacterium | 0.1290279 |
| Other | 0.0979779 |
| Haemophilus | 0.0837949 |
Other data manipulation functions
If you want to conduct your own analyses, the
get_long_data() function will prove convenient. The
get_stacked_data() function can be used for certain
visualizations that utilize a relative abundance aggregation
approach.
longdat <- get_long_data(dat_filt, "genus", log = TRUE, counts_to_CPM = TRUE)
stackeddat <- get_stacked_data(dat_filt, "genus", covariate_1 = "HIVStatus",
covariate_time = "timepoint")Visualizing Data
There are several plots by which we can visualize changes in relative abundance over time, accounting for a given covariate. In this case, we are interested in HIV exposure.
We’ll select a palette using paletteer.
this_palette <- c("#FED439", "#709AE1", "#8A9197", "#D2AF81", "#FD7446", "#D5E4A2", "#197EC0", "#F05C3B", "#46732E",
"#71D0F5", "#370335", "#075149", "#C80813", "#91331F", "#1A9993", "#FD8CC1")|>
magrittr::extract(seq_len(nrow(best_genus)))Alluvial plot
Alluvial diagrams illustrate individual taxa as stream fields that
change position at different time points. The height of a stream field
represents the relative abundance of that taxon. At a given time point,
stream fields are ranked from the highest to lowest abundance (top to
bottom). These can be constructed with plot_alluvial().
plot_alluvial(dat = dat_filt,
taxon_level = "genus",
covariate_1 = "HIVStatus",
covariate_time = "timepoint",
palette_input = this_palette,
subtitle = "Alluvial plot")
Spaghetti plot
We can create spaghetti or volatility plots to elucidate changes over
time on a sample level for a given taxon. This is advantageous as other
visualization methods are often aggregates of multiple samples and lack
granularity. These plots can be created with
plot_spaghetti().
p <- plot_spaghetti(dat = dat_filt,
covariate_time = "timepoint",
covariate_1 = "HIVStatus",
unit_var = "Subject",
taxon_level = "genus",
which_taxon = "Staphylococcus",
palette_input= this_palette,
title = "Spaghetti Plot",
subtitle = NULL) +
ggplot2::xlab("Infant Age (Days)") +
ggplot2::ylab("Relative Abundance (log CPM)")Stacked bar plot
Stacked bar plots are used here to visualize the relative abundance
of microbes at a given taxonomic level in each sample, represented as a
single bar, labeled by time point, and plotted within each HIV exposure
status group for separate mothers and infant comparisons. Use
plot_stacked_bar.
plot_stacked_bar(dat_filt, "genus",
"HIVStatus",
"timepoint",
palette_input = this_palette)
Stacked area chart
Stacked area charts are similar to stacked bar plots, but provide for
continuity between time points. These are created with
plot_stacked_area().
plot_stacked_area(dat_filt, "genus",
"HIVStatus",
"timepoint",
palette_input = this_palette)
Heatmap
With plot_heatmap, you can plot a heatmap of a specific
microbe to determine changes along one or more covariates.
this_taxon <- parsed$counts |>
animalcules::upsample_counts(parsed$tax, "genus") |>
animalcules::counts_to_logcpm()
p <- plot_heatmap(inputData = this_taxon,
annotationData = dplyr::select(parsed$sam, "timepoint", "HIVStatus", "pairing"),
name = "Data",
plot_title = "Example",
plottingColNames = NULL,
annotationColNames = NULL,
colList = list(),
scale = FALSE,
showColumnNames = FALSE,
showRowNames = FALSE,
colorSets = c("Set1", "Set2", "Set3", "Pastel1", "Pastel2", "Accent", "Dark2",
"Paired"),
choose_color = c("blue", "gray95", "red"),
split_heatmap = "none",
column_order = NULL
)
Longitudinal data analysis
For a concrete analysis of longitudinal microbiome data,
LegATo provides NMIT, both a paired and unpaired
multivariate Hotelling’s T-squared test, generalized estimating
equations (GEEs) and linear mixed effects models (LMEMs/LMMs).
Detailed explanations of each are available on the LegATo website:
- Nonparametric microbial interdependence test (NMIT)
- Modeling with GEEs and LMMs
- Hotelling’s T-Squared test
Otherwise a brief overview of each method is provided below.
Nonparametric microbial interdependence test (NMIT)
NMIT is a multivariate distance-based test intended to evaluate the association between key phenotypic variables and microbial interdependence. The test determines longitudinal sample similarity as a function of temporal microbial composition.
The authors thank Yilong Zhang for providing his code for adaptation into LegATo.
Citations:
Yilong Zhang, Sung Won Han, Laura M Cox, and Huilin Li. A multivariate distance-based analytic framework for microbial interdependence association test in longitudinal study. Genetic epidemiology, 41(8):769–778, 2017. doi:10.1002/gepi.22065.
dat_0.1 <- filter_animalcules_MAE(dat_cleaned, 0.1)
NMIT(dat_0.1, unit_var = "Subject", fixed_cov = "HIVStatus",
covariate_time = "timepoint",
method = "kendall", dist_type = "F",
heatmap = TRUE, classify = FALSE, fill_na = 0)
Hotelling’s T^2 Test
Hotelling’s T2 tests can be used to determine whether the microbiome profiles exhibit notable differences or trends across time and groups. We may then use t-tests to identify which genera contributed most to these differences.
For this example, HEU and HUU infants are designated as the two sampling units on which the relative abundances of the p most abundant genera will be measured.
For paired tests, we chose p = 6 variables to ensure that n < p so that singularity could be avoided and T2 could be properly computed, where n is the number of measurements in a sampling unit. Normality is met by using microbe abundances in log CPM units, which is calculated within the function.
test_hotelling_t2(dat = dat_filt,
test_index = which(dat_filt$MothChild == "Infant" &
dat_filt$timepoint == 6),
taxon_level = "genus",
# To avoid n < p, use top 5-6 species
num_taxa = 6,
paired = TRUE,
grouping_var = "HIVStatus",
pairing_var = "pairing")## $df1
## [1] 6
##
## $df2
## [1] 2
##
## $crit_F
## [1] 19.32953
##
## $F_stat
## [1] 2.690561
##
## $pvalue
## [1] 0.2955849Group comparisons can also be tested on unpaired data:
test_hotelling_t2(dat = dat_filt,
test_index = which(dat_filt$timepoint == 0),
taxon_level = "genus",
# To avoid n < p, use top 5-6 species
num_taxa = 6,
grouping_var = "HIVStatus",
unit_var = "Subject",
paired = FALSE)## $df1
## [1] 6
##
## $df2
## [1] 13
##
## $crit_F
## [1] 2.915269
##
## $F_stat
## [1] 0.8499803
##
## $pvalue
## [1] 0.5545602Modeling
Generalized Estimating Equations (GEEs)
Generalized estimating equations (GEEs) as described in Liang and Zeger (1986) and extended by Agresti (2002) have been widely used for modeling longitudinal data, and more recently for longitudinal microbiome data.
For each genus present in the microbial aggregate of samples, we model normalized log CPM relative taxon counts, estimating the effects of time point and HIV exposure status and their interaction, while accounting for the underlying structure of clusters formed by individual subjects.
output <- run_gee_model(dat_filt, unit_var = "Subject",
fixed_cov = c("HIVStatus", "timepoint"),
corstr = "ar1",
plot_out = FALSE,
plotsave_loc = ".",
plot_terms = NULL)
head(output) |> knitr::kable(caption = "GEE Outputs")| Coefficient | Coefficient Estimate | Standard Error | Statistic | Unadj p-value | Lower 95% CI | Upper 95% CI | Taxon | Adj p-value |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | 54987.561 | 45931.419 | 1.433208 | 0.2312416 | -35036.366 | 145011.49 | Dolosigranulum | 0.2697819 |
| (Intercept) | 4623.167 | 3270.381 | 1.998401 | 0.1574653 | -1786.661 | 11033.00 | Streptococcus | 0.2204514 |
| (Intercept) | 10838.983 | 4950.443 | 4.793900 | 0.0285607 | 1136.293 | 20541.67 | Moraxella | 0.0499812 |
| (Intercept) | 59424.348 | 26172.622 | 5.155068 | 0.0231786 | 8126.953 | 110721.74 | Staphylococcus | 0.0499812 |
| (Intercept) | 46510.999 | 20784.082 | 5.007831 | 0.0252329 | 5774.948 | 87247.05 | Corynebacterium | 0.0499812 |
| (Intercept) | 13394.894 | 4693.161 | 8.146062 | 0.0043155 | 4196.468 | 22593.32 | Other | 0.0302084 |
You can also create plots of the covariates, which will be saved to a folder specified by the user:
tempfolder <- tempfile()
dir.create(tempfolder)
# Trying out plotting
output <- run_gee_model(dat_filt, unit_var = "Subject",
taxon_level = "genus",
fixed_cov = c("HIVStatus", "timepoint"),
corstr = "ar1",
plot_out = TRUE,
plotsave_loc = tempfolder,
plot_terms = "timepoint",
width = 6, height = 4, units = "in", scale = 0.7)
unlink(tempfolder, recursive = TRUE)
Linear Mixed Models
Similarly, you can also run linear mixed-effects models:
output <- run_lmm_model(dat_filt, unit_var = "Subject",
taxon_level = "genus",
fixed_cov = c("timepoint", "HIVStatus"),
plot_out = FALSE,
plotsave_loc = ".",
plot_terms = NULL)## [1] "Dolosigranulum"
## [1] "Streptococcus"
## [1] "Moraxella"
## [1] "Staphylococcus"
## [1] "Corynebacterium"
## [1] "Other"
## [1] "Haemophilus"| effect | Coefficient | Coefficient Estimate | Standard Error | Statistic | df | Unadj p-value | Taxon | Adj p-value |
|---|---|---|---|---|---|---|---|---|
| fixed | (Intercept) | 54874.112 | 36388.018 | 1.508027 | 60.58909 | 0.1367446 | Dolosigranulum | 0.1914424 |
| fixed | (Intercept) | 5145.285 | 4970.582 | 1.035147 | 38.91848 | 0.3069896 | Streptococcus | 0.3581545 |
| fixed | (Intercept) | 10840.225 | 5796.478 | 1.870140 | 25.84537 | 0.0728319 | Moraxella | 0.1274558 |
| fixed | (Intercept) | 60627.051 | 16682.279 | 3.634219 | 44.96890 | 0.0007132 | Staphylococcus | 0.0049921 |
| fixed | (Intercept) | 45500.592 | 14615.001 | 3.113280 | 62.16246 | 0.0027961 | Corynebacterium | 0.0097863 |
| fixed | (Intercept) | 12702.733 | 5455.653 | 2.328362 | 23.11382 | 0.0289964 | Other | 0.0676583 |
Session Info
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: AlmaLinux 8.9 (Midnight Oncilla)
##
## Matrix products: default
## BLAS: /share/pkg.8/r/4.3.1/install/lib64/R/lib/libRblas.so
## LAPACK: /share/pkg.8/r/4.3.1/install/lib64/R/lib/libRlapack.so; LAPACK version 3.11.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] LegATo_0.99.0 BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.1 bitops_1.0-7
## [3] tibble_3.2.1 XML_3.99-0.16.1
## [5] rpart_4.1.19 rex_1.2.1
## [7] lifecycle_1.0.4 edgeR_3.42.4
## [9] doParallel_1.0.17 globals_0.16.2
## [11] lattice_0.21-8 MASS_7.3-60
## [13] insight_0.19.7 MultiAssayExperiment_1.28.0
## [15] backports_1.4.1 magrittr_2.0.3
## [17] limma_3.56.2 plotly_4.10.2
## [19] sass_0.4.8 rmarkdown_2.25
## [21] jquerylib_0.1.4 yaml_2.3.8
## [23] askpass_1.2.0 reticulate_1.30
## [25] minqa_1.2.6 RColorBrewer_1.1-3
## [27] multcomp_1.4-25 abind_1.4-5
## [29] zlibbioc_1.48.0 GenomicRanges_1.54.1
## [31] purrr_1.0.2 BiocGenerics_0.48.1
## [33] RCurl_1.98-1.14 TH.data_1.1-2
## [35] sandwich_3.1-0 circlize_0.4.15
## [37] GenomeInfoDbData_1.2.11 IRanges_2.36.0
## [39] S4Vectors_0.40.2 ggrepel_0.9.3
## [41] rmutil_1.1.10 inline_0.3.19
## [43] listenv_0.9.0 rentrez_1.2.3
## [45] vegan_2.6-4 umap_0.2.10.0
## [47] RSpectra_0.16-1 spatial_7.3-16
## [49] parallelly_1.36.0 pkgdown_2.0.7
## [51] permute_0.9-7 codetools_0.2-19
## [53] DelayedArray_0.28.0 DT_0.28
## [55] tidyselect_1.2.0 shape_1.4.6
## [57] ggeffects_1.3.4 farver_2.1.1
## [59] broom.mixed_0.2.9.4 lme4_1.1-35.1
## [61] stable_1.1.6 matrixStats_1.2.0
## [63] stats4_4.3.1 jsonlite_1.8.8
## [65] GetoptLong_1.0.5 ellipsis_0.3.2
## [67] survival_3.5-5 ggalluvial_0.12.5
## [69] iterators_1.0.14 emmeans_1.9.0
## [71] systemfonts_1.0.4 foreach_1.5.2
## [73] tools_4.3.1 ragg_1.2.5
## [75] Rcpp_1.0.12 glue_1.7.0
## [77] SparseArray_1.2.3 BiocBaseUtils_1.4.0
## [79] xfun_0.41 mgcv_1.8-42
## [81] DESeq2_1.40.1 MatrixGenerics_1.14.0
## [83] GenomeInfoDb_1.38.5 dplyr_1.1.4
## [85] numDeriv_2016.8-1.1 withr_3.0.0
## [87] BiocManager_1.30.22 timeSeries_4030.106.9000
## [89] fastmap_1.1.1 boot_1.3-28.1
## [91] fansi_1.0.6 shinyjs_2.1.0
## [93] openssl_2.1.1 animalcules_1.18.2
## [95] digest_0.6.34 R6_2.5.1
## [97] estimability_1.4.1 textshaping_0.3.6
## [99] colorspace_2.1-0 Cairo_1.6-2
## [101] modeest_2.4.0 utf8_1.2.4
## [103] tidyr_1.3.1 generics_0.1.3
## [105] data.table_1.14.10 httr_1.4.7
## [107] htmlwidgets_1.6.4 S4Arrays_1.2.0
## [109] pkgconfig_2.0.3 gtable_0.3.4
## [111] tsne_0.1-3.1 timeDate_4032.109
## [113] ComplexHeatmap_2.16.0 covr_3.6.2
## [115] XVector_0.42.0 furrr_0.3.1
## [117] htmltools_0.5.7 bookdown_0.34
## [119] clue_0.3-65 scales_1.3.0
## [121] Biobase_2.62.0 png_0.1-8
## [123] knitr_1.45 rstudioapi_0.15.0
## [125] reshape2_1.4.4 rjson_0.2.21
## [127] nloptr_2.0.3 statip_0.2.3
## [129] coda_0.19-4 nlme_3.1-162
## [131] zoo_1.8-12 cachem_1.0.8
## [133] GlobalOptions_0.1.2 stringr_1.5.1
## [135] parallel_4.3.1 fBasics_4022.94.9000
## [137] desc_1.4.3 pillar_1.9.0
## [139] grid_4.3.1 vctrs_0.6.5
## [141] xtable_1.8-4 cluster_2.1.4
## [143] evaluate_0.23 magick_2.7.4
## [145] mvtnorm_1.2-4 cli_3.6.2
## [147] locfit_1.5-9.8 compiler_4.3.1
## [149] rlang_1.1.3 crayon_1.5.2
## [151] labeling_0.4.3 forcats_1.0.0
## [153] plyr_1.8.9 fs_1.6.3
## [155] stringi_1.8.3 viridisLite_0.4.2
## [157] BiocParallel_1.36.0 TBSignatureProfiler_1.15.0
## [159] assertthat_0.2.1 lmerTest_3.1-3
## [161] munsell_0.5.0 lazyeval_0.2.2
## [163] geepack_1.3.9 Matrix_1.6-5
## [165] hms_1.1.3 stabledist_0.7-1
## [167] future_1.32.0 ggplot2_3.4.4
## [169] statmod_1.5.0 haven_2.5.4
## [171] SummarizedExperiment_1.32.0 highr_0.10
## [173] GUniFrac_1.8 broom_1.0.5
## [175] memoise_2.0.1 bslib_0.6.1
## [177] ape_5.7-1