数据分析:微生物数据的荟萃分析框架

介绍

Meta-analysis of fecal metagenomes reveals global microbial signatures that are specific for colorectal cancer提供了一种荟萃分析的框架,它主要基于常用的Wilcoxon rank-sum test和Blocked Wilcoxon rank-sum test 方法计算显著性,再使用分位数计算分组间的倍数变化,最后通过AUC判断物种的区分分组的能力。最后通过热图和森林图展示筛选到的在不同研究和荟萃分析均有差异的物种。

该框架可用于同类型的微生物荟萃分析。

加载R包

#| warning: false
#| message: falselibrary(tidyverse)
library(readr)
library(coin)
library(pROC)
library(RColorBrewer)
library(cowplot)# rm(list = ls())
options(stringsAsFactors = F)
options(future.globals.maxSize = 1000 * 1024^2)

导入数据

数据下载百度云盘链接: https://pan.baidu.com/s/1VS6S8p5s20vwZ6FyILYoaQ

提取码: g4y3

  • 物种表达谱数据

  • 样本分组信息

#| warning: false
#| message: falsefeat.all <- read.table("./data/meta-CRC-2019/feat_rel_crc.tsv", sep='\t', header=TRUE, stringsAsFactors = FALSE, check.names = FALSE, quote='') %>%as.matrix()meta <- read_tsv('./data/meta-CRC-2019/meta_crc.tsv', show_col_types = FALSE)
  • 其他参数(过滤和检验结果)
#| warning: false
#| message: falsealpha.meta <- 1e-05
alpha.single.study <- 0.005
mult.corr <- 'fdr'
pr.cutoff <- 0.05
log.n0 <- 1e-05
log.n0.func <- 1e-08
study.cols <- c('#2FBFBF', '#177254', '#F2CC30', '#74B347', '#8265CC')

数据预处理

  • 提出研究名称studies

  • 设置block分组,用于后续检验

#| warning: false
#| message: falsestudies <- meta %>% dplyr::pull(Study) %>% unique# block for colonoscopy and study as well
meta <- meta %>%dplyr::filter(!is.na(Sampling_rel_to_colonoscopy)) %>%dplyr::mutate(block = ifelse(Study != 'CN-CRC', Study, paste0(Study, '_', Sampling_rel_to_colonoscopy)))feat.all <- feat.all[, meta$Sample_ID]

荟萃分析

荟萃分析采用了Wilcoxon rank-sum test和Blocked Wilcoxon rank-sum test 两种方法对单个研究和合并所有研究做显著性检验。本次需要计Foldchange(FC)单个研究的pvalue + 所有研究的pvalue(p.val)单个研究和所有研究的AUC(aucs),以下是该代码的计算过程:

  • 先使用Wilcoxon rank-sum test计算每个研究的每个物种在case/control之间的显著性检验结果;

  • 再通过roc函数计算每个研究的每个物种在case/control之间的判别效果;

  • 接着通过分位数quantile计算每个研究的每个物种在case/control之间的倍数变化;

  • 然后通过Blocked Wilcoxon rank-sum test计算所有研究的荟萃差异检验结果;

  • 最后计算所有研究的平均倍数变化作为整体倍数变化和通过roc函数计算每个物种在case/control之间的判别效果。

