GMSB文章六:微生物SCFA关联分析

news2024/11/16 10:25:30

欢迎大家关注全网生信学习者系列:

  • 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

在这里插入图片描述

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

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

蒸汽架空管道中的关键守护者:滑动管托、导向管托与固定管托

蒸汽架空管道中的关键守护者&#xff1a;滑动管托、导向管托、固定管托与补偿器的重要角色在蒸汽架空管道系统中&#xff0c;每一个组件都扮演着不可或缺的角色&#xff0c;共同确保管道的安全、高效运行。今天&#xff0c;我们就来深入探讨滑动管托、导向管托、固定管托以及补…

信息安全时代,大学生是否有必要考取NISP证书?

在数字化浪潮席卷全球的今天&#xff0c;信息安全已成为国家、企业乃至个人都必须正视的重要议题。作为新时代的大学生&#xff0c;我们身处这个信息爆炸的时代&#xff0c;如何提升自己的信息安全素养&#xff0c;成为了一个值得深思的问题。而NISP(国家信息安全水平考试)证书…

单例模式(下)

文章目录 文章介绍步骤安排及单例讲解step1&#xff1a;注册单例类型&#xff08;main.cpp&#xff09;step2&#xff1a;定义类和私有构造函数&#xff08;keyboardinputmanager.h&#xff09;step3:&#xff08;keyboardinputmanager.cpp&#xff09;step4&#xff1a;在qml中…

云端智慧,赋能风电场:工业级控制系统云监控网关

风力发电场监控平台实现对风电场的运行状态和风机的实时数据进行监测、控制和管理&#xff0c;提高风电场的可靠性和运行效率&#xff0c;降低维护成本&#xff0c;实现智能化管理。 风机机组PLC、多功能仪表、无线测温、温度变送器、档位变送器、设备接入网关上传数据服务器。…

第1章 框架学习的基石与实战策略

第1章框架学习的基石与实战策略 1.1 框架学习的引路人&#xff1a;权威教程的重要性 在编程的世界里&#xff0c;掌握一个框架就像是装备了一套精良的工具&#xff0c;这不仅能显著提升开发速度&#xff0c;还能展现一个程序员的专业水平。对于那些刚刚踏入编程领域的初学者来…

在线开发、实时交互 | 三维天地低代码开发平台助力提高项目交付速度

1.什么是低代码开发平台? 低代码开发平台基于北京三维天地科技股份有限公司自研原生技术架构研发。三维天地作为国内知名的检验检测信息化领域软件开发服务商,拥有多项自主知识产权及自主研发核心技术,致力于为客户提供信息化整体解决方案及相关软件产品与服务。 三维天地低…

【Java面试场景题】如何优化系统架构设计来缓解流量压力提升并发性能?

一、问题解析 我会以直播互动为例&#xff0c;带你看看读多写多的情况下如何应对流量压力。- 一般来说&#xff0c;这种服务多数属于实时互动服务&#xff0c;因为时效性要求很高&#xff0c;导致很多场景下&#xff0c;我们无法用读缓存的方式来降低核心数据的压力。所以&…

B端系统:配置页面如何设计,这可是用户体验的关键的关键。

提升配置页面体验的十大原则 设计B端系统的配置页面时&#xff0c;用户体验确实是非常关键的。以下是一些设计原则和建议&#xff0c;可以帮助提高配置页面的用户体验&#xff1a; 简洁明了&#xff1a;配置页面应该尽量简洁明了&#xff0c;避免过多的复杂选项和信息。使用清…

基于先验知识引导的三域Transformer-GAN,直接从低计数正电子发射断层扫描图像重建| 文献速递-先进深度学习疾病诊断

Title 题目 Prior Knowledge-guided Triple-Domain Transformer-GAN for Direct PET Reconstruction from Low-Count Sinograms 基于先验知识引导的三域Transformer-GAN&#xff0c;用于直接从低计数正电子发射断层扫描图像重建 01 文献速递介绍 正电子发射断层扫描&…

动手学深度学习(Pytorch版)代码实践 -计算机视觉-38实战Kaggle比赛:图像分类 (CIFAR-10)

