Targeted stool metabolomics suggests exploratory catecholamine- and tryptophan-linked metabolic features in autism spectrum disorder.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches

[1] § Materials and methods › Statistical analysis › Machine learning ↔ statistical_analysis.Rmd, lines 953 to 1036 · score 0.77 · Brier score, predicted probabilities, positive class, Calibration, LOESS, smoothed · Report a mismatch
[2] § Results › Discriminative performance of metabolite features in classifying ASD ↔ statistical_analysis.Rmd, lines 491 to 587 · score 0.62 · cross validated, positive class, random forest, predicting ASD, ROC, AUC · Report a mismatch
[3] § Results › Model calibration performance ↔ statistical_analysis.Rmd, lines 953 to 1036 · score 0.62 · Brier score, predicted probabilities, calibration, LOESS, smoothed, curves · Report a mismatch
[4] § Materials and methods › Statistical analysis › Machine learning ↔ statistical_analysis.Rmd, lines 491 to 587 · score 0.60 · Random forest classifiers, cross validation, seed, trees, workflow, fold · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 1,040 lines · 30 KB · MIT · 4 matches

---
title: "Targeted ASD Stool Metabolomics Analysis"
author: "Kevin Liu"
date: "`r Sys.Date()`"
output:
pdf_document:
toc: true
keep_tex: true
---
```{r setup, include=FALSE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(scipen = 999)
library(tidyverse)
library(egg)
library(rstatix)
library(tidymodels)
library(yardstick)
library(ranger)
library(pROC)
library(patchwork)
library(scales)
library(broom)
library(ggpubr)
theme_set(theme_article() +
theme(aspect.ratio = 1))
```
```{r read_data, message=FALSE}
clin_data = read_csv("data/clin_data.csv")
avg_data = read_csv("data/data_averaged.csv") %>%
left_join(clin_data, by = "subject_id") %>%
relocate(subject_id, group, sex, age) %>%
mutate(group = factor(group, levels = c("NAC", "ASD")),
sex = factor(sex, levels = c("Female", "Male")))
norm_data = read_csv("data/data_normalized.csv") %>%
left_join(clin_data, by = "subject_id") %>%
relocate(subject_id, group, sex, age) %>%
mutate(group = factor(group, levels = c("NAC", "ASD")),
sex = factor(sex, levels = c("Female", "Male")))
name_map = read_csv("data/name_map.csv", na = c("", "NaN"))
levels(avg_data$group)
levels(norm_data$group)
```
```{r descriptive_stats}
# Figure S1
avg_data %>%
ggplot(aes(x = age, fill = group)) +
geom_histogram(data = subset(avg_data, group == "NAC"),
aes(y = after_stat(count)), binwidth = 1, boundary = 0, closed = "left") +
geom_histogram(data = subset(avg_data, group == "ASD"),
aes(y = -after_stat(count)), binwidth = 1, boundary = 0, closed = "left") +
scale_y_continuous(labels = abs, limits = c(-10, 10)) +
scale_fill_manual(values = c("NAC" = "#21918c", "ASD" = "#fde725")) +
labs(x = "Age (years)", y = "Number of Subjects", fill = "Group") +
theme(legend.position = "bottom") +
geom_hline(yintercept = 0, linewidth = 0.25) +
coord_flip()
age_tbl = avg_data %>%
group_by(group) %>%
summarise(
value = sprintf(
"%.2f ± %.2f, [%.2f, %.2f]",
mean(age, na.rm = TRUE),
sd(age, na.rm = TRUE),
min(age, na.rm = TRUE),
max(age, na.rm = TRUE)
),
.groups = "drop"
) %>%
pivot_wider(names_from = group, values_from = value) %>%
mutate(
Characteristic = "Age (mean ± SD, [min, max])")
sex_tbl = avg_data %>%
count(group, sex) %>%
group_by(group) %>%
mutate(value = sprintf("%d (%.1f%%)", n, 100 * n / sum(n))) %>%
ungroup() %>%
select(group, sex, value) %>%
pivot_wider(names_from = group, values_from = value) %>%
mutate(Characteristic = as.character(sex)) %>%
select(-sex)
wilcox.test(age ~ group, data = avg_data)
fisher.test(table(avg_data$group, avg_data$sex))
# Table 1
bind_rows(age_tbl, sex_tbl) %>%
select(Characteristic, NAC, ASD) %>%
write.csv("output_data/1_demographics.csv")
```
```{r get_mol_summary}
num_cols = avg_data %>%
select(where(is.numeric), -subject_id, -age) %>%
names()
# Table 2
avg_data %>%
select(group, all_of(num_cols)) %>%
pivot_longer(
cols = -group,
names_to = "variable",
values_to = "value"
) %>%
group_by(group, variable) %>%
summarise(
n = n(),
n_zero = sum(value == 0, na.rm = TRUE),
pct_zero = n_zero / n * 100,
n_nonzero = sum(value != 0 & !is.na(value)),
# stats excluding zeros
mean = mean(if_else(value == 0, NA_real_, value), na.rm = TRUE),
sd = sd(if_else(value == 0, NA_real_, value), na.rm = TRUE),
min = min(if_else(value == 0, NA_real_, value), na.rm = TRUE),
max = max(if_else(value == 0, NA_real_, value), na.rm = TRUE),
median = median(if_else(value == 0, NA_real_, value), na.rm = TRUE),
.groups = "drop"
) %>%
ungroup() %>%
mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
mutate(range = str_c("[", min, ", ", max, "]"),
mean_sd = str_c(mean, " (", sd, ")"),
zeros = str_c(n_zero, "/", n, " (", pct_zero, "%)")) %>%
select(group, variable, mean_sd, range, median, zeros) %>%
write_csv("output_data/2_mol_summary.csv")
```
```{r viz_data, fig.dpi=300, fig.width=8}
# Figure 1
avg_data %>%
pivot_longer(cols = starts_with("mol_"), names_to = "molecule") %>%
left_join(name_map %>% select(mol_name, compound_acronym), by = c("molecule" = "mol_name")) %>%
mutate(molecule = compound_acronym) %>%
select(-compound_acronym) %>%
ggplot(aes(x = group, y = value, fill = group)) +
geom_violin() +
geom_jitter(width = 0.2, alpha = 0.5) +
facet_wrap(~ molecule, scales = "free_y", ncol = 6) +
labs(x = "", y = "Concentration (nM)",
fill = "Group") +
theme(legend.position = "bottom", aspect.ratio = 1,
strip.text = element_text(size = 6)) +
scale_fill_manual(values = c("#21918c", "#fde725"))
```
```{r corr_matrix, fig.width=10}
num_vars = c("age", grep("^mol_", names(avg_data), value = TRUE))
name_lookup = setNames(c(name_map$compound_acronym, "Age"),
c(name_map$mol_name, "age"))
axis_levels = num_vars
corr_lower_by_group = function(df_group, group_label) {
mat = df_group %>%
select(all_of(num_vars)) %>%
as.matrix()
rc = Hmisc::rcorr(mat, type = "spearman")
R = rc$r
P = rc$P
vars = colnames(R)
idx = which(lower.tri(R), arr.ind = TRUE)
tibble(
Group = group_label,
Var1 = vars[idx[,1]],
Var2 = vars[idx[,2]],
r = as.numeric(R[idx]),
p = as.numeric(P[idx])
)
}
cor_df = bind_rows(
corr_lower_by_group(filter(avg_data, group == "ASD"), "ASD"),
corr_lower_by_group(filter(avg_data, group == "NAC"), "NAC")
) %>%
group_by(Group) %>%
mutate(
p_adj = p.adjust(p, method = "fdr"), # FDR within each group
p_cap = pmax(p_adj, 1e-300), # avoid -Inf for -log10
Var1 = factor(Var1, levels = axis_levels),
Var2 = factor(Var2, levels = axis_levels),
Significant = !is.na(p_adj) & p_adj < 0.05
) %>%
ungroup() %>%
mutate(
Var1_label = recode(Var1, !!!name_lookup),
Var2_label = recode(Var2, !!!name_lookup)
)
# Figure 5
cor_df %>%
ggplot(aes(x = Var1_label, y = Var2_label)) +
geom_point(aes(fill = r, size = -log10(p_cap)),
shape = 21, color = "grey60", stroke = 0.1, na.rm = TRUE) +
geom_point(data = filter(cor_df, Significant),
aes(size = -log10(p_cap)),
shape = 21, fill = NA, color = "black", stroke = 0.5, na.rm = TRUE) +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0, limits = c(-1, 1),
oob = scales::squish, name = "Spearman r") +
scale_size(range = c(1.5, 6), name = NULL, guide = "none") +
coord_fixed() +
facet_wrap(~ Group, ncol = 2) +
theme_light() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "right"
) +
labs(x = NULL, y = NULL)
# Table S5
cor_df %>%
filter(p_adj < 0.05) %>%
select(Group, Var1_label, Var2_label, r, p, p_adj) %>%
mutate(
r = signif(r, 3),
p = formatC(p, format = "e", digits = 2),
p_adj = formatC(p_adj, format = "e", digits = 2)
) %>%
write_csv("output_data/s5_groupwise_sigif_corr.csv")
```
```{r corr_scatterplots_all_significant_patchwork, fig.dpi=300, fig.width=10, fig.height=12}
# Supplementary Figure 2
# Scatterplots for all FDR-significant within-group correlations
cor_df_labeled = cor_df %>%
select(-any_of(c("Var1_label", "Var2_label"))) %>%
left_join(
name_map %>% select(mol_name, compound_acronym),
by = c("Var1" = "mol_name")
) %>%
rename(Var1_label = compound_acronym) %>%
left_join(
name_map %>% select(mol_name, compound_acronym),
by = c("Var2" = "mol_name")
) %>%
rename(Var2_label = compound_acronym) %>%
mutate(
Var1 = as.character(Var1),
Var2 = as.character(Var2),
Var1_label = if_else(Var1 == "age", "Age", Var1_label),
Var2_label = if_else(Var2 == "age", "Age", Var2_label)
)
panel_order = tribble(
~Group, ~Var1, ~Var2, ~order,
"ASD", "mol_trp", "mol_niacin", 1,
"NAC", "mol_trp", "mol_kyn", 2,
"ASD", "mol_trp", "mol_gaba", 3,
"ASD", "mol_trp", "mol_ne", 4,
"ASD", "mol_trp", "mol_aea", 5,
"NAC", "mol_trp", "mol_aea", 6,
"ASD", "mol_ne", "mol_aea", 7,
"ASD", "mol_niacin", "mol_gaba", 8,
"ASD", "mol_3nt", "mol_ua", 9,
"ASD", "mol_lca", "mol_dca", 10,
"NAC", "mol_lca", "mol_dca", 11
)
sig_corr_stats = cor_df_labeled %>%
filter(Significant == TRUE) %>%
inner_join(panel_order, by = c("Group", "Var1", "Var2")) %>%
arrange(order) %>%
mutate(
q_label = case_when(
p_adj < 0.001 ~ "q < 0.001",
TRUE ~ str_c("q = ", sprintf("%.3f", p_adj))
),
stat_label = str_c(
"Spearman r = ", sprintf("%.2f", r),
"\nFDR ", q_label
),
panel_title = str_c(Group, ": ", Var2_label, " vs ", Var1_label)
)
plot_sig_corr_pair = function(group_label, var1, var2, var1_label, var2_label,
stat_label, panel_title) {
plot_df = avg_data %>%
filter(group == group_label) %>%
transmute(
Group = group,
x = .data[[var2]],
y = .data[[var1]]
) %>%
filter(x > 0, y > 0)
label_df = plot_df %>%
summarise(
x = 10^(min(log10(x), na.rm = TRUE) + 0.05 * diff(range(log10(x), na.rm = TRUE))),
y = 10^(max(log10(y), na.rm = TRUE) - 0.05 * diff(range(log10(y), na.rm = TRUE)))
) %>%
mutate(label = stat_label)
ggplot(plot_df, aes(x = x, y = y)) +
geom_point(
aes(fill = Group),
shape = 21,
color = "black",
stroke = 0.2,
size = 1.8,
alpha = 0.8
) +
geom_smooth(
method = "loess",
se = FALSE,
linewidth = 0.5,
color = "black") +
geom_text(
data = label_df,
aes(x = x, y = y, label = label),
inherit.aes = FALSE,
hjust = 0,
vjust = 1,
size = 2.5,
lineheight = 0.95
) +
scale_x_log10() +
scale_y_log10() +
scale_fill_manual(values = c("NAC" = "#21918c", "ASD" = "#fde725")) +
labs(
title = panel_title,
x = str_c(var2_label, " (nM, log10)"),
y = str_c(var1_label, " (nM, log10)"),
fill = "Group"
) +
theme(
plot.title = element_text(size = 8),
axis.title = element_text(size = 7),
axis.text = element_text(size = 6),
aspect.ratio = 1
)
}
sig_corr_plots = pmap(
sig_corr_stats %>%
select(
group_label = Group,
var1 = Var1,
var2 = Var2,
var1_label = Var1_label,
var2_label = Var2_label,
stat_label,
panel_title
),
plot_sig_corr_pair
)
wrap_plots(sig_corr_plots, ncol = 3, guides = "collect") &
theme(
legend.position = "bottom",
legend.title = element_text(size = 9),
legend.text = element_text(size = 9)
)
```
```{r univar_tests}
# norm_data %>%
# pivot_longer(cols = starts_with("mol_"), names_to = "metabolite") %>%
# group_by(metabolite) %>%
# wilcox_test(value ~ group) %>%
# adjust_pvalue(method = "fdr") %>%
# arrange(p) %>%
# left_join(name_map %>% select(mol_name, compound_name), by = c("metabolite" = "mol_name"))
```
```{r adj_univar_tests}
adj_univar_lm = norm_data %>%
pivot_longer(cols = starts_with("mol_"), names_to = "molecule", values_to = "value") %>%
group_by(molecule) %>%
do(tidy(lm(value ~ group + age + sex, data = .))) %>%
ungroup() %>%
filter(term == "groupASD") %>%
mutate(p_adj = p.adjust(p.value, method = "fdr")) %>%
arrange(p.value) %>%
select(molecule, estimate, std.error, statistic, p.value, p_adj) %>%
left_join(name_map %>% select(mol_name, compound_acronym), by = c("molecule" = "mol_name")) %>%
mutate(molecule = compound_acronym) %>%
select(-compound_acronym)
```
```{r log2fc, fig.dpi=300}
log2fc_result = norm_data %>%
pivot_longer(cols = starts_with("mol_"), names_to = "molecule", values_to = "value") %>%
group_by(group, molecule) %>%
summarise(mean_value = mean(value, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(names_from = group, values_from = mean_value) %>%
mutate(
log2_fc = (ASD - NAC) * log2(10) # Convert log10 difference to log2 scale
) %>%
arrange(desc(log2_fc)) %>%
left_join(name_map %>% select(mol_name, compound_acronym), by = c("molecule" = "mol_name")) %>%
mutate(molecule = compound_acronym) %>%
select(-compound_acronym)
# Figure 2
log2fc_result %>%
mutate(molecule = fct_rev(fct_inorder(molecule))) %>%
ggplot(aes(x = molecule, y = log2_fc, fill = log2_fc > 0)) +
geom_col() +
coord_flip() +
scale_fill_manual(values = c("TRUE" = "#fde725", "FALSE" = "#21918c"),
labels = c("NAC > ASD", "ASD > NAC")) +
labs(
title = "Log2 Fold Change (ASD vs. NAC)",
x = NULL,
y = "Log2 Fold Change",
fill = "Direction"
) +
ylim(-2.25, 2.25) +
theme(legend.position = "bottom") +
geom_hline(yintercept = 1, linetype = "dashed", color = "gray") +
geom_hline(yintercept = -1, linetype = "dashed", color = "gray")
# Table S2
log2fc_result %>%
left_join(adj_univar_lm, by = "molecule") %>%
mutate(across(where(is.numeric), \(x) round(x, 4))) %>%
arrange(p.value) %>%
left_join(name_map %>% select(compound_acronym, compound_name), by = c("molecule" = "compound_acronym")) %>%
mutate(molecule = str_c(compound_name, " (", molecule, ")")) %>%
select(-compound_name) %>%
write_csv("output_data/s2_log2fc_adj_univar.csv")
```
```{r ne_sensitivity, message=FALSE, fig.width=5, fig.height=4}
# NE sensitivity check
# # NE distribution by group
# ggplot(avg_data, aes(x = group, y = mol_ne, color = group)) +
# geom_boxplot(outlier.shape = NA) +
# geom_jitter(width = 0.1, alpha = 0.7, size = 2) +
# labs(x = NULL, y = "Norepinephrine (nM)")
# Top NE values
ne_top = avg_data %>%
arrange(desc(mol_ne)) %>%
select(subject_id, group, sex, age, mol_ne) %>%
slice_head(n = 10)
ne_top
# IQR-based outlier threshold
Q1_ne = quantile(avg_data$mol_ne, 0.25, na.rm = TRUE)
Q3_ne = quantile(avg_data$mol_ne, 0.75, na.rm = TRUE)
IQR_ne = Q3_ne - Q1_ne
upper_bound_ne = Q3_ne + 1.5 * IQR_ne
upper_bound_ne
avg_data_ne = avg_data %>%
mutate(ne_outlier = mol_ne > upper_bound_ne)
# Full-data test
wilcox.test(mol_ne ~ group, data = avg_data_ne)
# Outlier-removed test
avg_data_ne_no_outliers = avg_data_ne %>%
filter(!ne_outlier)
wilcox.test(mol_ne ~ group, data = avg_data_ne_no_outliers)
# Mean and median by group
avg_data_ne %>%
group_by(group) %>%
summarise(
mean_ne = mean(mol_ne, na.rm = TRUE),
median_ne = median(mol_ne, na.rm = TRUE),
min_ne = min(mol_ne, na.rm = TRUE),
max_ne = max(mol_ne, na.rm = TRUE),
.groups = "drop"
)
```
```{r univar_rf}
rf_base_data = avg_data %>%
mutate(group = factor(group, levels = c("NAC", "ASD"))) # ASD is the positive class
set.seed(123)
folds = vfold_cv(rf_base_data, v = 5, strata = group)
mol_names = names(rf_base_data)[str_starts(names(rf_base_data), "mol_")]
get_cv_rf_auc = function(mol) {
df = rf_base_data %>% select(group, !!sym(mol))
rec = recipe(group ~ ., data = df) %>%
step_zv(all_predictors()) # Remove zero-variance predictors (just in case)
rf_spec = rand_forest(trees = 500) %>%
set_mode("classification") %>%
set_engine("ranger", importance = "none")
wf = workflow() %>%
add_recipe(rec) %>%
add_model(rf_spec)
res = fit_resamples(
wf,
resamples = folds,
metrics = metric_set(roc_auc),
control = control_resamples(save_pred = TRUE, verbose = FALSE)
)
auc_val = collect_metrics(res) %>%
filter(.metric == "roc_auc") %>%
pull(mean)
tibble(molecule = mol, auc = auc_val)
}
rf_auc_all = map_dfr(mol_names, get_cv_rf_auc) %>%
arrange(desc(auc))
# rf_auc_all %>%
# mutate(molecule = fct_reorder(molecule, auc)) %>%
# ggplot(aes(x = molecule, y = auc)) +
# geom_col(fill = "#21918c") +
# coord_flip() +
# labs(
# title = "Random Forest Classification Performance by Molecule",
# y = "Cross-Validated AUC",
# x = NULL
# ) +
# geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray")
ind_auc_results = map_dfr(mol_names, function(mol) {
df = avg_data %>% select(group, !!sym(mol))
rec = recipe(group ~ ., data = df) %>%
step_zv(all_predictors())
wf = workflow() %>%
add_recipe(rec) %>%
add_model(rand_forest(trees = 500) %>%
set_mode("classification") %>%
set_engine("ranger"))
res = fit_resamples(
wf, resamples = folds,
metrics = metric_set(roc_auc),
control = control_resamples(save_pred = TRUE)
)
preds = collect_predictions(res)
roc_obj = roc(response = preds$group, predictor = preds$.pred_ASD,
levels = c("NAC", "ASD"), direction = "<", ci = TRUE, boot.n = 2000)
tibble(
metabolite = mol,
auc = as.numeric(auc(roc_obj)),
ci_low = as.numeric(ci.auc(roc_obj)[1]),
ci_high = as.numeric(ci.auc(roc_obj)[3])
)
}) %>%
arrange(desc(auc)) %>%
mutate(
auc_ci = sprintf("%.3f [%.3f-%.3f]", auc, ci_low, ci_high)
) %>%
select(metabolite, auc_ci) %>%
left_join(name_map %>% select(mol_name, compound_name), by = c("metabolite" = "mol_name")) %>%
mutate(metabolite = compound_name) %>%
select(-compound_name)
# Table S3
ind_auc_results %>%
left_join(name_map %>% select(compound_acronym, compound_name), by = c("metabolite" = "compound_name")) %>%
mutate(metabolite = str_c(metabolite, " (", compound_acronym, ")")) %>%
select(-compound_acronym) %>%
write_csv("output_data/s3_univar_rf_auc.csv")
```
```{r univar_roc}
top_molecule = rf_auc_all %>%
slice_max(auc, n = 3) %>%
pull(molecule)
rf_top_data = avg_data %>%
mutate(group = factor(group, levels = c("NAC", "ASD"))) %>%
select(group, all_of(top_molecule))
rec = recipe(group ~ ., data = rf_top_data) %>%
step_zv(all_predictors())
rf_spec = rand_forest(trees = 500) %>%
set_mode("classification") %>%
set_engine("ranger")
wf = workflow() %>%
add_recipe(rec) %>%
add_model(rf_spec)
set.seed(123)
folds = vfold_cv(rf_top_data, v = 5, strata = group)
res = fit_resamples(
wf,
resamples = folds,
metrics = metric_set(roc_auc),
control = control_resamples(save_pred = TRUE)
)
preds = collect_predictions(res)
roc_obj = roc(
response = preds$group,
predictor = preds$.pred_ASD,
levels = c("NAC", "ASD"),
direction = "<",
ci = TRUE,
ci.alpha = 0.95,
boot.n = 2000,
stratified = FALSE
)
roc_df = tibble(
specificity = rev(roc_obj$specificities),
sensitivity = rev(roc_obj$sensitivities)
)
auc_val = auc(roc_obj)
ci_vals = ci.auc(roc_obj)
# ggplot(roc_df, aes(x = 1 - specificity, y = sensitivity)) +
# geom_line(color = "#21918c", linewidth = 1.2) +
# geom_abline(linetype = "dashed", color = "gray") +
# annotate(
# "text",
# x = 0.25,
# y = 0.25,
# label = sprintf("AUC = %.3f [%.3f–%.3f]", auc_val, ci_vals[1], ci_vals[3]),
# hjust = 0,
# size = 5
# ) +
# labs(
# title = paste0("ROC Curve with 95% CI for ", paste0(top_molecule, collapse = ", ")),
# x = "1 - Specificity",
# y = "Sensitivity"
# )
```
```{r multivar_rf_all, fig.dpi=300, fig.height=10}
set.seed(123)
folds = vfold_cv(avg_data, v = 5, strata = group)
top_mols = rf_auc_all %>%
slice_max(auc, n = 3) %>%
pull(molecule)
model_definitions = tribble(
~model_name, ~formula, ~panel_group, ~color,
"Age", group ~ age, "Demographics", "#21908CFF",
"Sex", group ~ sex, "Demographics", "#440154FF",
"Age + Sex", group ~ age + sex, "Demographics", "#FDE725FF",
"BH4", group ~ mol_bh4, "Metabolites", "#35B779FF",
"KYN", group ~ mol_kyn, "Metabolites", "#31688EFF",
"GABA", group ~ mol_gaba, "Metabolites", "#440154FF",
"BH4 + KYN + GABA", reformulate(top_mols, response = "group"), "Metabolites", "#26828EFF",
"Age + BH4 + KYN + GABA", reformulate(c("age", top_mols), response = "group"), "Combined", "#5DC863FF",
"Sex + BH4 + KYN + GABA", reformulate(c("sex", top_mols), response = "group"), "Combined", "#3B528BFF",
"Age + Sex + BH4 + KYN + GABA", reformulate(c("age", "sex", top_mols), response = "group"), "Combined", "#481567FF"
)
get_rf_roc = function(formula, model_label) {
rec = recipe(formula, data = avg_data) %>%
step_zv(all_predictors()) %>%
step_dummy(all_nominal_predictors())
wf = workflow() %>%
add_model(rand_forest(trees = 500, mode = "classification") %>%
set_engine("ranger")) %>%
add_recipe(rec)
res = fit_resamples(
wf, resamples = folds,
metrics = metric_set(roc_auc),
control = control_resamples(save_pred = TRUE)
)
preds = collect_predictions(res)
roc_obj = roc(
response = preds$group,
predictor = preds$.pred_ASD,
levels = c("NAC", "ASD"),
direction = "<",
ci = TRUE, boot.n = 2000
)
tibble(
model_name = model_label,
specificity = rev(roc_obj$specificities),
sensitivity = rev(roc_obj$sensitivities),
auc = as.numeric(auc(roc_obj)),
ci_low = as.numeric(ci.auc(roc_obj)[1]),
ci_high = as.numeric(ci.auc(roc_obj)[3])
)
}
roc_results = model_definitions %>%
mutate(roc = map2(formula, model_name, get_rf_roc)) %>%
unnest(roc, names_sep = "_") %>%
rename(
specificity = roc_specificity,
sensitivity = roc_sensitivity,
auc = roc_auc,
ci_low = roc_ci_low,
ci_high = roc_ci_high
)
roc_labeled = roc_results %>%
group_by(model_name) %>%
mutate(
label = sprintf("%s\nAUC = %.3f [%.3f-%.3f]", model_name, first(auc), first(ci_low), first(ci_high))
) %>%
ungroup()
label_order = roc_labeled %>%
distinct(model_name, label) %>%
arrange(factor(model_name, levels = model_definitions$model_name)) %>%
pull(label)
roc_labeled = roc_labeled %>%
mutate(label = factor(label, levels = label_order))
plot_roc_group = function(data, panel) {
df = filter(data, panel_group == panel)
ggplot(df, aes(x = 1 - specificity, y = sensitivity, color = label)) +
geom_line(linewidth = 1) +
geom_abline(linetype = "dashed", color = "gray") +
scale_color_manual(values = setNames(df$color, df$label)) +
labs(
x = "1-Specificity",
y = "Sensitivity",
color = NULL
) +
coord_equal()
}
p_demo = plot_roc_group(roc_labeled, "Demographics")
p_metab = plot_roc_group(roc_labeled, "Metabolites")
p_combined = plot_roc_group(roc_labeled, "Combined")
# Figure 3
p_demo / p_metab / p_combined +
plot_annotation(tag_levels = "A", tag_prefix = "(", tag_suffix = ")")
```
```{r multivar_rf_metrics}
set.seed(123)
folds = vfold_cv(avg_data, v = 5, strata = group)
compute_model_metrics = function(formula, label) {
rec = recipe(formula, data = avg_data) %>%
step_zv(all_predictors()) %>%
step_dummy(all_nominal_predictors())
rf_spec = rand_forest(trees = 500) %>%
set_mode("classification") %>%
set_engine("ranger")
wf = workflow() %>%
add_model(rf_spec) %>%
add_recipe(rec)
res = fit_resamples(
wf,
resamples = folds,
metrics = metric_set(roc_auc, accuracy, sens, spec, precision, recall, f_meas),
control = control_resamples(save_pred = TRUE)
)
preds = collect_predictions(res)
roc_obj = roc(
response = preds$group,
predictor = preds$.pred_ASD,
levels = c("NAC", "ASD"),
direction = "<",
ci = TRUE,
boot.n = 2000
)
auc_val = as.numeric(auc(roc_obj))
ci_vals = as.numeric(ci.auc(roc_obj))
metric_vals = preds %>%
summarise(
accuracy = accuracy_vec(group, .pred_class),
sens = sens_vec(group, .pred_class, event_level = "second"),
spec = spec_vec(group, .pred_class, event_level = "second"),
precision = precision_vec(group, .pred_class, event_level = "second"),
recall = recall_vec(group, .pred_class, event_level = "second"),
f_meas = f_meas_vec(group, .pred_class, event_level = "second")
)
tibble(
model = label,
auc = auc_val,
auc_ci_low = ci_vals[1],
auc_ci_high = ci_vals[3]
) %>%
bind_cols(metric_vals)
}
formulas = list(
"Age" = group ~ age,
"Sex" = group ~ sex,
"Age + Sex" = group ~ age + sex,
"BH4" = group ~ mol_bh4,
"KYN" = group ~ mol_kyn,
"GABA" = group ~ mol_gaba,
"BH4 + KYN + GABA" = reformulate(top_mols, response = "group"),
"Sex + BH4 + KYN + GABA" = reformulate(c("sex", top_mols), response = "group"),
"Age + BH4 + KYN + GABA" = reformulate(c("age", top_mols), response = "group"),
"Age + Sex + BH4 + KYN + GABA" = reformulate(c("age", "sex", top_mols), response = "group")
)
metrics_all = imap_dfr(formulas, compute_model_metrics) %>%
mutate(across(where(is.numeric), ~ round(.x, 3)),
auc = paste0(auc, " [", auc_ci_low, "-", auc_ci_high, "]")) %>%
select(-auc_ci_low, -auc_ci_high)
```
```{r rf_model_comparison}
# Model comparison to assess whether adding more metabolites
# improves classification performance beyond the top 3 features.
set.seed(123)
folds = vfold_cv(avg_data, v = 5, strata = group)
# Select metabolite sets
top3_mols = rf_auc_all %>%
slice_max(auc, n = 3, with_ties = FALSE) %>%
pull(molecule)
top5_mols = rf_auc_all %>%
slice_max(auc, n = 5, with_ties = FALSE) %>%
pull(molecule)
all_mols = names(avg_data)[str_starts(names(avg_data), "mol_")]
feature_labels = list(
"Top 3 metabolites" = str_c(top3_mols, collapse = ", "),
"Top 3 + Trp" = str_c(c(top3_mols, "mol_trp"), collapse = ", "),
"Top 5 metabolites" = str_c(top5_mols, collapse = ", "),
"All 18 metabolites" = "All metabolites"
)
# Reuse the same modeling framework as in Table S4
compute_model_metrics_compare = function(formula, label) {
rec = recipe(formula, data = avg_data) %>%
step_zv(all_predictors()) %>%
step_dummy(all_nominal_predictors())
rf_spec = rand_forest(trees = 500) %>%
set_mode("classification") %>%
set_engine("ranger")
wf = workflow() %>%
add_model(rf_spec) %>%
add_recipe(rec)
res = fit_resamples(
wf,
resamples = folds,
metrics = metric_set(roc_auc),
control = control_resamples(save_pred = TRUE)
)
preds = collect_predictions(res) %>%
mutate(
.pred_class = factor(
if_else(.pred_ASD >= 0.5, "ASD", "NAC"),
levels = c("NAC", "ASD")
)
)
roc_obj = roc(
response = preds$group,
predictor = preds$.pred_ASD,
levels = c("NAC", "ASD"),
direction = "<",
ci = TRUE,
boot.n = 2000
)
auc_val = as.numeric(auc(roc_obj))
ci_vals = as.numeric(ci.auc(roc_obj))
metric_vals = preds %>%
summarise(
accuracy = accuracy_vec(group, .pred_class),
sens = sens_vec(group, .pred_class, event_level = "second"),
spec = spec_vec(group, .pred_class, event_level = "second"),
precision = precision_vec(group, .pred_class, event_level = "second"),
f_meas = f_meas_vec(group, .pred_class, event_level = "second")
)
tibble(
model = label,
auc = auc_val,
auc_ci_low = ci_vals[1],
auc_ci_high = ci_vals[3]
) %>%
bind_cols(metric_vals)
}
# Define comparison models
compare_formulas = list(
"Top 3 metabolites" = reformulate(top3_mols, response = "group"),
"Top 3 + Trp" = reformulate(c(top3_mols, "mol_trp"), response = "group"),
"Top 5 metabolites" = reformulate(top5_mols, response = "group"),
"All 18 metabolites" = reformulate(all_mols, response = "group")
)
# Run comparison
rf_model_compare = imap_dfr(compare_formulas, compute_model_metrics_compare) %>%
mutate(
features = feature_labels[model] %>% unlist()
) %>%
mutate(
across(where(is.numeric), \(x) round(x, 3)),
auc = paste0(auc, " [", auc_ci_low, "-", auc_ci_high, "]")
) %>%
select(model, features, everything(), -auc_ci_low, -auc_ci_high)
rf_model_compare
# Optional export for supplement
rf_model_compare %>%
write_csv("output_data/s6_rf_model_comparison.csv")
```
```{r multivar_rf_calibration, fig.dpi=300, fig.height=10}
get_calibration_metrics = function(formula, model_label) {
rec = recipe(formula, data = avg_data) %>%
step_zv(all_predictors()) %>%
step_dummy(all_nominal_predictors())
wf = workflow() %>%
add_model(rand_forest(trees = 500, mode = "classification") %>% set_engine("ranger")) %>%
add_recipe(rec)
res = fit_resamples(
wf, resamples = folds,
control = control_resamples(save_pred = TRUE)
)
preds = collect_predictions(res) %>%
mutate(group = factor(group, levels = c("NAC", "ASD"))) # ASD = positive class
brier = mean((if_else(preds$group == "ASD", 1, 0) - preds$.pred_ASD)^2)
cal_curve_raw = preds %>%
transmute(
pred = .pred_ASD,
obs = as.integer(group == "ASD"),
model_name = model_label
)
list(
brier = tibble(model_name = model_label, brier_score = brier),
curve_raw = cal_curve_raw
)
}
set.seed(123)
cal_results = model_definitions %>%
mutate(cal = map2(formula, model_name, get_calibration_metrics))
brier_all = map_dfr(cal_results$cal, "brier")
curve_all_raw = map_dfr(cal_results$cal, "curve_raw")
curve_all_annotated = curve_all_raw %>%
left_join(model_definitions %>% select(model_name, panel_group, color), by = "model_name") %>%
left_join(brier_all, by = "model_name") %>%
mutate(
label = sprintf("%s\nBrier = %.3f", model_name, brier_score)
)
label_order_brier = model_definitions %>%
left_join(brier_all, by = "model_name") %>%
mutate(label = sprintf("%s\nBrier = %.3f", model_name, brier_score)) %>%
pull(label)
curve_all_annotated = curve_all_annotated %>%
mutate(label = factor(label, levels = label_order_brier))
plot_calibration_panel = function(data, panel_name) {
df = filter(data, panel_group == panel_name)
ggplot(df, aes(x = pred, y = obs, color = label)) +
geom_smooth(method = "loess", formula = y ~ x, se = FALSE, span = 1, linewidth = 1) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
scale_color_manual(values = setNames(df$color, df$label)) +
labs(
x = "Predicted Probability (ASD)",
y = "Observed Proportion (ASD)",
color = NULL
) +
coord_fixed(xlim = c(0, 1), ylim = c(0, 1))
}
p_cal_demo = plot_calibration_panel(curve_all_annotated, "Demographics")
p_cal_metab = plot_calibration_panel(curve_all_annotated, "Metabolites")
p_cal_combined = plot_calibration_panel(curve_all_annotated, "Combined")
# Figure 4
p_cal_demo / p_cal_metab / p_cal_combined +
plot_annotation(tag_levels = "A", tag_prefix = "(", tag_suffix = ")")
# Table S4
metrics_all %>%
left_join(brier_all %>% mutate(brier_score = signif(brier_score, 3)), by = c("model" = "model_name")) %>%
write_csv("output_data/s4_multivar_rf_perf_brier_scores.csv")
```
```{r session-info}
sessionInfo()
```

