文章MSM_metagenomics(三):Alpha多样性分析

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

  • WX公zhong号:生信学习者
  • Xiao hong书:生信学习者
  • 知hu:生信学习者
  • CDSN:生信学习者2

介绍

本教程使用基于R的函数来估计微生物群落的香农指数和丰富度,使用MetaPhlAn profile数据。估计结果进一步进行了可视化,并与元数据关联,以测试统计显著性。

数据

大家通过以下链接下载数据:

  • 百度网盘链接:https://pan.baidu.com/s/1f1SyyvRfpNVO3sLYEblz1A
  • 提取码: 请关注WX公zhong号_生信学习者_后台发送 复现msm 获取提取码

R 包

  • SummarizedExperiment
  • mia
  • ggpubr
  • ggplot2
  • lfe

Alpha diversity estimation and visualization

使用alpha_diversity_funcs.R计算alpha多样性和可视化。

  • 代码
SE_converter <- function(md_rows, tax_starting_row, mpa_md) {# SE_converter function is to convery metadata-wedged mpa table into SummarisedExperiment structure.# md_rows: a vector specifying the range of rows indicating metadata.# tax_starting_row: an interger corresponding to the row where taxonomic abundances start.# mpa_md: a metaphlan table wedged with metadata, in the form of dataframe.md_df <- mpa_md[md_rows,] # extract metadata part from mpa_md tabletax_df <- mpa_md[tax_starting_row: nrow(mpa_md),] # extract taxonomic abundances part from mpa_md table### convert md_df to a form compatible with SummarisedExperiment ### SE_md_df <- md_df[, -1]rownames(SE_md_df) <- md_df[, 1]SE_md_df <- t(SE_md_df)### convert md_df to a form compatible with SummarisedExperiment ###### prep relab values in a form compatible with SummarisedExperiment ###SE_relab_df <- tax_df[, -1]rownames(SE_relab_df) <- tax_df[, 1]col_names <- colnames(SE_relab_df)SE_relab_df[, col_names] <- apply(SE_relab_df[, col_names], 2, function(x) as.numeric(as.character(x)))### prep relab values in a form compatible with SummarisedExperiment ###SE_tax_df <- tax_df[, 1:2]rownames(SE_tax_df) <- tax_df[, 1]SE_tax_df <- SE_tax_df[-2]colnames(SE_tax_df) <- c("species")SE_data <- SummarizedExperiment::SummarizedExperiment(assays = list(relative_abundance = SE_relab_df),colData = SE_md_df,rowData = SE_tax_df)SE_data
}est_alpha_diversity <- function(se_data) {# This function is to estimate alpha diversity (shannon index and richness) of a microbiome and output results in dataframe.# se_data: the SummarizedExperiment data structure containing metadata and abundance values.se_data <- se_data |>mia::estimateRichness(abund_values = "relative_abundance", index = "observed")se_data <- se_data |>mia::estimateDiversity(abund_values = "relative_abundance", index = "shannon")se_alpha_div <- data.frame(SummarizedExperiment::colData(se_data))se_alpha_div
}make_boxplot <- function(df, xlabel, ylabel, font_size = 11, jitter_width = 0.2, dot_size = 1, font_style = "Arial", stats = TRUE, pal = NULL) {# This function is to create a boxplot using categorical data.# df: The dataframe containing microbiome alpha diversities, e.g. `shannon` and `observed` with categorical metadata.# xlabel: the column name one will put along x-axis.# ylabel: the index estimate one will put along y-axis.# font_size: the font size, default: [11]# jitter_width: the jitter width, default: [0.2]# dot_size: the dot size inside the boxplot, default: [1]# font_style: the font style, default: `Arial`# pal: a list of color codes for pallete, e.g. c(#888888, #eb2525). The order corresponds the column order of boxplot.# stats: wilcox rank-sum test. default: TRUEif (stats) {nr_group = length(unique(df[, xlabel])) # get the number of groupsif (nr_group == 2) {group_pair = list(unique(df[, xlabel]))ggpubr::ggboxplot(data = df, x = xlabel, y = ylabel, color = xlabel,palette = pal, ylab = ylabel, xlab = xlabel,add = "jitter", add.params = list(size = dot_size, jitter = jitter_width)) +ggpubr::stat_compare_means(comparisons = group_pair, exact = T, alternative = "less") +ggplot2::stat_summary(fun.data = function(x) data.frame(y = max(df[, ylabel]), label = paste("Mean=",mean(x))), geom="text") +ggplot2::theme(text = ggplot2::element_text(size = font_size, family = font_style))} else {group_pairs = my_combn(unique((df[, xlabel])))ggpubr::ggboxplot(data = df, x = xlabel, y = ylabel, color = xlabel,palette = pal, ylab = ylabel, xlab = xlabel,add = "jitter", add.params = list(size = dot_size, jitter = jitter_width)) +ggpubr::stat_compare_means() + ggpubr::stat_compare_means(comparisons = group_pairs, exact = T, alternative = "greater") +ggplot2::stat_summary(fun.data = function(x) data.frame(y= max(df[, ylabel]), label = paste("Mean=",mean(x))), geom="text") +ggplot2::theme(text = ggplot2::element_text(size = font_size, family = font_style))}} else {ggpubr::ggboxplot(data = df, x = xlabel, y = ylabel, color = xlabel,palette = pal, ylab = ylabel, xlab = xlabel,add = "jitter", add.params = list(size = dot_size, jitter = jitter_width)) +ggplot2::theme(text = ggplot2::element_text(size = font_size, family = font_style))}
}my_combn <- function(x) {combs <- list()comb_matrix <- combn(x, 2)for (i in 1: ncol(comb_matrix)) {combs[[i]]  <- comb_matrix[,i]}combs
}felm_fixed <- function(data_frame, f_factors, t_factor, measure) {# This function is to perform fixed effect linear modeling# data_frame: a dataframe containing measures and corresponding effects  # f_factors: a vector of header names in the dataframe which represent fixed effects# t_factors: test factor name in the form of string# measure: the measured values in column, e.g., shannon or richness
#   all_factors <- c(t_factor, f_factors)
#   for (i in all_factors) {
#     vars <- unique(data_frame[, i])
#     lookup <- setNames(seq_along(vars) -1, vars)
#     data_frame[, i] <- lookup[data_frame[, i]]
#   }
#   View(data_frame)str1 <- paste0(c(t_factor, paste0(f_factors, collapse = " + ")), collapse = " + ")str2 <- paste0(c(measure, str1), collapse = " ~ ")felm_stats <- lfe::felm(eval(parse(text = str2)), data = data_frame)felm_stats
}