#| warning: false
#| message: falsep.val <- matrix(NA, nrow = nrow(feat.all), ncol = length(studies)+1, dimnames = list(row.names(feat.all), c(studies, 'all')))
fc <- p.val
aucs.mat <- p.val
aucs.all <- vector('list', nrow(feat.all))cat("Calculating effect size for every feature...\n")
pb <- txtProgressBar(max = nrow(feat.all), style = 3)# caluclate wilcoxon test and effect size for each feature and study
for (f in row.names(feat.all)) {# for each studyfor (s in studies) {x <- feat.all[f, meta %>% dplyr::filter(Study == s) %>% dplyr::filter(Group=='CRC') %>% dplyr::pull(Sample_ID)]y <- feat.all[f, meta %>% dplyr::filter(Study==s) %>% dplyr::filter(Group=='CTR') %>% dplyr::pull(Sample_ID)]# Wilcoxon: 对单个研究的单个物种检验p.val[f, s] <- wilcox.test(x, y, exact=FALSE)$p.value# AUC:评估每个物种区分分组的能力aucs.all[[f]][[s]] <- c(roc(controls=y, cases=x, direction='<', ci=TRUE, auc=TRUE)$ci)aucs.mat[f, s] <- c(roc(controls=y, cases=x, direction='<', ci=TRUE, auc=TRUE)$ci)[2]# FC:使用10分位数计算每个物种的相对丰度再计算Foldchange结果q.p <- quantile(log10(x+log.n0), probs = seq(.1, .9, .05))q.n <- quantile(log10(y+log.n0), probs = seq(.1, .9, .05))fc[f, s] <- sum(q.p - q.n)/length(q.p)}# calculate effect size for all studies combined# Wilcoxon + blocking factor:计算所有研究混合在一起的检验结果d <- data.frame(y = feat.all[f,], x = meta$Group, block = meta$block) %>%dplyr::mutate(x = factor(x),block = factor(block))p.val[f, 'all'] <- coin::pvalue(wilcox_test(y ~ x | block, data = d))# other metricsx <- feat.all[f, meta %>% dplyr::filter(Group=='CRC') %>% dplyr::pull(Sample_ID)]y <- feat.all[f, meta %>% dplyr::filter(Group=='CTR') %>% dplyr::pull(Sample_ID)]# FC: 取所有样本的平均FC结果fc[f, 'all'] <- mean(fc[f, studies])# AUC:合并数据集每个物种区分不同分组样本的能力aucs.mat[f, 'all'] <- c(roc(controls=y, cases=x, direction='<', ci=TRUE, auc=TRUE)$ci)[2]# progressbarsetTxtProgressBar(pb, (pb$getVal()+1))
}
cat('\n')# multiple hypothesis correction
p.adj <- data.frame(apply(p.val, MARGIN=2, FUN=p.adjust, method=mult.corr),check.names = FALSE)

查看结果

查看上述荟萃分析的结果

#| warning: false
#| message: falsehead(p.adj)head(aucs.mat)head(fc)

画图

文章给出的图分成两部分,上部分是热图形式,下半部是森林图。

  • 热图: 展示不同研究显著差异的物种
#| warning: false
#| message: falsespecies.heatmap <- rownames(p.adj)[which(p.adj$all < alpha.single.study)]fc.sign <- sign(fc)
fc.sign[fc.sign == 0] <- 1p.val.signed <- -log10(p.adj[species.heatmap,"all", drop=FALSE]) * fc.sign[species.heatmap, 'all']top.markers <- rownames(p.val.signed[is.infinite(p.val.signed$all) , , drop=FALSE])
p.val.signed[top.markers, 'all'] <- 100 + aucs.mat[top.markers, 'all']species.heatmap.orderd <- rownames(p.val.signed[order(p.val.signed$all), , drop=FALSE])# take only those
fc.mat.plot <- fc[species.heatmap.orderd, ] %>% as.data.frame()
p.vals.plot <- p.adj[species.heatmap.orderd, ]# ##############################################################################
# prepare plotting# colorscheme for fc heatmap 
mx <- max(abs(range(fc.mat.plot, na.rm=TRUE)))
mx <- ifelse(round(mx, digits = 1) < mx, round(mx, digits = 1) + 0.1, round(mx, digits = 1))
brs <- seq(-mx, mx, by=0.05)
num.col.steps <- length(brs) - 1
n <- floor(0.45*num.col.steps)
col.hm <- c(rev(colorRampPalette(brewer.pal(9, 'Blues'))(n)),rep('#FFFFFF', num.col.steps-2*n),colorRampPalette(brewer.pal(9, 'Reds'))(n))
# color scheme for pval heatmap
alpha.breaks <- c(1e-06, 1e-05, 1e-04, 1e-03, 1e-02, 1e-01)
p.vals.bin <- data.frame(apply(p.vals.plot, 2, FUN=.bincode, breaks = c(0, alpha.breaks, 1), include.lowest = TRUE),check.names = FALSE)
p.val.greys <- c(paste0('grey', round(seq(from=10, to=80, length.out = length(alpha.breaks)))), 'white')
names(p.val.greys) <- as.character(1:7)# function to plot both into a grid
plot.single.study.heatmap <- function(x) {# x = "FR-CRC"df.plot <- tibble(species = factor(rownames(p.vals.plot), levels = rev(rownames(p.vals.plot))),p.vals = as.factor(p.vals.bin[[x]]),fc = fc.mat.plot[[x]])g1 <- df.plot %>% ggplot(aes(x = species, y = 1, fill = fc)) + geom_tile() + theme_minimal() + theme(axis.text = element_blank(),axis.ticks = element_blank(),axis.title = element_blank(), panel.grid = element_blank(),panel.background = element_rect(fill=NULL, colour='black'),plot.margin = unit(c(0, 0, 0, 0), 'cm')) + scale_y_continuous(expand = c(0, 0)) + scale_fill_gradientn(colours=col.hm, limits=c(-mx, mx), guide=FALSE)g2 <- df.plot %>% ggplot(aes(x=species, y=1, fill=p.vals)) +geom_tile() + theme_minimal() + theme(axis.text = element_blank(),axis.ticks = element_blank(),axis.title = element_blank(), panel.grid = element_blank(),panel.background = element_rect(fill=NULL, colour='black'),plot.margin = unit(c(0, 0, 0, 0), 'cm')) + scale_y_continuous(expand = c(0, 0)) + scale_fill_manual(values=p.val.greys, na.value='white', guide=FALSE)g.return <- plot_grid(g2, g1, ncol = 1, rel_heights = c(0.25, 0.75))return(g.return)
}# ##############################################################################
# plot# p.value histogram
g1 <- tibble(species = factor(rownames(p.vals.plot), levels = rev(rownames(p.vals.plot))),p.vals = -log10(p.vals.plot$all),colour = p.vals > 5) %>% ggplot(aes(x = species, y = p.vals, fill = colour)) + geom_bar(stat = 'identity') + theme_classic() + xlab('Gut microbial species') + ylab('-log10(q-value)') + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.ticks.x = element_blank(),axis.text.x = element_blank(),panel.background = element_rect(fill = NULL, color = 'black')) + scale_y_continuous(limits = c(0, 15), expand = c(0, 0)) + scale_x_discrete(position = 'top') + scale_fill_manual(values = c('lightgrey', 'darkgrey'), guide = FALSE)g.lst <- lapply(studies, plot.single.study.heatmap)pl1 <- plot_grid(g1, g.lst[[1]], g.lst[[2]], g.lst[[3]], g.lst[[4]],g.lst[[5]], ncol = 1, align = 'v', rel_heights = c(0.3, rep(0.12, 5)))pl1