38实战Kaggle比赛&#xff1a;图像分类 (CIFAR-10) 比赛链接&#xff1a;CIFAR-10 - Object Recognition in Images | Kaggle 导入包 import os import glob import pandas as pd import numpy as np import torch import torchvision from torch.utils.data import Dataset…

【应届应知应会】Linux常用指令

SueWakeup 个人主页&#xff1a;SueWakeup 系列专栏&#xff1a;学习技术栈 个性签名&#xff1a;保留赤子之心也许是种幸运吧 本文封面由 凯楠&#x1f4f8;友情提供 目录 文件与目录管理 目录操作命令&#xff1a; ls [选项] [目录或文件] mkdir 文件操作命令&#xf…

MacOS java多版本安装与管理

安装sdkman curl -s "https://get.sdkman.io" | bashsource "$HOME/.sdkman/bin/sdkman-init.sh"sdk version正常出现sdkman版本号就安装成功了 # 安装java # 安装java8 sdk install java 8.0.412.fx-zulu建议和上述一样安装 fx-zulu 的jdk&#xff0c…

谷歌SEO在外贸推广中的应用效果如何?

谷歌SEO在外贸推广中非常有效。通过优化网站&#xff0c;可以提高在搜索结果中的排名&#xff0c;这意味着更多的潜在客户会看到你的产品和服务。 一个高排名的网站能带来更多自然流量&#xff0c;不需要花费广告费用。这种流量通常质量较高&#xff0c;因为用户是主动搜索相关…

【仿真建模-anylogic】Scale解析

Author&#xff1a;赵志乾 Date&#xff1a;2024-06-27 Declaration&#xff1a;All Right Reserved&#xff01;&#xff01;&#xff01; 1. 应用场景 Scale是比例尺&#xff0c;用于长度单位和像素之间的换算&#xff0c;anylogic默认为每个agent生成一个scale&#xff0c;…

1.iptables

iptables 防火墙iptables工作流程iptables表与链filter表nat表 防火墙 防火墙开源iptables、firewalld管理控制网络流量、封端口、封IP、nat、&#xff08;snat、dnat&#xff09;映射 共享上网硬件防火墙思科、华三等、深信服、路由器内置防火墙保护内部网络、检测和阻挡恶意…

ISO 50001能源管理体系:激活绿色动能和共塑可持续发展

在当今全球化加速和工业化水平不断提高的背景下&#xff0c;能源消费呈现出前所未有的增长趋势。然而&#xff0c;能源资源的有限性、能源价格的波动以及能源消费对环境造成的影响&#xff0c;尤其是温室气体排放导致的全球气候变化问题&#xff0c;已经成为全球关注的焦点。为…

C++之STL(十二)

1、容器适配器 #include <iostream> #include <stack> #include <list> #include <queue> #include <functional> #include <iterator>using namespace std;int main() {// 栈&#xff08;先进后出filo&#xff09;stack<int, list<…

Linux-笔记 嵌入式gdb远程调试

目录 前言 实现 1、内核配置 2、GDB移植 3、准备调试程序 4、开始调试 前言 gdb调试器是基于命令行的GNU项目调试器&#xff0c;通过gdb工具我们可以实现许多调试手段&#xff0c;同时gdb支持多种语言&#xff0c;兼容性很强。 在桌面 Linux 系统&#xff08;如 Ubuntu、Cent…

基于SpringBoot校园一卡通系统设计和实现(源码+LW+调试文档+讲解等)

&#x1f497;博主介绍&#xff1a;✌全网粉丝10W,CSDN作者、博客专家、全栈领域优质创作者&#xff0c;博客之星、平台优质作者、专注于Java、小程序技术领域和毕业项目实战✌&#x1f497; Java精品实战案例《600套》 2023-2025年最值得选择的Java毕业设计选题大全&#xff1…

Zabbix对接Elasticsearch(ES)数据库(未成功)

0.需求分析 不管zabbix的后端数据库是oracle还是mysql&#xff0c;当zabbix监控的量级达到了一定程度后&#xff0c;那么对数据库的性能是一个非常严峻的挑战。特别是对历史数据的查询&#xff0c;将会变得非常非常的慢&#xff0c;别告诉我可以建索引优化&#xff0c;当量级达…