This pipeline is still under development and has not yet been peer reviewed. We appreciate any comments or feedback to help improve this platform, email bealabgithub@gmail.com
LIMON-metab is an expansion of LIMON that allows users to investigate longitudinal changes in microbial metabolic omics networks. This package allows for covariate correction, individual network inferences, and networks per time. In this model, temporal metabolite abundances undergo vsn normalization, are fitted to a user specified linear mixed model to remove covariates,and network inference with Spearman correlation, and finally estimation of individual network properties longitudinally using Linear Interpolation to Obtain Network Estimates for Single Samples (LIONESS). This approach allows users to remove the random effects of repeated samples and sample covariates, return networks per time point, identify interaction changes between each time point, and finally return individual networks and network characteristics per sample/time point. In doing so, LIMON-metab provides a platform to identify the relationship between network centralities and sample features of interest overtime.
You can install LIMON using devtools:
install.packages("devtools")
library(devtools)
install_github("salverna/LIMONmetab")In this script, we will be running the full LIMONmetab pipeline on an infant dataset published in 2019. The metabolomics data were pulled from this summarized article, specifically reference the study by He et al 2019
library(tidyverse)
library(igraph)
library(NBZIMM)
library(SpiecEasi)
library(LIMON)
library(here)
library(lme4)
library(Matrix)
library(MASS)
library(matrixcalc)
library(gridExtra)
library(devtools)
library(reshape2)
library(ggpubr)
library(broom)
library(ggnewscale)
library(coin)
library(parallel)
library(vsn)
library(LIMONmetab)
library(hexbin)
library(ggeffects)Metabolite Data
data("He_2019_metabolites")
metabolite_data <- He_2019_metabolites %>% column_to_rownames("Sample")
# Remove those with the lowest 20% abundance
threshold <- quantile(colSums(metabolite_data), probs = 0.20)
filt_metabolite <- metabolite_data[, colSums(metabolite_data) > threshold]Mean and variance of the data
vsn::meanSdPlot(as.matrix(filt_metabolite))
Metadata
filter down and create some binary
data("He_2019_metadata")
metadata <- He_2019_metadata %>%
dplyr::select(Subject, X, Study.Group, Age, Gender, diet) %>%
group_by(Subject) %>%
# Filter to those with at least three visits
filter(n() >= 4) %>%
ungroup() %>%
column_to_rownames("X") %>%
# add a binary column for Gender
mutate(Gender_binary = case_when(
Gender == "Male" ~ 1,
Gender == "Female" ~ 0)) %>%
# add a binary column for Diet (Formula or Breastmilk)
mutate(Diet_binary = case_when(diet == "Standard infant formula" ~ 0,
diet == "Experimental infant formula" ~ 0,
diet == "Breast milk" ~ 1))# Ensure sample counts and metadata are in the same order
common_samples <- intersect(rownames(filt_metabolite), rownames(metadata))
filt_metabolite <- filt_metabolite[common_samples, ]
metadata <- metadata[common_samples, ]
# Make the Object
Lm_obj <- LIMONm_Obj(Metab = filt_metabolite,
SampleData = metadata)Lm_obj2 <- LIMONm_VSN(Lm_obj, make_pos = FALSE)Re check the mean and variance after normalization
vsn::meanSdPlot(as.matrix(Lm_obj2[["Normalized_Metab"]]))
Distribution Fitting and check
# Set seed
set.seed(12345)
# Fit the distribution/remove covariates
#################################################################################
Lm_obj3 <- LIMONm_DistrFit(Obj = Lm_obj2,
Time = "Age",
Subject = "Subject",
Covariates = c("Gender_binary"),
model = "Gender_binary",
distribution = "LMM")#infer network
Lm_obj4 <- LIMONm_NetInf_Time(Obj = Lm_obj3)# Print Networks
Lm_obj5 <- LIMONm_Edges_Networks(Lm_obj4, threshold = 0.5, vertex.size = 8,
vertex.label.cex = 8, vertex.label.color = "black")
Individualized Networks
# individual Networks
Lm_obj5 <- LIMONm_IndNet(Obj = Lm_obj4)Extract edges - using a stronger threshold
Lm_obj6 <- LIMONm_IndEdges(Lm_obj5, threshold = 0.5)Extract Centralities - using a stronger threshold
Lm_obj7 <- LIMONm_Centralities(Lm_obj6, threshold = 0.5)Next identify which edges are associated with diet group using a multinomial logistic regression. Similar to LIMON, this will only run if there are at least 30 samples for tha time point that have that edge
Lm_obj8 <- LIMONm_StatNodes(Lm_obj7, time = "Age",
dependent = "diet", pval = 0.05, method = "multinom")
# Optional to save output
#saveRDS(Lm_obj8, here("Output", "He_metabolites.rds"))The top edge was 4-Hydroxyphenyllactate-Glycine at time 6. Lets take a look at this edge at 6 months and 12 months by diet group
4-Hydroxyphenyllactate-Glycine time 6
Edge_Table <- Lm_obj8[["Merged_Metabolite_Edge_Table"]]
# list to store all of the model data in
node_results <- list()
# Filter data to a timepoint of interest
edge_data_time_full <- Edge_Table %>% filter(Edge_Table[["Age"]] == 6)
# Filter to an interaction of interest
edge_data_time <- edge_data_time_full %>% filter(Interaction == "4-Hydroxyphenyllactate-Glycine")
dependent <- "diet"
# Run the model
formula <- as.formula(paste(dependent, "~Edge_weight"))
set.seed(12345)
model <- nnet::multinom(formula, data = edge_data_time, trace = FALSE)
model_summary <- broom::tidy(model, conf.int = TRUE)
model_summary$Interaction <- "4-Hydroxyphenyllactate-Glycine"
model_summary$Time_Level <- 6
model_summary$Model_Type <- "Multinomial"
model_summary$Model_SampleSize <- stats::nobs(model)
node_results[[paste("4-Hydroxyphenyllactate-Glycine", 6, sep = "_")]] <- model_summary
# get probabilities
set.seed(12345)
predict(model, newdata = data.frame(Edge_weight = 2), type = 'probs')
#> Breast milk Experimental infant formula
#> 0.983191587 0.001152812
#> Standard infant formula
#> 0.015655600
set.seed(12345)
prob_effects <- ggeffect(model, terms = "Edge_weight[-2:2,by=0.5]")
# Graph probabilities
ggplot(prob_effects) +
aes(x = x, y = predicted, fill = response.level, color = response.level) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 1/2) +
labs(x = '4-Hydroxyphenyllactate-Glycine Edge Weight', y = 'Predicted Probability') +
ylim(c(0, 1)) +
theme_classic() +
scale_color_manual(
values = c("#bc80bd", "#ffed6f", "#ccebc5"),
name = "Infant Diet",
labels = c("Breast Milk" = "Breast.milk",
"Experimental infant formula" = "Experimental.infant.formula",
"Standard infant formula" = "Standard.infant.formula") ) +
scale_fill_manual(
values = c("#bc80bd", "#ffed6f", "#ccebc5"),
name = "Infant Diet",
labels = c("Breast Milk" = "Breast.milk",
"Experimental infant formula" = "Experimental.infant.formula",
"Standard infant formula" = "Standard.infant.formula") ) +
theme(
axis.text.x = element_text(color = "black", family = "Arial", size = 11),
axis.text.y = element_text(color = "black", family = "Arial", size = 11))
# Check model
set.seed(12345)
car::Anova(model)
#> # weights: 6 (2 variable)
#> initial value 46.141716
#> final value 44.559658
#> converged
#> Analysis of Deviance Table (Type II tests)
#>
#> Response: diet
#> LR Chisq Df Pr(>Chisq)
#> Edge_weight 26.247 2 1.998e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14-Hydroxyphenyllactate-Glycine time 12
Edge_Table <- Lm_obj8[["Merged_Metabolite_Edge_Table"]]
# list to store all of the model data in
node_results <- list()
# Filter data to a timepoint of interest
edge_data_time_full <- Edge_Table %>% filter(Edge_Table[["Age"]] == 12)
# Filter to an interaction of interest
edge_data_time <- edge_data_time_full %>% filter(Interaction == "4-Hydroxyphenyllactate-Glycine")
dependent <- "diet"
# Run the model
formula <- as.formula(paste(dependent, "~Edge_weight"))
set.seed(12345)
model <- nnet::multinom(formula, data = edge_data_time, trace = FALSE)
model_summary <- broom::tidy(model, conf.int = TRUE)
model_summary$Interaction <- "4-Hydroxyphenyllactate-Glycine"
model_summary$Time_Level <- 12
model_summary$Model_Type <- "Multinomial"
model_summary$Model_SampleSize <- stats::nobs(model)
node_results[[paste("4-Hydroxyphenyllactate-Glycine", 12, sep = "_")]] <- model_summary
# get probabilities
set.seed(12345)
predict(model, newdata = data.frame(Edge_weight = 2), type = 'probs')
#> Breast milk Experimental infant formula
#> 0.0334799 0.4254905
#> Standard infant formula
#> 0.5410296
set.seed(12345)
prob_effects <- ggeffect(model, terms = "Edge_weight[-2:2,by=0.5]")
# Graph probabilities
ggplot(prob_effects) +
aes(x = x, y = predicted, fill = response.level, color = response.level) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 1/2) +
labs(x = '4-Hydroxyphenyllactate-Glycine Edge Weight', y = 'Predicted Probability') +
ylim(c(0, 1)) +
theme_classic() +
scale_color_manual(
values = c("#bc80bd", "#ffed6f", "#ccebc5"),
name = "Infant Diet",
labels = c("Breast Milk" = "Breast.milk",
"Experimental infant formula" = "Experimental.infant.formula",
"Standard infant formula" = "Standard.infant.formula") ) +
scale_fill_manual(
values = c("#bc80bd", "#ffed6f", "#ccebc5"),
name = "Infant Diet",
labels = c("Breast Milk" = "Breast.milk",
"Experimental infant formula" = "Experimental.infant.formula",
"Standard infant formula" = "Standard.infant.formula") ) +
theme(
axis.text.x = element_text(color = "black", family = "Arial", size = 11),
axis.text.y = element_text(color = "black", family = "Arial", size = 11))
# Check model
set.seed(12345)
car::Anova(model)
#> # weights: 6 (2 variable)
#> initial value 36.254206
#> final value 33.895925
#> converged
#> Analysis of Deviance Table (Type II tests)
#>
#> Response: diet
#> LR Chisq Df Pr(>Chisq)
#> Edge_weight 7.8691 2 0.01955 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LIMÓN-metab identified a stronger positive interaction between 4-Hydroxyphenyllactate and Glycine as being highly predictive of a breast milk diet (BF) at the 6 month period. At 12 months, it was more likely to appear in formula fed infants. The original study researchers found the stool metabolomics between the the two formula groups were very similar while breast fed infants had a more heterogenous stool metabolome. Interestingly, their original analysis of the stool metabolome showed higher levels of Hydroxyphenyllactate in BF infants at 4 and 6 months compared to formula fed (FF) infants, but no difference in its concentration at 12 months. LIMÓN- metab identifying this metabolites interaction with Glycine being predictive of BF at 6 months and FF at 12 months may provide a deeper understanding of the complexity of metabolic alterations in early childhood development.