SECOM Tutorial
knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = NA,
fig.width = 6.25, fig.height = 5)
library(ANCOMBC)
library(tidyverse)1. Introduction
Sparse Estimation of Correlations among Microbiomes (SECOM) (Lin, Eggesbø, and Peddada 2022) is a methodology that aims to detect both linear and nonlinear relationships between a pair of taxa within an ecosystem (e.g., gut) or across ecosystems (e.g., gut and tongue). SECOM corrects both sample-specific and taxon-specific biases and obtains a consistent estimator for the correlation matrix of microbial absolute abundances while maintaining the underlying true sparsity. For more details, please refer to the SECOM paper.
2. Installation
Download package.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ANCOMBC")Load the package.
3. Example Data
The HITChip Atlas dataset contains genus-level microbiota profiling with HITChip for 1006 western adults with no reported health complications, reported in (Lahti et al. 2014). The dataset is available via the microbiome R package (Lahti et al. 2017) in phyloseq (McMurdie and Holmes 2013) format.
data(atlas1006, package = "microbiome")
# Subset to baseline
pseq = phyloseq::subset_samples(atlas1006, time == 0)
# Re-code the bmi group
meta_data = microbiome::meta(pseq)
meta_data$bmi = recode(meta_data$bmi_group,
obese = "obese",
severeobese = "obese",
morbidobese = "obese")
# Note that by default, levels of a categorical variable in R are sorted
# alphabetically. In this case, the reference level for `bmi` will be
# `lean`. To manually change the reference level, for instance, setting `obese`
# as the reference level, use:
meta_data$bmi = factor(meta_data$bmi, levels = c("obese", "overweight", "lean"))
# You can verify the change by checking:
# levels(meta_data$bmi)
# Create the region variable
meta_data$region = recode(as.character(meta_data$nationality),
Scandinavia = "NE", UKIE = "NE", SouthEurope = "SE",
CentralEurope = "CE", EasternEurope = "EE",
.missing = "unknown")
phyloseq::sample_data(pseq) = meta_data
# Subset to lean, overweight, and obese subjects
pseq = phyloseq::subset_samples(pseq, bmi %in% c("lean", "overweight", "obese"))
# Discard "EE" as it contains only 1 subject
# Discard subjects with missing values of region
pseq = phyloseq::subset_samples(pseq, ! region %in% c("EE", "unknown"))
print(pseq)phyloseq-class experiment-level object
otu_table() OTU Table: [ 130 taxa and 873 samples ]
sample_data() Sample Data: [ 873 samples by 12 sample variables ]
tax_table() Taxonomy Table: [ 130 taxa by 3 taxonomic ranks ]4. Run SECOM on a Single Ecosystem
4.1 Run secom functions using the phyloseq object
set.seed(123)
# Linear relationships
res_linear = secom_linear(data = list(pseq), taxa_are_rows = TRUE,
tax_level = "Phylum",
aggregate_data = NULL, meta_data = NULL, pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), method = "pearson",
soft = FALSE, thresh_len = 20, n_cv = 10,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)
# Nonlinear relationships
res_dist = secom_dist(data = list(pseq), taxa_are_rows = TRUE,
tax_level = "Phylum",
aggregate_data = NULL, meta_data = NULL, pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), R = 1000,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)4.2 Visualizations
Pearson correlation with thresholding
corr_linear = res_linear$corr_th
cooccur_linear = res_linear$mat_cooccur
# Filter by co-occurrence
overlap = 10
corr_linear[cooccur_linear < overlap] = 0
df_linear = data.frame(get_upper_tri(corr_linear)) %>%
rownames_to_column("var1") %>%
pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
filter(!is.na(value)) %>%
mutate(value = round(value, 2))
tax_name = sort(union(df_linear$var1, df_linear$var2))
df_linear$var1 = factor(df_linear$var1, levels = tax_name)
df_linear$var2 = factor(df_linear$var2, levels = tax_name)
heat_linear_th = df_linear %>%
ggplot(aes(var2, var1, fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white", na.value = "grey",
midpoint = 0, limit = c(-1,1), space = "Lab",
name = NULL) +
scale_x_discrete(drop = FALSE) +
scale_y_discrete(drop = FALSE) +
geom_text(aes(var2, var1, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL, title = "Pearson (Thresholding)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1,
face = "italic"),
axis.text.y = element_text(size = 12, face = "italic"),
strip.text.x = element_text(size = 14),
strip.text.y = element_text(size = 14),
legend.text = element_text(size = 12),
plot.title = element_text(hjust = 0.5, size = 15),
panel.grid.major = element_blank(),
axis.ticks = element_blank(),
legend.position = "none") +
coord_fixed()
heat_linear_th
Pearson correlation with p-value filtering
corr_linear = res_linear$corr_fl
cooccur_linear = res_linear$mat_cooccur
# Filter by co-occurrence
overlap = 10
corr_linear[cooccur_linear < overlap] = 0
df_linear = data.frame(get_upper_tri(corr_linear)) %>%
rownames_to_column("var1") %>%
pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
filter(!is.na(value)) %>%
mutate(value = round(value, 2))
tax_name = sort(union(df_linear$var1, df_linear$var2))
df_linear$var1 = factor(df_linear$var1, levels = tax_name)
df_linear$var2 = factor(df_linear$var2, levels = tax_name)
heat_linear_fl = df_linear %>%
ggplot(aes(var2, var1, fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white", na.value = "grey",
midpoint = 0, limit = c(-1,1), space = "Lab",
name = NULL) +
scale_x_discrete(drop = FALSE) +
scale_y_discrete(drop = FALSE) +
geom_text(aes(var2, var1, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL, title = "Pearson (Filtering)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1,
face = "italic"),
axis.text.y = element_text(size = 12, face = "italic"),
strip.text.x = element_text(size = 14),
strip.text.y = element_text(size = 14),
legend.text = element_text(size = 12),
plot.title = element_text(hjust = 0.5, size = 15),
panel.grid.major = element_blank(),
axis.ticks = element_blank(),
legend.position = "none") +
coord_fixed()
heat_linear_fl
Distance correlation with p-value filtering
corr_dist = res_dist$dcorr_fl
cooccur_dist = res_dist$mat_cooccur
# Filter by co-occurrence
overlap = 10
corr_dist[cooccur_dist < overlap] = 0
df_dist = data.frame(get_upper_tri(corr_dist)) %>%
rownames_to_column("var1") %>%
pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
filter(!is.na(value)) %>%
mutate(value = round(value, 2))
tax_name = sort(union(df_dist$var1, df_dist$var2))
df_dist$var1 = factor(df_dist$var1, levels = tax_name)
df_dist$var2 = factor(df_dist$var2, levels = tax_name)
heat_dist_fl = df_dist %>%
ggplot(aes(var2, var1, fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white", na.value = "grey",
midpoint = 0, limit = c(-1,1), space = "Lab",
name = NULL) +
scale_x_discrete(drop = FALSE) +
scale_y_discrete(drop = FALSE) +
geom_text(aes(var2, var1, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL, title = "Distance (Filtering)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1,
face = "italic"),
axis.text.y = element_text(size = 12, face = "italic"),
strip.text.x = element_text(size = 14),
strip.text.y = element_text(size = 14),
legend.text = element_text(size = 12),
plot.title = element_text(hjust = 0.5, size = 15),
panel.grid.major = element_blank(),
axis.ticks = element_blank(),
legend.position = "none") +
coord_fixed()
heat_dist_fl
4.3 Run secom functions using the tse object
One can also run SECOM function using the
TreeSummarizedExperiment object.
tse = mia::makeTreeSummarizedExperimentFromPhyloseq(atlas1006)
tse = tse[, tse$time == 0]
tse$bmi = recode(tse$bmi_group,
obese = "obese",
severeobese = "obese",
morbidobese = "obese")
tse = tse[, tse$bmi %in% c("lean", "overweight", "obese")]
tse$bmi = factor(tse$bmi, levels = c("obese", "overweight", "lean"))
tse$region = recode(as.character(tse$nationality),
Scandinavia = "NE", UKIE = "NE", SouthEurope = "SE",
CentralEurope = "CE", EasternEurope = "EE",
.missing = "unknown")
tse = tse[, ! tse$region %in% c("EE", "unknown")]
set.seed(123)
# Linear relationships
res_linear = secom_linear(data = list(tse), taxa_are_rows = TRUE,
assay_name = "counts", tax_level = "Phylum",
aggregate_data = NULL, meta_data = NULL, pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), method = "pearson",
soft = FALSE, thresh_len = 20, n_cv = 10,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)
# Nonlinear relationships
res_dist = secom_dist(data = list(tse), taxa_are_rows = TRUE,
assay_name = "counts", tax_level = "Phylum",
aggregate_data = NULL, meta_data = NULL, pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), R = 1000,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)4.4 Run secom functions by directly providing the abundance and metadata
Please ensure that you provide the following: 1) The abundance data at its lowest possible taxonomic level. 2) The aggregated data at the desired taxonomic level; if no aggregation is performed, it can be the same as the original abundance data. 3) The sample metadata.
abundance_data = microbiome::abundances(pseq)
aggregate_data = microbiome::abundances(microbiome::aggregate_taxa(pseq, "Phylum"))
meta_data = microbiome::meta(pseq)
set.seed(123)
# Linear relationships
res_linear = secom_linear(data = list(abundance_data),
taxa_are_rows = TRUE,
aggregate_data = list(aggregate_data),
meta_data = list(meta_data),
pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), method = "pearson",
soft = FALSE, thresh_len = 20, n_cv = 10,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)
# Nonlinear relationships
res_dist = secom_dist(data = list(abundance_data),
taxa_are_rows = TRUE,
aggregate_data = list(aggregate_data),
meta_data = list(meta_data),
pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), R = 1000,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)5. Run SECOM on Multiple Ecosystems
5.1 Data manipulation
To compute correlations whithin and across different ecosystems, one needs to make sure that there are samples in common across these ecosystems.
# Select subjects from "CE" and "NE"
pseq1 = phyloseq::subset_samples(pseq, region == "CE")
pseq2 = phyloseq::subset_samples(pseq, region == "NE")
phyloseq::sample_names(pseq1) = paste0("Sample-", seq_len(phyloseq::nsamples(pseq1)))
phyloseq::sample_names(pseq2) = paste0("Sample-", seq_len(phyloseq::nsamples(pseq2)))
print(pseq1)phyloseq-class experiment-level object
otu_table() OTU Table: [ 130 taxa and 578 samples ]
sample_data() Sample Data: [ 578 samples by 12 sample variables ]
tax_table() Taxonomy Table: [ 130 taxa by 3 taxonomic ranks ]phyloseq-class experiment-level object
otu_table() OTU Table: [ 130 taxa and 181 samples ]
sample_data() Sample Data: [ 181 samples by 12 sample variables ]
tax_table() Taxonomy Table: [ 130 taxa by 3 taxonomic ranks ]5.2 Run secom functions using the phyloseq object
set.seed(123)
# Linear relationships
res_linear = secom_linear(data = list(CE = pseq1, NE = pseq2),
taxa_are_rows = TRUE,
tax_level = c("Phylum", "Phylum"),
aggregate_data = NULL, meta_data = NULL, pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), method = "pearson",
soft = FALSE, thresh_len = 20, n_cv = 10,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)
# Nonlinear relationships
res_dist = secom_dist(data = list(CE = pseq1, NE = pseq2),
taxa_are_rows = TRUE,
tax_level = c("Phylum", "Phylum"),
aggregate_data = NULL, meta_data = NULL, pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), R = 1000,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)5.3 Visualizations
Pearson correlation with thresholding
corr_linear = res_linear$corr_th
cooccur_linear = res_linear$mat_cooccur
# Filter by co-occurrence
overlap = 10
corr_linear[cooccur_linear < overlap] = 0
df_linear = data.frame(get_upper_tri(corr_linear)) %>%
rownames_to_column("var1") %>%
pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
filter(!is.na(value)) %>%
mutate(var2 = gsub("\\...", " - ", var2),
value = round(value, 2))
tax_name = sort(union(df_linear$var1, df_linear$var2))
df_linear$var1 = factor(df_linear$var1, levels = tax_name)
df_linear$var2 = factor(df_linear$var2, levels = tax_name)
txt_color = ifelse(grepl("CE", tax_name), "#1B9E77", "#D95F02")
heat_linear_th = df_linear %>%
ggplot(aes(var2, var1, fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
na.value = "grey", midpoint = 0, limit = c(-1,1),
space = "Lab", name = NULL) +
scale_x_discrete(drop = FALSE) +
scale_y_discrete(drop = FALSE) +
geom_text(aes(var2, var1, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL, title = "Pearson (Thresholding)") +
theme_bw() +
geom_vline(xintercept = 6.5, color = "blue", linetype = "dashed") +
geom_hline(yintercept = 6.5, color = "blue", linetype = "dashed") +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1,
face = "italic", color = txt_color),
axis.text.y = element_text(size = 12, face = "italic",
color = txt_color),
strip.text.x = element_text(size = 14),
strip.text.y = element_text(size = 14),
legend.text = element_text(size = 12),
plot.title = element_text(hjust = 0.5, size = 15),
panel.grid.major = element_blank(),
axis.ticks = element_blank(),
legend.position = "none") +
coord_fixed()
heat_linear_th
Pearson correlation with p-value filtering
corr_linear = res_linear$corr_th
cooccur_linear = res_linear$mat_cooccur
# Filter by co-occurrence
overlap = 10
corr_linear[cooccur_linear < overlap] = 0
df_linear = data.frame(get_upper_tri(corr_linear)) %>%
rownames_to_column("var1") %>%
pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
filter(!is.na(value)) %>%
mutate(var2 = gsub("\\...", " - ", var2),
value = round(value, 2))
tax_name = sort(union(df_linear$var1, df_linear$var2))
df_linear$var1 = factor(df_linear$var1, levels = tax_name)
df_linear$var2 = factor(df_linear$var2, levels = tax_name)
txt_color = ifelse(grepl("CE", tax_name), "#1B9E77", "#D95F02")
heat_linear_fl = df_linear %>%
ggplot(aes(var2, var1, fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
na.value = "grey", midpoint = 0, limit = c(-1,1),
space = "Lab", name = NULL) +
scale_x_discrete(drop = FALSE) +
scale_y_discrete(drop = FALSE) +
geom_text(aes(var2, var1, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL, title = "Pearson (Filtering)") +
theme_bw() +
geom_vline(xintercept = 6.5, color = "blue", linetype = "dashed") +
geom_hline(yintercept = 6.5, color = "blue", linetype = "dashed") +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1,
face = "italic", color = txt_color),
axis.text.y = element_text(size = 12, face = "italic",
color = txt_color),
strip.text.x = element_text(size = 14),
strip.text.y = element_text(size = 14),
legend.text = element_text(size = 12),
plot.title = element_text(hjust = 0.5, size = 15),
panel.grid.major = element_blank(),
axis.ticks = element_blank(),
legend.position = "none") +
coord_fixed()
heat_linear_fl
Distance correlation with p-value filtering
corr_dist = res_dist$dcorr_fl
cooccur_dist = res_dist$mat_cooccur
# Filter by co-occurrence
overlap = 10
corr_dist[cooccur_dist < overlap] = 0
df_dist = data.frame(get_upper_tri(corr_dist)) %>%
rownames_to_column("var1") %>%
pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
filter(!is.na(value)) %>%
mutate(var2 = gsub("\\...", " - ", var2),
value = round(value, 2))
tax_name = sort(union(df_dist$var1, df_dist$var2))
df_dist$var1 = factor(df_dist$var1, levels = tax_name)
df_dist$var2 = factor(df_dist$var2, levels = tax_name)
txt_color = ifelse(grepl("CE", tax_name), "#1B9E77", "#D95F02")
heat_dist_fl = df_dist %>%
ggplot(aes(var2, var1, fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
na.value = "grey", midpoint = 0, limit = c(-1,1),
space = "Lab", name = NULL) +
scale_x_discrete(drop = FALSE) +
scale_y_discrete(drop = FALSE) +
geom_text(aes(var2, var1, label = value), color = "black", size = 4) +
labs(x = NULL, y = NULL, title = "Distance (Filtering)") +
theme_bw() +
geom_vline(xintercept = 6.5, color = "blue", linetype = "dashed") +
geom_hline(yintercept = 6.5, color = "blue", linetype = "dashed") +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1,
face = "italic", color = txt_color),
axis.text.y = element_text(size = 12, face = "italic",
color = txt_color),
strip.text.x = element_text(size = 14),
strip.text.y = element_text(size = 14),
legend.text = element_text(size = 12),
plot.title = element_text(hjust = 0.5, size = 15),
panel.grid.major = element_blank(),
axis.ticks = element_blank(),
legend.position = "none") +
coord_fixed()
heat_dist_fl
5.4 Run secom functions using the tse object
One can also run SECOM function using the
TreeSummarizedExperiment object.
# Select subjects from "CE" and "NE"
tse1 = tse[, tse$region == "CE"]
tse2 = tse[, tse$region == "NE"]
# Rename samples to ensure there is an overlap of samples between CE and NE
colnames(tse1) = paste0("Sample-", seq_len(ncol(tse1)))
colnames(tse2) = paste0("Sample-", seq_len(ncol(tse2)))
set.seed(123)
# Linear relationships
res_linear = secom_linear(data = list(CE = tse1, NE = tse2),
taxa_are_rows = TRUE,
assay_name = c("counts", "counts"),
tax_level = c("Phylum", "Phylum"),
aggregate_data = NULL, meta_data = NULL, pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), method = "pearson",
soft = FALSE, thresh_len = 20, n_cv = 10,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)
# Nonlinear relationships
res_dist = secom_dist(data = list(CE = tse1, NE = tse2),
taxa_are_rows = TRUE,
assay_name = c("counts", "counts"),
tax_level = c("Phylum", "Phylum"),
aggregate_data = NULL, meta_data = NULL, pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), R = 1000,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)5.5 Run secom functions by directly providing the abundance and metadata
Please ensure that you provide the following: 1) The abundance data at its lowest possible taxonomic level. 2) The aggregated data at the desired taxonomic level; if no aggregation is performed, it can be the same as the original abundance data. 3) The sample metadata.
ce_idx = which(meta_data$region == "CE")
ne_idx = which(meta_data$region == "NE")
abundance_data1 = abundance_data[, ce_idx]
abundance_data2 = abundance_data[, ne_idx]
aggregate_data1 = aggregate_data[, ce_idx]
aggregate_data2 = aggregate_data[, ne_idx]
meta_data1 = meta_data[ce_idx, ]
meta_data2 = meta_data[ne_idx, ]
sample_size1 = ncol(abundance_data1)
sample_size2 = ncol(abundance_data2)
colnames(abundance_data1) = paste0("Sample-", seq_len(sample_size1))
colnames(abundance_data2) = paste0("Sample-", seq_len(sample_size2))
rownames(meta_data1) = paste0("Sample-", seq_len(sample_size1))
rownames(meta_data2) = paste0("Sample-", seq_len(sample_size2))
set.seed(123)
# Linear relationships
res_linear = secom_linear(data = list(CE = abundance_data1, NE = abundance_data2),
taxa_are_rows = TRUE,
aggregate_data = list(aggregate_data1, aggregate_data2),
meta_data = list(meta_data1, meta_data2),
pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), method = "pearson",
soft = FALSE, thresh_len = 20, n_cv = 10,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)
# Nonlinear relationships
res_dist = secom_dist(data = list(CE = abundance_data1, NE = abundance_data2),
taxa_are_rows = TRUE,
aggregate_data = list(aggregate_data1, aggregate_data2),
meta_data = list(meta_data1, meta_data2),
pseudo = 0,
prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
wins_quant = c(0.05, 0.95), R = 1000,
thresh_hard = 0.3, max_p = 0.005, n_cl = 2)Session information
R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
[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: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] doRNG_1.8.6 rngtools_1.5.2 foreach_1.5.2 DT_0.33
[5] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[9] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[13] ggplot2_3.5.1 tidyverse_2.0.0 ANCOMBC_2.9.0 rmarkdown_2.29
loaded via a namespace (and not attached):
[1] sys_3.4.3 rstudioapi_0.17.1 jsonlite_1.8.9
[4] magrittr_2.0.3 TH.data_1.1-2 farver_2.1.2
[7] nloptr_2.1.1 zlibbioc_1.52.0 vctrs_0.6.5
[10] multtest_2.63.0 minqa_1.2.8 base64enc_0.1-3
[13] htmltools_0.5.8.1 energy_1.7-12 haven_2.5.4
[16] Rhdf5lib_1.29.0 cellranger_1.1.0 Formula_1.2-5
[19] rhdf5_2.51.0 sass_0.4.9 bslib_0.8.0
[22] htmlwidgets_1.6.4 plyr_1.8.9 sandwich_3.1-1
[25] rootSolve_1.8.2.4 zoo_1.8-12 cachem_1.1.0
[28] buildtools_1.0.0 igraph_2.1.1 lifecycle_1.0.4
[31] iterators_1.0.14 pkgconfig_2.0.3 Matrix_1.7-1
[34] R6_2.5.1 fastmap_1.2.0 GenomeInfoDbData_1.2.13
[37] rbibutils_2.3 digest_0.6.37 Exact_3.3
[40] numDeriv_2016.8-1.1 colorspace_2.1-1 S4Vectors_0.45.2
[43] crosstalk_1.2.1 Hmisc_5.2-0 vegan_2.6-8
[46] labeling_0.4.3 fansi_1.0.6 timechange_0.3.0
[49] mgcv_1.9-1 httr_1.4.7 compiler_4.4.2
[52] proxy_0.4-27 bit64_4.5.2 withr_3.0.2
[55] doParallel_1.0.17 gsl_2.1-8 htmlTable_2.4.3
[58] backports_1.5.0 MASS_7.3-61 biomformat_1.35.0
[61] permute_0.9-7 gtools_3.9.5 CVXR_1.0-15
[64] gld_2.6.6 tools_4.4.2 foreign_0.8-87
[67] ape_5.8 nnet_7.3-19 glue_1.8.0
[70] rhdf5filters_1.19.0 nlme_3.1-166 grid_4.4.2
[73] Rtsne_0.17 checkmate_2.3.2 cluster_2.1.6
[76] reshape2_1.4.4 ade4_1.7-22 generics_0.1.3
[79] microbiome_1.29.0 gtable_0.3.6 tzdb_0.4.0
[82] class_7.3-22 data.table_1.16.2 lmom_3.2
[85] hms_1.1.3 utf8_1.2.4 XVector_0.47.0
[88] BiocGenerics_0.53.3 pillar_1.9.0 splines_4.4.2
[91] lattice_0.22-6 survival_3.7-0 gmp_0.7-5
[94] bit_4.5.0 tidyselect_1.2.1 maketools_1.3.1
[97] Biostrings_2.75.1 knitr_1.49 gridExtra_2.3
[100] phyloseq_1.51.0 IRanges_2.41.1 stats4_4.4.2
[103] xfun_0.49 expm_1.0-0 Biobase_2.67.0
[106] stringi_1.8.4 UCSC.utils_1.3.0 yaml_2.3.10
[109] boot_1.3-31 evaluate_1.0.1 codetools_0.2-20
[112] cli_3.6.3 rpart_4.1.23 DescTools_0.99.58
[115] Rdpack_2.6.2 munsell_0.5.1 jquerylib_0.1.4
[118] Rcpp_1.0.13-1 GenomeInfoDb_1.43.1 readxl_1.4.3
[121] parallel_4.4.2 lme4_1.1-35.5 Rmpfr_1.0-0
[124] mvtnorm_1.3-2 lmerTest_3.1-3 scales_1.3.0
[127] e1071_1.7-16 crayon_1.5.3 rlang_1.1.4
[130] multcomp_1.4-26 References
Lahti, Leo, Jarkko Salojärvi, Anne Salonen, Marten Scheffer, and Willem
M De Vos. 2014. “Tipping Elements in the Human Intestinal
Ecosystem.” Nature Communications 5 (1): 1–10.
Lahti, Leo, Sudarshan Shetty, T Blake, J Salojarvi, et al. 2017.
“Tools for Microbiome Analysis in r.” Version 1:
10013.
Lin, Huang, Merete Eggesbø, and Shyamal Das Peddada. 2022. “Linear
and Nonlinear Correlation Estimators Unveil Undescribed Taxa
Interactions in Microbiome Data.” Nature Communications
13 (1): 1–16.
McMurdie, Paul J, and Susan Holmes. 2013. “Phyloseq: An r Package
for Reproducible Interactive Analysis and Graphics of Microbiome Census
Data.” PloS One 8 (4): e61217.