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])
sampleidsubjidvisitvisit_numstatusage
F-1DUMMY1v110nc58
F-2DUMMY2v110nc60
F-3DUMMY3v110nc58
F-4DUMMY4v110nc72
F-5DUMMY5v110nc57
F-6DUMMY6v110nc59

数据预处理

对不同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])
subjidstatusleu3pleu2precept_analabx_use
DUMMY1nc40.535.95yes
DUMMY2nc53.927.76yes
DUMMY3nc30.942.53no
DUMMY4nc50.319.37no
DUMMY5nc42.838.84no
DUMMY6nc36.946.74yes

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.value
  ind_p <- clme_obj$p.value.ind
  est_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_gap
    end_p_pos <- max(est_mean) + p_gap
  } else if (est_mean[2] > est_mean[length(est_mean)]) {
    start_p_pos <- max(est_mean) + p_gap
    end_p_pos <- min(est_mean) + p_gap
  } else {
    if (trend == "increase") {
      start_p_pos <- est_mean[2] + p_gap
      end_p_pos <- max(est_mean) + 2 * p_gap
    } else {
      start_p_pos <- max(est_mean) + 2 * p_gap
      end_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)
}
  1. P-value is obtained by linear regression adjusting for antibiotics usage.

  2. P-values were not adjusted for multiple comparisons.

限制参数

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

Primary grouping: Acetate

短链脂肪酸AcetatePrimary 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

短链脂肪酸ButyratePrimary 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

短链脂肪酸PropionatePrimary 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

短链脂肪酸ValeratePrimary 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

短链脂肪酸AcetateSecondary 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

短链脂肪酸ButyrateSecondary 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

短链脂肪酸PropionateSecondary 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

短链脂肪酸ValerateSecondary 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

在这里插入图片描述

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/749737.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

@城规人快来抄作业!转GIS开发月薪12000+

从性价比极低的时薪&#xff0c;到相对稳定的月薪过万&#xff0c;我做对了哪些事情&#xff1f; 今天分享的是城乡规划专业的L拿到GIS开发高薪offer的故事。 初识新中地 该同学是城乡规划专业本科&#xff0c;下面称他为L同学。 L同学是今年夏天在网络上了解了GIS开发和新…

Kafka入门到精通(四)-SpringBoot+Kafka

一丶IDEA创建一个空项目 二丶添加相关依赖 <dependencies><dependency><groupId>org.springframework.boot</groupId><artifactId>spring-boot-starter-web</artifactId></dependency><dependency><groupId>org.springf…

MySQL改密

这里写目录标题 更改登录密码&#xff1a;有权限账号能登录mysql中&#xff1a;有权限账号不能登录mysql中&#xff1a;mysql5.6版本命令mysql5.7版本命令修改密码8.0版本改完后&#xff1a; mysql登录不上了本机安装了5.6后&#xff0c;又安装了mysql8.0 更改登录密码&#xf…

双麒麟系统!RK3588+银河麒麟/开放麒麟,全国产让您的产品更具竞争力

01 银河麒麟嵌入式系统介绍 银河麒麟嵌入式操作系统V10 SP1是为物联网及工业互联网场景设计的安全实时系统&#xff0c;基于Linux内核&#xff0c;采用“分域虚拟化 多域隔离”架构&#xff0c;结合了Linux的丰富生态和RTOS的硬实时能力。 该系统支持主流嵌入式芯片&#x…

“数字政协”平台如何提高政协工作效率?正宇软件助力建设!

随着信息技术的飞速发展&#xff0c;数字化已成为推动各行各业转型升级的重要力量。在政协工作中&#xff0c;数字政协平台的建设与运用&#xff0c;正成为提高政协工作效率、促进民主协商的重要手段。本文将从数字政协平台的功能特点、优势分析以及实践应用等方面&#xff0c;…

【Android】【Compose】Compose里面的Row和Column的简单使用

内容 Row和Column的简单使用方式和常用属性含义 Row 在 Jetpack Compose 中&#xff0c;Row 是一种用于在水平方向排列子元素的布局组件。它类似于传统 Android 中的 LinearLayout&#xff0c;但更加灵活和强大。 Row的代码 Composable inline fun Row(modifier: Modifier…

小九首度回应与小水分手传闻揭秘

#小九首度回应&#xff01;与小水分手传闻揭秘#近日&#xff0c;泰国娱乐圈掀起了一股热议的狂潮&#xff01;传闻中的“金童玉女”组合——“小水”平采娜与“小九”NINE疑似分手的消息&#xff0c;如同巨石投入平静的湖面&#xff0c;激起了千层浪花。而在这股狂潮中&#xf…

高效同步的PWM升压DC/DC转换器 SD6201/SD6201-AF

SD6201是高效同步的PWM升压DC/DC转换器优化为介质提供高效的解决方案电力系统。这些设备在输入电压介于0.9V和4.4V之间&#xff0c;带有1.4MHz固定频率切换。这些功能通过允许使用小型、薄型电感器以及陶瓷电容器。自动PWM/PFM轻负载下的模式切换可节省电力提高了效率。电压在2…