加载一个包含元数据和分类群丰度的合并MetaPhlAn profile文件

mpa_df <- data.frame(read.csv("./data/merged_abundance_table_species_sgb_md.tsv", header = TRUE, sep = "\t"))
sampleP057P054P052P049
sexual_orientationMSMMSMMSMNon-MSM
HIV_statusnegativepositivepositivenegative
ethnicityCaucasianCaucasianCaucasianCaucasian
antibiotics_6monthYesNoNoNo
BMI_kg_m2_WHOObeseClassIOverweightNormalOverweight
Methanomassiliicoccales_archaeon0.00.00.00.01322
Methanobrevibacter_smithii0.00.00.00.19154

接下来,我们将数据框转换为SummarizedExperiment数据结构,以便使用SE_converter函数继续分析,该函数需要指定三个参数:

  • md_rows: a vector specifying the range of rows indicating metadata. Note: 1-based.
  • tax_starting_row: an interger corresponding to the row where taxonomic abundances start.
  • mpa_md: a metaphlan table wedged with metadata, in the form of dataframe.
SE <- SE_converter(md_rows = 1:5,tax_starting_row = 6, mpa_md = mpa_df)SE                   class: SummarizedExperiment
dim: 1676 139
metadata(0):
assays(1): relative_abundance
rownames(1676): Methanomassiliicoccales_archaeon|t__SGB376GGB1567_SGB2154|t__SGB2154 ... Entamoeba_dispar|t__EUK46681Blastocystis_sp_subtype_1|t__EUK944036
rowData names(1): species
colnames(139): P057 P054 ... KHK16 KHK11
colData names(5): sexual_orientation HIV_status ethnicityantibiotics_6month BMI_kg_m2_WHO

