GMSB文章六：微生物SCFA關聯分析

歡迎大家關注全網生信學習者系列：

WX公zhong號：生信學習者
Xiao hong書：生信學習者
知hu：生信學習者
CDSN：生信學習者2

介紹

微生物短鏈脂肪酸（SCFAs）是由腸道微生物發酵膳食纖維、抗性淀粉、低聚糖等碳水化合物產生的一類有機脂肪酸，主要包括乙酸、丙酸和丁酸。SCFAs在腸道健康、免疫調節、能量代謝等方面發揮重要作用，并且與多種疾病狀態有關。

通過約束線性混合效應模型（CLME）研究SCFA與排序暴露分組之間的單調遞增關系，這種方法通常用于分析具有層次結構的數據，能夠處理數據中的非獨立性和非正態分布問題。在研究SCFA與某些因素（如飲食、疾病狀態等）的關系時，CLME可以用來評估這些因素對SCFA水平的影響，并檢測它們之間是否存在單調遞增或遞減的關系。

加載R包

library(readr)
library(tidyverse) 
library(CLME)
library(magrittr)
library(qwraps2)
library(ggprism)
library(ggsci)
library(ggpubr)
library(rstatix)
library(jtools)
library(kableExtra)

導入數據

大家通過以下鏈接下載數據：

百度網盤鏈接：https://pan.baidu.com/s/1fz5tWy4mpJ7nd6C260abtg
提取碼: 請關注WX公zhong號_生信學習者_后臺發送 復現gmsb 獲取提取碼

df_merge <- read_csv("./data/GMSB-data/df_merge.csv", show_col_types = FALSE)head(df_merge[, 1:6])

sampleid	subjid	visit	visit_num	status	age
F-1	DUMMY1	v1	10	nc	58
F-2	DUMMY2	v1	10	nc	60
F-3	DUMMY3	v1	10	nc	58
F-4	DUMMY4	v1	10	nc	72
F-5	DUMMY5	v1	10	nc	57
F-6	DUMMY6	v1	10	nc	59

數據預處理

對不同visit進行分組中，提取它們的SCFA短鏈脂肪酸

df_v1 <- df_merge %>%dplyr::filter(visit == "v1",group1 != "missing")
df_v2 <- df_merge %>%dplyr::filter(visit == "v2",group1 != "missing")df_v1$group1 <- recode(df_v1$group1,`g1` = "G1", `g2` = "G2",`g3` = "G3", `g4` = "G4")
df_v1$group2 <- recode(df_v1$group2,`g1` = "G1", `g2` = "G2",`g3` = "G3", `g4` = "G4", `g5` = "G5")df_v1 <- df_v1 %>%dplyr::transmute(subjid, status, leu3p, leu2p, recept_anal, abx_use, group1, group2, acetate_v1 = acetate, propionate_v1 = propionate, butyrate_v1 = butyrate, valerate_v1 = valerate)
df_v2 <- df_v2 %>%dplyr::transmute(subjid, acetate_v2 = acetate, propionate_v2 = propionate, butyrate_v2 = butyrate, valerate_v2 = valerate) df_v1 <- df_v1 %>%dplyr::left_join(df_v2, by = "subjid")df_v1$group1 <- factor(df_v1$group1, levels = c("G1", "G2", "G3", "G4"), ordered = TRUE)
df_v1$group2 <- factor(df_v1$group2, levels = c("G1", "G2", "G3", "G4", "G5"), ordered = TRUE)head(df_v1[, 1:6])

subjid	status	leu3p	leu2p	recept_anal	abx_use
DUMMY1	nc	40.5	35.9	5	yes
DUMMY2	nc	53.9	27.7	6	yes
DUMMY3	nc	30.9	42.5	3	no
DUMMY4	nc	50.3	19.3	7	no
DUMMY5	nc	42.8	38.8	4	no
DUMMY6	nc	36.9	46.7	4	yes

visit1