结果:差异物种在不同研究和整合数据的倍数变化和显著性结果

  • 最上图是物种在整合数据的显著性结果(adjustedPvalue);

  • 接下来的热图是物种在单个研究的显著性结果(上半部)和倍数变化(下半部:红色是富集在CRC,蓝色是CTRL);

  • 森林图: 不同物种在每个研究区分case/control的能力 (通过alpha.meta更严格筛选)

#| warning: false
#| message: false# select and order
marker.set <- rownames(p.val.signed)[abs(p.val.signed$all) > -log10(alpha.meta)]
p.val.signed.red <- p.val.signed[marker.set, ,drop=FALSE]
marker.set.orderd <- rev(rownames(p.val.signed.red[order(p.val.signed.red$all),,drop=FALSE]))# extract those from the auc list
df.plot <- tibble()
for (i in marker.set.orderd){for (s in studies){temp <- aucs.all[[i]][[s]]df.plot <- bind_rows(df.plot, tibble(species=i, study=s,low=temp[1], auc=temp[2], high=temp[3]))}
}df.plot <- df.plot %>% dplyr::mutate(species = factor(species, levels = marker.set.orderd)) %>% dplyr::mutate(study = factor(study, levels = studies))# plot everything
pl2 <- df.plot %>% ggplot(aes(x = study, y = auc)) + geom_linerange(aes(ymin = low, ymax = high), color = 'lightgrey') + geom_point(pch = 23, aes(fill = study)) + facet_grid(~species, scales = 'free_x', space = 'free') + theme_minimal() + scale_y_continuous(limits=c(0, 1)) + theme(panel.grid.major.x = element_blank(),axis.ticks.x = element_blank(),axis.text.x = element_blank(),strip.text = element_text(angle=90, hjust=0)) + scale_fill_manual(values = study.cols, guide = FALSE) + ylab('AUROC') + xlab('Gut microbial species')pl2

  • 合并图: 最后文章呈现的图是经过修改的
#| warning: false
#| message: falsecowplot::plot_grid(pl1, pl2, ncol = 1)

总结