statistical_analysis.Rmd at commit f07a3a4, under MIT · at the source

Overview

Authors: Kevin Liu¹, Huixi Li², Shaohan Zhang³, Muya Xi¹, Junru Zhu^2,4, Jonathan Chen³, William Xu³, Alexander Xie¹, Alexandros Makriyannis², Jason J. Guo^2,4, Xue-Jun Kong^1,3,5

Fetal-Neonatal Neuroimaging and Developmental Science Center, Division of Newborn Medicine, Boston Children’s Hospital, Boston, MA, United States
Department of Pharmaceutical Sciences, Center for Drug Discovery, Northeastern University, Boston, MA, United States
Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charlestown, MA, United States
Department of Chemistry and Chemical Biology, Barnett Institute for Chemical and Biological Analysis, Northeastern University, Boston, MA, United States
Department of Medicine and Psychiatry, Beth Israel Deaconess Medical Center, Boston, MA, United States

Institutions: Boston Children's Hospital (United States); Northeastern University (United States); Massachusetts General Hospital (United States); Athinoula A. Martinos Center for Biomedical Imaging (United States); Beth Israel Deaconess Medical Center (United States)

Journal: Frontiers in neuroscience, volume 20, article 1858005

Dates: received 16 April 2026; accepted 22 June 2026; published online 3 July 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.3389/fnins.2026.1858005 · PMID 42490967 · PMCID PMC13375968 · OpenAlex W7167221654