visit1群體的SCFA的數據分布情況

df_v1 %>% dplyr::select(acetate_v1, propionate_v1, butyrate_v1, valerate_v1) %>%pastecs::stat.desc() %>% kable() %>% kable_styling()

在這里插入圖片描述

結果：不同的SCFA的數據描述結果，比如acetate_v1最大和最小值都要高于其他的。

visit2

visit2群體的SCFA的數據分布情況

df_v1 %>% dplyr::select(acetate_v2, propionate_v2, butyrate_v2, valerate_v2) %>%pastecs::stat.desc() %>% kable() %>% kable_styling()

在這里插入圖片描述

結果：不同的SCFA的數據描述結果，比如acetate_v1最大和最小值都要高于其他的。

SCFA ~ groups: box plots

線性回歸校正abx_use因素比較不同分組的差異

plot_box <- function(group_lab, group_var, key_var, title, y.position, step.increase) {df_fig <- df_v1 %>%transmute(group = get(group_var), value = get(key_var), abx_use)df_fig$group <- factor(df_fig$group, ordered = FALSE)lm_fit <- lm(value ~ group + abx_use, data = df_fig)df_p <- data.frame(group1 = group_lab[1],group2 = group_lab[-1],p = summary(lm_fit)$coef[grepl("group", names(coef(lm_fit))), "Pr(>|t|)"]) %>%mutate(p = round(p, 2))bxp <- ggboxplot(df_fig, x = "group", y = "value", color = "group",add = "jitter", xlab = FALSE, ylab = FALSE, title = title) +stat_pvalue_manual(df_p, y.position = y.position, step.increase = step.increase, label = "p") +scale_color_npg(name = "Group") +theme(axis.text.x = element_blank(),axis.ticks.x = element_blank(),strip.background = element_rect(fill = "white"),legend.position = "bottom",plot.title = element_text(hjust = 0.5))return(bxp)
}

Primary group

Primary group: 按照頻率分組

G1: # receptive anal intercourse = 0
G2: # receptive anal intercourse = 1
G3: # receptive anal intercourse = 2 - 5
G4: # receptive anal intercourse = 6 +

bxp_list <- list()
var_list <- c("acetate_v1", "butyrate_v1", "propionate_v1", "valerate_v1")
title_list <- c("Acetate", "Butyrate", "Propionate", "Valerate")
y_pos_list <- c(6, 0.25, 0.35, 0.07)
step_list <- c(0.1, 0.1, 0.1, 0.1)for (i in seq_along(var_list)) {bxp <- plot_box(group_lab = c("G1", "G2", "G3", "G4"),group_var = "group1",key_var = var_list[i], title = title_list[i],y.position = y_pos_list[i], step.increase = step_list[i])bxp_list[[i]] <- bxp
}ggarrange(bxp_list[[1]], bxp_list[[2]], bxp_list[[3]], bxp_list[[4]],common.legend = TRUE)

在這里插入圖片描述

結果：上述短鏈脂肪酸在分組間均不存在顯著差異。

Secondary group

Secondary group: 按照頻率分組2

G1: # receptive anal intercourse = 0
G2: # receptive anal intercourse = 1
G3: # receptive anal intercourse = 2 - 3
G4: # receptive anal intercourse = 4 - 8
G5: # receptive anal intercourse = 9 +

bxp_list <- list()
var_list <- c("acetate_v1", "butyrate_v1", "propionate_v1", "valerate_v1")
title_list <- c("Acetate", "Butyrate", "Propionate", "Valerate")
y_pos_list <- c(6, 0.25, 0.35, 0.07)
step_list <- c(0.1, 0.1, 0.1, 0.1)for (i in seq_along(var_list)) {bxp <- plot_box(group_lab = c("G1", "G2", "G3", "G4", "G5"),group_var = "group2",key_var = var_list[i], title = title_list[i],y.position = y_pos_list[i], step.increase = step_list[i])bxp_list[[i]] <- bxp
}ggarrange(bxp_list[[1]], bxp_list[[2]], bxp_list[[3]], bxp_list[[4]],common.legend = TRUE)