在进行荟萃分析时,本研究采用了一种特定的统计方法——Blocked Wilcoxon rank-sum test,以评估和整合不同研究中的case/control物种的显著性结果。该方法特别适用于处理微生物数据这类稀疏性数据集,因为它能够在计算两组之间的倍数变化时有效避免零值过多的问题。通过使用分位数方法,研究者能够更准确地估计和比较不同组之间的差异,从而提高了分析结果的可靠性和有效性。

对于类似类型的研究,研究者可以采用与本研究相似的分析框架进行荟萃分析。这包括以下几个关键步骤:

  • 数据的收集与整理:确保收集到的数据是高质量的,并且适合进行荟萃分析。
  • 选择合适的统计方法:根据数据的特点选择合适的统计检验方法,如Blocked Wilcoxon rank-sum test,以确保分析的准确性。
  • 数据处理:对于稀疏数据,采用分位数方法来处理零值过多的问题,以提高分析的稳健性。
  • 结果的整合与解释:将不同研究的结果进行整合,并采用适当的统计方法来评估整体的显著性。

通过遵循这样的框架,研究者可以对类似主题的研究进行系统性地分析和比较,从而为该领域的研究提供更深入的见解。

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

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

相关文章

Kafka核心知识点整理,值得收藏!

消息队列应用场景 提高系统性能&#xff1a;通过异步处理减少响应时间。削峰/限流&#xff1a;应对高并发场景。降低系统耦合性&#xff1a;解耦生产者和消费者。 消息队列对比 Kafka&#xff1a;高吞吐量&#xff0c;适合日志收集和传输&#xff0c;适合大型公司。RocketMQ…

SpringBoot启动命令过长

Error running DromaraApplication: Command line is too long. Shorten command line for DromaraApplication or also for Spring Boot default configuration?

线上环境服务器CPU飙升排查

前因 收到线上服务器CPU使用率100%的告警信息。 环境 jdk1.8CentOS Linux &#xff1b;CentOS Linux 排查 查看服务器CPU使用率 果然cpu已经达到了100%了 命令 top 使用arthas工具 使用方式 arthas 执行命令java -jar arthas-boot.jar 然后执行命令 thread 看到有两个…

【linux】Shell脚本三剑客之awk命令的详细用法攻略

✨✨ 欢迎大家来到景天科技苑✨✨ 🎈🎈 养成好习惯,先赞后看哦~🎈🎈 🏆 作者简介:景天科技苑 🏆《头衔》:大厂架构师,华为云开发者社区专家博主,阿里云开发者社区专家博主,CSDN全栈领域优质创作者,掘金优秀博主,51CTO博客专家等。 🏆《博客》:Python全…

【基础教程】Tutorial on Pytorch 结合官方基础文档和个人经验

参考与前言 此教程首次书写于2021年12月份 至 2022年4月份间不断补充&#xff1b;阅读本文时可以对着代码运行查看 官方网址&#xff1a;https://pytorch.org/tutorials/ 【基本从这里翻译而来 更简洁版碎碎念】https://pytorch.org/tutorials/beginner/blitz/cifar10_tutori…

vue3+element-plus 实现动态菜单和动态路由的渲染

在 Vue.js 中&#xff0c;使用 Vue Router 管理路由数据&#xff0c;并将其用于渲染 el-menu&#xff08;Element UI 的菜单组件&#xff09;通常涉及以下几个步骤&#xff1a; 定义路由元数据&#xff1a; 在你的路由配置中&#xff0c;为每个路由项添加 meta 字段&#xff0c…

无人机10公里WiFi图传摄像模组,飞睿智能超清远距离无线监控,智能安防新潮流

在这个科技日新月异的时代&#xff0c;我们对影像的捕捉和传播有了更高的要求。从传统的有线传输到无线WiFi图传&#xff0c;每一次技术的飞跃都为我们带来了全新的视觉体验。今天&#xff0c;我们要探讨的&#xff0c;正是一款具有划时代意义的科技产品——飞睿智能10公里WiFi…

自学网络安全,从小白到大神的破茧之路!

在当今数字化高速发展的时代&#xff0c;网络安全已经成为了至关重要的领域。无论是个人的隐私保护&#xff0c;还是企业、国家的关键信息资产维护&#xff0c;都离不开网络安全的有力保障。出于对这一领域的浓厚兴趣以及对未来职业发展的清晰规划&#xff0c;我毅然决然地踏上…

0722_驱动1 字符设备驱动框架