武汉星起航:挂牌上海股权托管交易中心,亚马逊影响力再掀波澜

在全球化日益加深的今天&#xff0c;跨境电商行业正迎来前所未有的发展机遇。而在这个风起云涌的时代&#xff0c;武汉星起航电子商务有限公司以其卓越的实力和前瞻性的战略眼光&#xff0c;成功在上海股权托管交易中心挂牌展示&#xff0c;正式登陆资本市场&#xff0c;这一重…

CSS的媒体查询:响应式布局的利器

关于CSS的媒体查询 CSS媒体查询是CSS层叠样式表(Cascading Style Sheets)中的一个核心功能&#xff0c;它使得开发者能够根据不同的设备特性和环境条件来应用不同的样式规则。这是实现响应式网页设计的关键技术&#xff0c;确保网站或应用能够在多种设备上&#xff0c;包括桌面…

提升用户转化率秘诀!Xinstall的H5拉起应用技术让您领先一步!

在移动互联网时代&#xff0c;App的推广和运营面临着诸多挑战。其中&#xff0c;H5页面如何高效、便捷地拉起应用&#xff0c;成为了一个亟待解决的问题。今天&#xff0c;我们就来谈谈如何利用Xinstall品牌&#xff0c;轻松解决这一痛点&#xff0c;提升用户体验&#xff0c;助…

CentOS 7.9 CDH6.3.2集群生产环境实战部署指南

一、环境准备 1、系统环境&#xff1a; # cat /etc/os-release 2、准备工作&#xff1a; 部署资源分配 节点centos 7.9&#xff08;生产&#xff09;节点规划Postgresql部署组件备注pgsql32c、128G、2TB国产数据库Postgresql&#xff08;翰高&#xff09;可根据实际情况调整…

启动台出现agent app的解决方法~

启动台出现agent app的解决方法&#xff5e; 如果用了战网&#xff0c;Battle.net&#xff0c;在卸载后有一个agent app&#xff0c;启动台删除不掉&#xff0c;应用程序里面没有&#xff0c;怎么办呢&#xff1f; 解决方法&#xff1a;找到这个app所在位置&#xff0c;可以通…

Facebook之梦:数字社交的无限可能

在当今数字化和全球化的时代&#xff0c;社交网络已经成为人们日常生活不可或缺的一部分。作为全球最大的社交平台之一&#xff0c;Facebook不仅连接了数十亿用户&#xff0c;还深刻影响了我们的社交方式、文化交流和信息传播。然而&#xff0c;Facebook所代表的不仅仅是一个网…

深入理解 Dubbo:分布式服务框架的核心原理与实践

目录 Dubbo 概述Dubbo 的架构Dubbo 的关键组件 服务提供者&#xff08;Provider&#xff09;服务消费者&#xff08;Consumer&#xff09;注册中心&#xff08;Registry&#xff09;监控中心&#xff08;Monitor&#xff09;调用链追踪&#xff08;Trace&#xff09; Dubbo 的…

【Java】字节数组 pcm 与 wav 格式互转 (附原理概述)

前言 最近实现了一个文字转语音的功能&#xff0c;语音引擎返回的是pcm格式的数据。需要转化成wav格式前端才能播放。本文首先会给出解决方案&#xff0c;后续会讲背后的原理。 场景 git 仓库 https://github.com/ChenghanY/pcm-wav-converter 1. pcm wav 转化工具类 入参和…

MES管理系统的实施难点以及解决方案

随着智能制造的浪潮席卷全球&#xff0c;MES管理系统成为了众多制造企业提升竞争力的关键武器。MES管理系统以其强大的功能&#xff0c;能够有效连接企业的上层ERP系统与底层自动化设备&#xff0c;实现生产过程的实时监控与优化。然而&#xff0c;实施MES管理系统并非一帆风顺…

Linux通用系统高危漏洞(CVE-2024-1086)修复案例

一、漏洞描述 2024年3月28日&#xff0c; Linux kernel权限提升漏洞&#xff08;CVE-2024-1086&#xff09;的PoC/EXP在互联网上公开&#xff0c;该漏洞的CVSS评分为7.8&#xff0c;目前漏洞细节已经公开披露&#xff0c;美国网络安全与基础设施安全局&#xff08;CISA&#x…

springboot框架使用Netty依赖中解码器的作用及实现详解

在项目开发 有需求 需要跟硬件通信 也没有mqtt 作为桥接 也不能http 请求 api 所以也不能 json字符串这么爽传输 所以要用tcp 请求 进行数据交互 数据还是16进制的 写法 有帧头 什么的 对于这种物联网的这种对接 我的理解就是 我们做的工作就像翻译 把这些看不懂的 字节流 变成…

深圳技术大学oj C : 生成r子集

Description 输出给定序列按字典序的 &#xfffd; 组合&#xff0c;按照所有 &#xfffd; 个元素出现与否的 01 标记串 &#xfffd;&#xfffd;&#xfffd;&#xfffd;−1,...,&#xfffd;1 的字典序输出. 此处01串的字典序指&#xff1a;先输入的数字对应低位&#x…