在這里插入圖片描述

結果：短鏈脂肪酸在分組間的分布。

Acetate在G1和G5間顯著差異；
Valerate在G1和G5間顯著差異

SCFA ~ groups: increasing trend

約束線性混合效應模型（CLME, Constrained Linear Mixed Effects Models）：用于評估性暴露組和血漿炎癥細胞因子水平之間的單調遞增趨勢，同時校正細菌抗生素使用情況。

plot_clme用于繪制趨勢圖

plot_clme <- function(clme_obj, group, y_min, y_max, p_gap, ann_pos, trend = "increase", ...) {overall_p <- clme_obj$p.valueind_p <- clme_obj$p.value.indest_mean <- clme_obj$theta[group]est_se <- sqrt(diag(clme_obj$cov.theta))[group]df_fig <- data.frame(x = group, y = est_mean, err = est_se)if (est_mean[2] < est_mean[length(est_mean)]) {start_p_pos <- est_mean[2] + p_gapend_p_pos <- max(est_mean) + p_gap} else if (est_mean[2] > est_mean[length(est_mean)]) {start_p_pos <- max(est_mean) + p_gapend_p_pos <- min(est_mean) + p_gap} else {if (trend == "increase") {start_p_pos <- est_mean[2] + p_gapend_p_pos <- max(est_mean) + 2 * p_gap} else {start_p_pos <- max(est_mean) + 2 * p_gapend_p_pos <- min(est_mean) + p_gap}}df_p <- data.frame(group1 = group[seq_len(length(group) - 1)],group2 = group[-1],x = group[-1],label = paste0("p = ", round(ind_p, 3)),y.position = seq.int(from = start_p_pos, to = end_p_pos, length.out = length(group) - 1))df_ann <- data.frame(x = group[1], y = ann_pos,fill = "white",label = paste0("Trend p-value = ", round(overall_p, 3)))fig <- df_fig %>%ggplot(aes(x = x, y = y)) +geom_bar(stat = "identity", color = "black", aes(fill = x)) + geom_errorbar(aes(ymin = y - 1.96 * err, ymax = y + 1.96 * err), width = .2, position = position_dodge(.9)) +add_pvalue(df_p,xmin = "group1",xmax = "group2",label = "label",y.position = "y.position",remove.bracket = FALSE, ...) +geom_label(data = df_ann, aes(x = x, y = y, label = label), size = 4, vjust = -0.5, hjust = 0, color = "black") +ylim(y_min, y_max) +theme_bw()return(fig)
}

P-value is obtained by linear regression adjusting for antibiotics usage.
P-values were not adjusted for multiple comparisons.

限制參數

cons <- list(order = "simple", decreasing = FALSE, node = 1)

Primary grouping: Acetate

短鏈脂肪酸Acetate與Primary grouping的相關關系

fit1 <- clme(acetate_v1 ~ group1 + abx_use, data = df_v1, constraints = cons, seed = 123)
summ_fit1 <- summary(fit1)fig_ace_primary <- plot_clme(summ_fit1, group = c("G1", "G2", "G3", "G4"), y_min = 0, y_max = 4, p_gap = 0.7, ann_pos = 3.5) +scale_fill_npg(name = NULL) +labs(x = NULL, y = "Acetate")fig_ace_primary

在這里插入圖片描述

結果：acetate在不同組有顯著差異的單調遞增的趨勢（G3和G4組存在顯著差異）。

Primary grouping: Butyrate

短鏈脂肪酸Butyrate與Primary grouping的相關關系

fit2 <- clme(butyrate_v1 ~ group1 + abx_use, data = df_v1, constraints = cons, seed = 123)
summ_fit2 <- summary(fit2)fig_but_primary <- plot_clme(summ_fit2, group = c("G1", "G2", "G3", "G4"), y_min = 0, y_max = 0.18, p_gap = 0.04, ann_pos = 0.16) +scale_fill_npg(name = NULL) +labs(x = NULL, y = "Butyrate")fig_but_primary