Open access: gold, a free copy (OpenAlex)

Status: code verified (how the code is found and verified)

Categories: genetics / omics (modality), autism (population), clinical / translational (subfield)

Methods: Connectivity, Statistics, Machine learning

Keywords: autism spectrum disorder, biomarkers, catecholamines, gut-brain axis, machine learning, microbiome, stool metabolomics, tryptophan metabolism

Topic: Autism Spectrum Disorder Research (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: Boston Children's Hospital (92436); MGH (233263)

Citations: not cited yet (Europe PMC); 62 references in the paper

Abstract

Background: Gut-brain axis dysregulation and microbiome-linked metabolic alterations have been implicated in autism spectrum disorder (ASD), but the contribution of gut-derived neuroactive metabolites remains incompletely characterized.

Methods: We conducted a cross-sectional case-control study of 59 participants (32 ASD, 27 controls) and quantified 18 stool metabolites related to catecholamine synthesis, inhibitory neurotransmission, and tryptophan-linked NAD+-precursor metabolism using targeted liquid chromatography-tandem mass spectrometry. Group differences were assessed using fold-change analysis and linear models adjusted for age and sex. Random forest models evaluated classification performance, and within-group Spearman correlations were used to examine metabolic relationships.

Results: Norepinephrine showed the largest increase in ASD, whereas dopamine and tetrahydrobiopterin exhibited nominal group differences that did not remain significant after correction for multiple testing. A three-metabolite panel comprising tetrahydrobiopterin, γ-aminobutyric acid, and kynurenine showed exploratory discrimination between groups (area under the receiver operating characteristic curve = 0.750, 95% confidence interval 0.622–0.878), but this performance requires external validation. Correlation analysis revealed conserved bile acid coupling in both groups. In controls, tryptophan was positively associated with kynurenine, whereas this relationship was not observed in ASD. Instead, ASD samples showed broader associations between tryptophan and metabolites linked to neurotransmission and NAD+-precursor metabolism.

Conclusion: Stool metabolite profiling revealed altered organization of tryptophan- and catecholamine-linked metabolic associations in ASD and identified a small metabolite panel with exploratory discriminative potential. These findings provide a foundation for future studies examining gut-derived neuroactive metabolites in ASD and their relationship to gut-brain axis biology.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

kevinliu-bmb/targeted-asd-stool-metabolomics

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: f07a3a45dec086833da3a29b2dcc84d758b69361, 4 July 2026

Languages: R (1)

Size: 4 files, 1 script

Software Heritage: not archived

Found in: the text, “Statistical analysis”

Holds: README, license file, 1 notebook

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: broom (1 file), ggpubr (1 file), patchwork (1 file), pROC (1 file), rstatix (1 file), tidyverse (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

3 files

statistical_analysis.Rmd, R, 1,040 lines, 4 matches
LICENSE, License, 21 lines
README.md, Text, 4 lines

At the source: github.com/kevinliu-bmb/targeted-asd-stool-metabolomics

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 11 authors, 8 keywords, 2 funders, 61 references.

Cite

This paper

Liu, K., Li, H., Zhang, S., Xi, M., Zhu, J., Chen, J., Xu, W., Xie, A., Makriyannis, A., Guo, J. J., & Kong, X.-J. (2026). Targeted stool metabolomics suggests exploratory catecholamine- and tryptophan-linked metabolic features in autism spectrum disorder. Frontiers in neuroscience, 20, 1858005. https://doi.org/10.3389/fnins.2026.1858005

BibTeX

@article{liu2026targeted, author = {Liu, Kevin and Li, Huixi and Zhang, Shaohan and Xi, Muya and Zhu, Junru and Chen, Jonathan and Xu, William and Xie, Alexander and Makriyannis, Alexandros and Guo, Jason J. and Kong, Xue-Jun}, title = {{Targeted stool metabolomics suggests exploratory catecholamine- and tryptophan-linked metabolic features in autism spectrum disorder}}, journal = {Frontiers in neuroscience}, year = {2026}, month = jul, volume = {20}, pages = {1858005}, publisher = {Frontiers Media SA}, issn = {1662-4548}, doi = {10.3389/fnins.2026.1858005}, url = {https://doi.org/10.3389/fnins.2026.1858005}, pmid = {42490967}, pmcid = {PMC13375968} }

RIS

TY - JOUR AU - Liu, Kevin AU - Li, Huixi AU - Zhang, Shaohan AU - Xi, Muya AU - Zhu, Junru AU - Chen, Jonathan AU - Xu, William AU - Xie, Alexander AU - Makriyannis, Alexandros AU - Guo, Jason J. AU - Kong, Xue-Jun TI - Targeted stool metabolomics suggests exploratory catecholamine- and tryptophan-linked metabolic features in autism spectrum disorder T2 - Frontiers in neuroscience J2 - Front Neurosci PY - 2026 DA - 2026/07/03 VL - 20 SP - 1858005 SN - 1662-4548 PB - Frontiers Media SA DO - 10.3389/fnins.2026.1858005 UR - https://doi.org/10.3389/fnins.2026.1858005 LA - en ER -

CSL-JSON

{ "id": "10.3389/fnins.2026.1858005", "type": "article-journal", "title": "Targeted stool metabolomics suggests exploratory catecholamine- and tryptophan-linked metabolic features in autism spectrum disorder", "container-title": "Frontiers in neuroscience", "author": [ { "family": "Liu", "given": "Kevin" }, { "family": "Li", "given": "Huixi" }, { "family": "Zhang", "given": "Shaohan" }, { "family": "Xi", "given": "Muya" }, { "family": "Zhu", "given": "Junru" }, { "family": "Chen", "given": "Jonathan" }, { "family": "Xu", "given": "William" }, { "family": "Xie", "given": "Alexander" }, { "family": "Makriyannis", "given": "Alexandros" }, { "family": "Guo", "given": "Jason J." }, { "family": "Kong", "given": "Xue-Jun" } ], "container-title-short": "Front Neurosci", "volume": "20", "page": "1858005", "DOI": "10.3389/fnins.2026.1858005", "PMID": "42490967", "PMCID": "PMC13375968", "ISSN": "1662-4548", "publisher": "Frontiers Media SA", "URL": "https://doi.org/10.3389/fnins.2026.1858005", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 3 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-025-58459-1 [code]: Relationships between brain activity, tryptophan-related gut metabolites, and autism symptomatology.
Journal: Nature communications
In common: autism, genetics / omics, 7 references
[2] doi:10.1016/j.xcrm.2026.102655 [code]: An integrative multi-omics approach identifies microbiome alterations linked to pathological and behavioral features in autism spectrum disorder.
Journal: Cell reports. Medicine
In common: pROC, patchwork, ggpubr, 1 other tool, autism, genetics / omics, 2 references
[3] doi:10.1126/sciadv.aec9291 [code]: Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks.
Journal: Science advances
In common: pROC, broom, rstatix, 3 other tools, autism
[4] doi:10.1038/s41467-025-56829-3 [code]: Machine learning-based meta-analysis reveals gut microbiome alterations associated with Parkinson's disease.
Journal: Nature communications
In common: pROC, broom, rstatix, 3 other tools, clinical / translational
[5] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: pROC, broom, rstatix, 3 other tools, genetics / omics
[6] doi:10.1186/s12888-026-08178-8 [code]: Machine learning model to identify gut microbiome-derived metabolites as potential biomarkers of autism spectrum disorder: a pilot study.
Journal: BMC psychiatry
In common: autism, clinical / translational, genetics / omics, 4 references
[7] doi:10.1111/jsr.70309 [code]: An Electroencephalographic Study of Sleep Spindle and Infraslow Oscillation in Children With Autism Spectrum Disorder.
Journal: Journal of sleep research
In common: broom, patchwork, ggpubr, 1 other tool, autism, 2 references
[8] doi:10.1038/s41467-025-64328-8 [code]: Integrative brain omics approach highlights sn-1 lysophosphatidylethanolamine in Alzheimer's dementia.
Journal: Nature communications
In common: pROC, broom, rstatix, 2 other tools, clinical / translational, genetics / omics
[9] doi:10.1038/s43587-025-01016-8 [code]: Organ-specific proteomic aging clocks predict disease and longevity across diverse populations.
Journal: Nature aging
In common: pROC, broom, patchwork, 2 other tools, clinical / translational, genetics / omics
[10] doi:10.1158/2159-8290.cd-25-0074 [code]: Detection of Brain Cancer Using Genome-wide Cell-free DNA Fragmentomes.
Journal: Cancer discovery
In common: pROC, broom, patchwork, 2 other tools, clinical / translational, genetics / omics

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown

<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.3389/fnins.2026.1858005.svg)](https://openscicode.org/badge/doi/10.3389/fnins.2026.1858005)
<!-- oscr-badge:end -->

reStructuredText

.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.3389/fnins.2026.1858005.svg
   :target: https://openscicode.org/badge/doi/10.3389/fnins.2026.1858005
   :alt: OSCR

.. oscr-badge:end

HTML

<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.3389/fnins.2026.1858005"><img src="https://openscicode.org/badge/doi/10.3389/fnins.2026.1858005.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

kevinliu-bmb/targeted-asd-stool-metabolomics: README.md on GitHub

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.