江苏龙网

  1. 主页 > 生活百科 > >

R语言完美重现STAMP结果图

R语言

写在前面
之前写过一篇关于统计学软件STAMP的教程，使用STAMP对微生物群落数据进行统计学分析还是挺方便的，尤其是对R语言不是很熟悉的朋友来说，图形化的界面和相对简单的操作还是挺友好的 。

  • STAMP——微生物组间差异统计分析 简明教程 中文帮助文档
我们通常使用的STAMP的结果主要就是两组数据之间差异性检验的被称作Extended error bar(扩展柱状图)的图像 。
由于STAMP的结果图相对固定，可修改的图像参数有限，经常会遇到一些问题，比如靶标物种或功能基因名字过程就会导致显示不全，在与其它图像拼接成一副图的时候也会出现字号太小导致看不清楚的问题 。
正好前几天在群里有人询问了这个图有没有其它的绘制办法，今天就给大家带来一个使用ggplot2绘制Extended error bar的方法 。
R语言完美重现STAMP结果图

文章插图
 
数据准备这里我将使用一套同一环境位点水体和沉积物16S扩增子测序的PICRUSt功能预测结果作为示例 。
选择的是KEGG L2水平的功能预测的相对丰度数据 。
绘图的数据文件有两个，一个是丰度数据，另一个是样本分组数据 。
R语言完美重现STAMP结果图

文章插图
R语言完美重现STAMP结果图

文章插图
后台回复“STAMP”获取示例文件和完整代码 。
首先将数据导入R环境，我是首先过滤掉了平均丰度低于1%的功能分类 。
data <- read.table("KEGG_L2.txt",header = TRUE,row.names = 1,sep = "t")group <- read.table("group.txt",header = FALSE,sep = "t")library(tidyverse)data <- data*100data <- data %>% filter(Apply(data,1,mean) > 1)??如果使用不同分类学水平的微生物数据或着更深层次的功能注释数据，由于物种或功能基因种类较多，可能会导致结果中具有差异的数目特别多，比如几十个差异物种放在一副图里基本上是不可能看清的，这种时候就要对数据进行过滤，去除低丰度的物种或基因，具体的过滤标准请根据自身数据情况自行确定 。
 
统计学检验在绘制Extended error bar之前首先要对数据进行差异显著性检验，以选出丰度在不同组间具有显著差异的物种或功能基因，这里以两组数据为例，使用的检验方式通常为t-test和Wilcox秩和检验 。
当分析数据符合正态分布时，使用t-test，如不符合正态分布，则使用Wilcox秩和检验 。
 
t-test首先对数据进行调整，构建用于t-test的数据框 。
data <- t(data)data1 <- data.frame(data,group$V2)colnames(data1) <- c(colnames(data),"Group")data1$Group <- as.factor(data1$Group)由于R语言的t-test一次只能分析一列数据，在网上搜到了一个批量进行t-test的方法，感觉是最简便的了 。
【R语言完美重现STAMP结果图】首先使用select_if选择格式为数字列，然后使用map_df分别对每一个列进行t-test，最后使用broom：tidy将结果整合在tidy的数据框中 。
diff <- data1 %>%select_if(is.numeric) %>%map_df(~ broom::tidy(t.test(. ~ Group,data = https://www.isolves.com/sh/zs/2020-07-23/data1)), .id = 'var')最后对t-test的p值进行校正，保留校正后p值小于0.05的数据 。
diff$p.value <- p.adjust(diff$p.value,"bonferroni")diff <- diff %>% filter(p.value < 0.05) 
秩和检验秩和检验和上面的t-test一样，只需要把代码中的t.test换成wilcox.test就可以了 。
diff1 <- data1 %>%select_if(is.numeric) %>%map_df(~ broom::tidy(wilcox.test(. ~ Group,data = https://www.isolves.com/sh/zs/2020-07-23/data1)), .id = 'var')diff1$p.value <- p.adjust(diff1$p.value,"bonferroni")diff1 <- diff %>% filter(p.value < 0.05) 
数据可视化这个Extended error bar的结果图整体分为两个部分，左侧为组建物种或基因丰度平均值的比较条形图，右侧为组间平均丰度及其95%置信区间的散点图 。
画图的思路是首先分别绘制左右两幅图，之后再拼接在一起，因此需要首先构建绘制两幅图所需的绘图文件 。
 
绘图数据获取对于左侧的组间丰度均值比较条形图，我们首先根据差异性检验的结果从原始的丰度数据文件中提取具有显著差异的列，之后按照分组计算其组内平均丰度，再转换成ggplot绘图所用的长格式数据框 。
abun.bar <- data1[,c(diff$var,"Group")] %>%gather(variable,value,-Group) %>%group_by(variable,Group) %>%summarise(Mean = mean(value))


推荐阅读


上一篇：什么是真正的岩茶 教你如何选购

下一篇：宋文化时期的茶盏,吉州窑