接下来,我们可以使用est_alpha_diversity函数来估计每个宏基因组样本的香农指数和丰富度。

alpha_df <- est_alpha_diversity(se_data = SE)
alpha_df
sexual_orientationHIV_statusethnicityantibiotics_6monthBMI_kg_m2_WHOobservedshannon
P057MSMnegativeCaucasianYesObeseClassI1343.1847
P054MSMpositiveCaucasianNoOverweight1412.1197
P052MSMpositiveCaucasianNoNormal1522.5273

为了比较不同组之间的alpha多样性差异,我们可以使用make_boxplot函数,并使用参数:

  • df: The dataframe containing microbiome alpha diversities, e.g. shannon and observed with categorical metadata.
  • xlabel: the column name one will put along x-axis.
  • ylabel: the index estimate one will put along y-axis.
  • font_size: the font size, default: [11]
  • jitter_width: the jitter width, default: [0.2]
  • dot_size: the dot size inside the boxplot, default: [1]
  • font_style: the font style, default: Arial
  • pal: a list of color codes for pallete, e.g. c(#888888, #eb2525). The order corresponds the column order of boxplot.
  • stats: wilcox rank-sum test. default: TRUE
shannon <- make_boxplot(df = alpha_df,xlabel = "sexual_orientation",ylabel = "shannon",stats = TRUE,pal = c("#888888", "#eb2525"))richness <- make_boxplot(df = alpha_df,xlabel = "sexual_orientation", ylabel = "observed",stats = TRUE,pal = c("#888888", "#eb2525"))
multi_plot <- ggpubr::ggarrange(shannon, richness, ncol = 2)
ggplot2::ggsave(file = "shannon_richness.svg", plot = multi_plot, width = 4, height = 5)

请添加图片描述

通过固定效应线性模型估计关联的显著性

在宏基因组分析中,除了感兴趣的变量(例如性取向)之外,通常还需要处理多个变量(例如HIV感染和抗生素使用)。因此,在测试微生物群落矩阵(例如香农指数或丰富度)与感兴趣的变量(例如性取向)之间的关联时,控制这些混杂效应非常重要。在这里,我们使用基于固定效应线性模型的felm_fixed函数,该函数实现在R包lfe 中,以估计微生物群落与感兴趣变量之间的关联显著性,同时控制其他变量的混杂效应。

  • data_frame: The dataframe containing microbiome alpha diversities, e.g. shannon and observed with multiple variables.
  • f_factors: A vector of variables representing fixed effects.
  • t_factor: The variable of interest for testing.
  • measure: The header indicating microbiome measure, e.g. shannon or richness
lfe_stats <- felm_fixed(data_frame = alpha_df,f_factors = c(c("HIV_status", "antibiotics_6month")),t_factor = "sexual_orientation",measure = "shannon")
summary(lfe_stats)
Residuals:Min      1Q  Median      3Q     Max 
-2.3112 -0.4666  0.1412  0.5200  1.4137 Coefficients:Estimate Std. Error t value Pr(>|t|)    
(Intercept)            3.62027    0.70476   5.137 9.64e-07 ***
sexual_orientationMSM  0.29175    0.13733   2.125   0.0355 *  
HIV_statuspositive    -0.28400    0.14658  -1.937   0.0548 .  
antibiotics_6monthNo  -0.10405    0.67931  -0.153   0.8785    
antibiotics_6monthYes  0.01197    0.68483   0.017   0.9861    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.6745 on 134 degrees of freedom
Multiple R-squared(full model): 0.07784   Adjusted R-squared: 0.05032 
Multiple R-squared(proj model): 0.07784   Adjusted R-squared: 0.05032 
F-statistic(full model):2.828 on 4 and 134 DF, p-value: 0.02725 
F-statistic(proj model): 2.828 on 4 and 134 DF, p-value: 0.02725

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

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

相关文章

签约喜报 | Smartbi朋友圈又添新朋友啦~

近期&#xff0c;一系列业界翘楚如国际精密集团、惠达卫浴、华天科技、中国人寿等新老朋友纷纷携手Smartbi&#xff0c;共同探索数据驱动业务的新路径、新思路。 Smartbi数10年专注于商业智能BI与大数据分析软件与服务&#xff0c;为各行各业提供提供一站式商业智能平台&#x…

视频转换器在线哪个好?让视频播放不受格式限制

在日常的视频观看中&#xff0c;我们可能会遇到视频格式与设备不兼容的问题&#xff0c;导致无法顺畅播放。这就像是缺少了播放的钥匙&#xff0c;让人无法享受视频内容。 面对视频格式不兼容的挑战&#xff0c;选择合适的转换工具至关重要。但不用担心&#xff0c;本文将分享…

通义千问2(Qwen2)大语言模型在PAI-QuickStart的微调、评测与部署实践

Qwen2&#xff08;通义千问2&#xff09;是阿里云最近推出的开源大型语言模型系列&#xff0c;相比2月推出的Qwen1.5&#xff0c;Qwen2实现了整体性能的代际飞跃&#xff0c;大幅提升了代码、数学、推理、指令遵循、多语言理解等能力。其中&#xff0c;Qwen2系列包含5个尺寸的预…

好像也没那么失望!SD3玩起来,Stable Diffusion 3工作流商业及广告设计(附安装包)

今天基于SD3 base 工作流来尝试进行下广告设计&#xff0c;这要是一配上设计文案&#xff0c;视觉感就出来了。下面来看看一些效果展示~ SD3 Medium模型及ComfyUI工作流下载地址&#xff1a;文末领取&#xff01; 1.清凉夏日——西瓜音乐会 提示词&#xff1a; a guitar wi…

停止游戏中的循环扣血显示

停止游戏中循环扣血并显示的具体实现方式会依赖于你的代码结构和游戏的逻辑。通常情况下&#xff0c;你可以通过以下方式来实现停止循环扣血和显示&#xff1a; 1、问题背景 在使用 Python 代码为游戏开发一个生命值条时&#xff0c;遇到了一个问题。代码使用了循环来减少生命…

Studio One软件最新版下载及详细安装教程

Studio One 6是一款功能丰富、专业级的音乐制作软件&#xff0c;它具备灵活的工作流程和高效的团队协作能力&#xff0c;能帮助用户实现高质量的音乐创作和制作。 智能模板更快的启动&#xff0c;全新的智能模板为你手头的任务提供了必要的工具集&#xff0c;包括基本录制、混音…

【前端】HTML5基础

目录 0 参考1 网页1.1 什么是网页1.2 什么是HTML1.3 网页的形成 2 浏览器2.1 常用的浏览器2.2 浏览器内核 3 Web标准3.1 为什么需要Web标准3.2 Web标准的构成 4 HTML 标签4.1 HTML语法规范4.1.1 基本语法概述4.1.2 标签关系4.1.2.1 包含关系4.1.2.2 并列关系 4.2 HTML基本结构标…

03.VisionMaster 机器视觉 位置修正 工具

VisionMaster 机器视觉 位置修正 工具 官方解释&#xff1a;位置修正是一个辅助定位、修正目标运动偏移、辅助精准定位的工具。可以根据模板匹配结果中的匹配点和匹配框角度建立位置偏移的基准&#xff0c;然后再根据特征匹配结果中的运行点和基准点的相对位置偏移实现ROI检测…

远程医疗服务包含哪些服务内容?

在当今数字化时代&#xff0c;远程医疗服务正在迅速崛起&#xff0c;成为医疗保健领域的一项重要创新。通过远程医疗服务&#xff0c;患者可以足不出户就能获得医疗服务。那么远程医疗究竟能提供哪些服务呢?下面我们就来看看。 1. 远程咨询 远程咨询是远程医疗服务的基础&…

内网安全【2】-域防火墙

1.判断什么时候用代理 2.判断什么时候用隧道 3.判断出网和不出网协议 4.如何使用代理建立节点并连接 5.如何使用隧道技术封装协议上线 6.判断哪些代理或隧道情况选择放弃 代理技术&#xff1a;解决网络通讯不通的问题(利用跳板机建立节点后续操作)&#xff08;网络设置导…

C# WinForm —— 35 StatusStrip 介绍

1. 简介 状态栏 StatusStrip&#xff0c;默认在软件的最下方&#xff0c;用于显示系统时间、版本、进度条、账号、角色信息、操作位置信息等 可以在状态栏中添加的控件类型有&#xff1a;StatusLabel、ProgressBar、DropDownButton、SplitButton 2. 属性 属性解释(Name)控…

注册中心理论学习

注册中心介绍 注册中心&#xff08;也称为服务注册中心或服务发现服务&#xff09;是微服务架构中的一个关键组件&#xff0c;它负责服务的注册与发现。在微服务体系中&#xff0c;服务实例的数量和位置是动态变化的&#xff0c;注册中心提供了一个集中的地方来存储这些信息&a…

前端问题整理

Vue vue mvvm&#xff08;Model-View-ViewModel&#xff09;架构模式原理 Model 是数据层&#xff0c;即 vue 实例中的数据View 是视图层&#xff0c; 即 domViewModel&#xff0c;即连接Model和Vue的中间层&#xff0c;Vue实例就是ViewModelViewModel 负责将 Model 的变化反映…

测试基础13:测试用例设计方法-错误推断、因果图判定表

课程大纲 1、错误推测法 靠主观经验和直觉来推测可能容易出现问题的功能或场景&#xff0c;设计相关测试用例进行验证。 2、因果图&判定表 2.1定义 因果图和判定表是分析和表达多逻辑条件下&#xff0c;执行不同操作的情况的工具。 &#xff08;因果图和判定表配合使用&a…

ETL可视化工具 DataX -- 简介( 一)

引言 DataX 系列文章&#xff1a; ETL可视化工具 DataX – 安装部署 ( 二) 1.1 DataX 1.1.1 Data X概览 DataX 是阿里云DataWorks数据集成的开源版本&#xff0c;在阿里巴巴集团内被广泛使用的离线数据同步工具/平台。DataX 实现了包括 MySQL、Oracle、OceanBase、SqlServ…

(el-Transfer)操作(不使用 ts):Element-plus 中 Select 组件动态设置 options 值需求的解决过程

Ⅰ、Element-plus 提供的Select选择器组件与想要目标情况的对比&#xff1a; 1、Element-plus 提供Select组件情况&#xff1a; 其一、Element-ui 自提供的Select代码情况为(示例的代码)&#xff1a; // Element-plus 提供的组件代码: <template><div class"f…

上位机图像处理和嵌入式模块部署(h750和市场上的开发板)

【 声明&#xff1a;版权所有&#xff0c;欢迎转载&#xff0c;请勿用于商业用途。 联系信箱&#xff1a;feixiaoxing 163.com】 目前在电商网站上面&#xff0c;关于h750的开发板很多。一种是某原子和某火出品的板子&#xff0c;这一类的板子就是做的比较大&#xff0c;功能比…

【LeetCode刷题】前缀和解决问题:560.和为k的子数组

【LeetCode刷题】Day 16 题目1&#xff1a;560.和为k的子数组思路分析&#xff1a;思路1&#xff1a;前缀和 哈希表 题目1&#xff1a;560.和为k的子数组 思路分析&#xff1a; 问题1&#xff1a;怎样找到数组所有子数组&#xff1f; 方式一&#xff1a;暴力枚举出来&#x…

AI模型部署:一文搞定Triton Inference Server的常用基础配置和功能特性

前言 Triton Inference Server是由NVIDIA提供的一个开源模型推理框架&#xff0c;在前文《AI模型部署&#xff1a;Triton Inference Server模型部署框架简介和快速实践》中对Triton做了简介和快速实践&#xff0c;本文对Triton的常用配置和功能特性做进一步的汇总整理&#xf…

STM32单片机选型方法

一.STM32单片机选型方法 1.首先要确定需求&#xff1a; 性能需求&#xff1a;根据应用的复杂度和性能要求&#xff0c;选择合适的CPU性能和主频。 内存需求&#xff1a;确定所需的内存大小&#xff0c;包括RAM和Flash存储空间。 外设需求&#xff1a;根据应用所需的功能&…