在這里插入圖片描述

Primary grouping: Propionate

短鏈脂肪酸Propionate與Primary grouping的相關關系

fit3 <- clme(propionate_v1 ~ group1 + abx_use, data = df_v1, constraints = cons, seed = 123)
summ_fit3 <- summary(fit3)fig_pro_primary <- plot_clme(summ_fit3, group = c("G1", "G2", "G3", "G4"),y_min = 0, y_max = 0.25, p_gap = 0.05, ann_pos = 0.22) +scale_fill_npg(name = NULL) +labs(x = NULL, y = "Propionate")fig_pro_primary

在這里插入圖片描述

Primary grouping: Valerate

短鏈脂肪酸Valerate與Primary grouping的相關關系

fit4 <- clme(valerate_v1 ~ group1 + abx_use, data = df_v1, constraints = cons, seed = 123)
summ_fit4 <- summary(fit4)fig_val_primary <- plot_clme(summ_fit4, group = c("G1", "G2", "G3", "G4"),y_min = 0, y_max = 0.053, p_gap = 0.009, ann_pos = 0.048) +scale_fill_npg(name = NULL) +labs(x = NULL, y = "Valerate")fig_val_primary

在這里插入圖片描述

Secondary grouping: Acetate

短鏈脂肪酸Acetate與Secondary grouping的相關關系

fit1 <- clme(acetate_v1 ~ group2 + abx_use, data = df_v1, constraints = cons, seed = 123)
summ_fit1 <- summary(fit1)fig_ace_secondary <- plot_clme(summ_fit1, group = c("G1", "G2", "G3", "G4", "G5"), y_min = 0, y_max = 4, p_gap = 0.7, ann_pos = 3.5) +scale_fill_npg(name = NULL) +labs(x = NULL, y = "Acetate")fig_ace_secondary

在這里插入圖片描述

Secondary grouping: Butyrate

短鏈脂肪酸Butyrate與Secondary grouping的相關關系

fit2 <- clme(butyrate_v1 ~ group2 + abx_use, data = df_v1, constraints = cons, seed = 123)
summ_fit2 <- summary(fit2)fig_but_secondary <- plot_clme(summ_fit2, group = c("G1", "G2", "G3", "G4", "G5"), y_min = 0, y_max = 0.18, p_gap = 0.04, ann_pos = 0.16) +scale_fill_npg(name = NULL) +labs(x = NULL, y = "Butyrate")fig_but_secondary

在這里插入圖片描述

Secondary grouping: Propionate

短鏈脂肪酸Propionate與Secondary grouping的相關關系

fit3 <- clme(propionate_v1 ~ group2 + abx_use, data = df_v1, constraints = cons, seed = 123)
summ_fit3 <- summary(fit3)fig_pro_secondary <- plot_clme(summ_fit3, group = c("G1", "G2", "G3", "G4", "G5"),y_min = 0, y_max = 0.25, p_gap = 0.05, ann_pos = 0.22) +scale_fill_npg(name = NULL) +labs(x = NULL, y = "Propionate")fig_pro_secondary

在這里插入圖片描述

Secondary grouping: Valerate

短鏈脂肪酸Valerate與Secondary grouping的相關關系

fit4 <- clme(valerate_v1 ~ group2 + abx_use, data = df_v1, constraints = cons, seed = 123)
summ_fit4 <- summary(fit4)fig_val_secondary <- plot_clme(summ_fit4, group = c("G1", "G2", "G3", "G4", "G5"),y_min = 0, y_max = 0.053, p_gap = 0.009, ann_pos = 0.048) +scale_fill_npg(name = NULL) +labs(x = NULL, y = "Valerate")fig_val_secondary

在這里插入圖片描述