一、字符设备驱动框架 字符设备驱动按照字节流进行访问&#xff0c;并且只能顺序访问 设备号一共有32位&#xff0c;主设备号&#xff08;高12位&#xff09;次设备号&#xff08;低20位&#xff09; 二、注册/注销字符设备驱动API接口 2.1、注册字符设备驱动(入口) #include &…

树莓派智能家居中枢

一个先进的枢纽&#xff0c;使智能家居系统更智能、更可定制、更易于控制 Homey Pro由树莓派 Compute Module 4 供电,Homey Pro 为用户提供了一个单一界面,用于控制和监控来自不同品牌的所有智能家居设备。它完全在本地网络上运行,而不是依赖云端,从而实现了最低的延迟、最高的…

深入理解CSS中的变量(应用篇)

在现代Web开发中,换肤功能已经成为提升用户体验的重要手段之一。通过使用CSS变量和JavaScript,我们可以轻松实现动态换肤功能。本文将介绍如何动态生成和应用CSS变量来实现换肤效果。 1. 定义基础CSS变量 首先,我们在CSS中定义一些基础的CSS变量。这些变量将用于存储不同主…

测试面试宝典(三十三)—— 接口测试有没有测试出什么问题?

在之前的接口测试工作中&#xff0c;确实发现了一些问题。比如&#xff0c;在对某关键业务接口进行测试时&#xff0c;发现当输入的参数值超出正常范围时&#xff0c;接口没有按照预期返回错误提示&#xff0c;而是出现了系统崩溃的情况。 还有一次&#xff0c;在测试一个数据…

关于任务栏设置闪退的一种解决方法

昨天我用了下360安全卫士&#xff0c;就发现我的任务栏变成了小任务栏模式。我想把它关掉&#xff0c;于是想要打开任务栏设置&#xff0c;然后就发现任务栏设置闪退。 我去晚上找了许多方法&#xff0c;像是卸载360&#xff0c;用指令修补破损文件&#xff0c;重启电脑等等&am…

【PostgreSQL教程】PostgreSQL 删除表格

博主介绍:✌全网粉丝20W+,CSDN博客专家、Java领域优质创作者,掘金/华为云/阿里云/InfoQ等平台优质作者、专注于Java技术领域✌ 技术范围:SpringBoot、SpringCloud、Vue、SSM、HTML、Nodejs、Python、MySQL、PostgreSQL、大数据、物联网、机器学习等设计与开发。 感兴趣的可…

docker发布镜像到自己远程私有仓库

1、登录docker hub创建自己的仓库地址&#xff1a;https://hub.docker.com/repository/create 输入仓库名称 2.构建镜像 略过。。。。请自己查找别的资料&#xff0c;此篇文章只讲述镜像推送到远程 3.推送 假设你已经构建了一个镜像 web/online-editor:latest&#xff0c;现…

推荐系统三十六式学习笔记:工程篇.常见架构25|Netflix个性化推荐架构

目录 架构的重要性经典架构1.数据流2.在线层3.离线层4.近线层 简化架构总结 你是否曾经觉得算法就是推荐系统的全部&#xff0c;即便不是全部&#xff0c;至少也是嫡长子&#xff0c;然而实际上&#xff0c;工程实现才是推荐系统的骨架。如果没有好的软件实现&#xff0c;算法不…

构建安全的单点登录(SSO)系统:PHP框架实现指南

引言 单点登录&#xff08;SSO&#xff09;是一种允许用户使用单一身份验证来访问多个相关但独立的系统或应用程序的技术。在企业环境中&#xff0c;SSO可以极大地提高效率&#xff0c;减少密码管理的复杂性&#xff0c;并增强安全性。本文将详细介绍如何在PHP框架中实现SSO&a…

Python面试题:使用Matplotlib和Seaborn进行数据可视化

使用Matplotlib和Seaborn进行数据可视化是数据分析中非常重要的一部分。以下示例展示了如何使用这两个库来创建各种图表&#xff0c;包括基本的线图、柱状图、散点图和高级的分类数据可视化图表。 安装 Matplotlib 和 Seaborn 如果你还没有安装这两个库&#xff0c;可以使用以…

树 ----- 基础学习

树 树&#xff1a;n&#xff08;n>0&#xff09;个结点的有限集合。n 0 ,空树。 在任意一个非空树中&#xff0c; 1&#xff0c;有且仅有一个特定的根结点 2&#xff0c;当n>1 时&#xff0c;其余结点可分为m个互不相交的有限集合T1,T2,T3.。。。。Tm&#xff0c;其中每…