RDA环境因子群落结构统计检验可视化

环境因子的筛选及数据的转化方面请参阅宏基因组公众号之前的推文，本文主要侧重统计分析与可视化

看到师兄文章里的图自己可能用到，想复现一下，于是就尝试了一下，顺便写个推文记录，在黄静同学的帮助下完成

百度云链接：https://pan.baidu.com/s/1vNmxPV5kw51Ek_oXMmNycw
提取码：no73

因公众号文章不可修改，如以上链接失效，或想获取代码的更新版，请在“宏基因组”公众号后台回复本文关键字“rdaenv”获取最新下载地址。

rm(list=ls())
library(pacman)
p_load(ggplot2,patchwork,vegan,geosphere,psych,corrplot,permute,lattice,ggpubr,RColorBrewer,tidyverse,graphics)
data=read.csv("sum_c.csv",row.names = 1)#读入物种(以Phylum水平为例)矩阵表
head(data,n=3)
env=read.csv("env.csv",row.names = 1)#读入环境因子数据(示例为随机数)
head(env,n=3)
print(decorana(t(data)))
#DCA分析，根据Axis lengths行的第一个值选择排序分析模型
#Axis Lengths >4.0-CCA(基于单峰模型，典范对应分析)；
#如果在3.0-4.0之间-RDA/CCA均可；
#如果小于3.0-RDA(基于线性模型，冗余分析)
B.rda=rda(t(data),env[-1],scale = T)#RDA分析，如果没有环境因子参数，就是PCA分析
#提取样本得分
B.rda.data=data.frame(B.rda$CCA$u[,1:2],env$Treat,rep(c("Mar","Apr","May","Jun","Jul","Aug"),each = 3))#为了仿师兄的图添加的
colnames(B.rda.data)=c("RDA1","RDA2","group","Month")
head(B.rda.data,n=3)
#提取物种得分
B.rda.spe=data.frame(B.rda$CCA$v[,1:2])
B.rda.spe=as.data.frame(B.rda.spe)
B.rda.spe$Species<-rownames(B.rda.spe)
head(B.rda.spe,n=3)
#提取环境因子得分
B.rda.env <- B.rda$CCA$biplot[,1:2]
B.rda.env <- as.data.frame(B.rda.env)
head(B.rda.env,n=3)
#带有环境因子，物种信息且进行不同月份不同处理标记的RDA图(仿师兄)
yanse<-c("darkolivegreen3","gold","dodgerblue4","darkseagreen","chartreuse4","darkorange","burlywood2","brown3","#984EA3","cyan3")
p1=ggplot(data=B.rda.data,aes(RDA1,RDA2))+geom_point(aes(color=Month,fill=Month,shape=group),size=5)+scale_color_manual(values=yanse)+scale_shape_manual(values = c(21,22,23))+scale_fill_manual(values = yanse)+geom_point(data=B.rda.spe,aes(RDA1,RDA2),pch=8,size=5)+geom_text(data=B.rda.spe,aes(x=B.rda.spe[,1],y=B.rda.spe[,2],label=Species),size=5.5,colour="black",hjust=0.5,vjust=1)+labs(title = "B RDA plot",x=paste("RDA1",round(B.rda$CCA$eig[1]/sum(B.rda$CCA$eig)*100,2)," %"),y=paste("RDA2",round(B.rda$CCA$eig[2]/sum(B.rda$CCA$eig)*100,2)," %"))+geom_hline(yintercept = 0,lty=3)+geom_vline(xintercept = 0,lty=3)+geom_segment(data=B.rda.env,aes(x=0,y=0,xend=B.rda.env[,1],yend=B.rda.env[,2]),colour="blue",size=0.8,arrow=arrow(angle = 35,length=unit(0.3,"cm")))+geom_text(data=B.rda.env,aes(x=B.rda.env[,1],y=B.rda.env[,2],label=rownames(B.rda.env)),size=6.5,colour="blue", hjust=(1-sign(B.rda.env[,1]))/2,angle=(180/pi)*atan(B.rda.env[,2]/B.rda.env[,1]))+ggprism::theme_prism()
p1

#统计
B.sum=summary(B.rda)
B.sum$constr.chi/B.sum$tot.chi #constrained表示环境因子对群落结构差异的解释度
B.sum$unconst.chi/B.sum$tot.chi#unconstrained表示环境因子对群落结构不能解释的部分
#环境因子对群落结构差异解释量的饼图绘制
cor_data=data.frame(row.names = c("Explained","Unexplained"),B=c(B.sum$constr.chi/B.sum$tot.chi,B.sum$unconst.chi/B.sum$tot.chi))
cor_data$group=rownames(cor_data)
head(cor_data,n=3)
cor_data <- data.frame(cor_data)
cor_data=arrange(cor_data,B)
head(cor_data,n=3)
labs<-paste0(cor_data$group,"\n(",round(cor_data$B/sum(cor_data$B)*100,2),"%)")
pie(cor_data$B,labels=labs,init.angle = 90,col=brewer.pal(nrow(cor_data),"Reds"),boder="black")
#在R中手动导出，右侧出图区Export-PDF

#anova.cca检验
B.perm=permutest(B.rda,permu=999) # permu=999是表示置换999次
B.perm #是做环境因子整体与群落结构差异的相关性(解释量)，anova.cca {vegan}
#envfit检验   envfit函数跟mantel(见下文)的功能是一样的
B.ef=envfit(B.rda,env[-1],permu=999)#是做每一个环境因子与群落结构差异的相关性(解释量)
B.ef$vectors$r#R2值
B.ef$vectors$pvals#P值
#每一个环境因子对群落结构差异解释量的柱形图绘制的数据整理
cor_com=data.frame(tax=rownames(B.rda.env),B.r=B.ef$vectors$r,B.p=B.ef$vectors$pvals)
cor_com=arrange(cor_com,B.r)
head(cor_com,n=3)
cor_com[c(3)]=cor_com[c(3)]>0.05
head(cor_com,n=3)
#将p<0.05标记为FALSE，p>0.05标记为TRUE，使用此数据绘制柱形图，将其可视化
#envfit检验可视化
cor_com$tax = factor(cor_com$tax,order = T,levels = row.names(cor_com))#按R2值排序
p2 <- ggplot(cor_com, aes(x =tax, y = B.r),size=2) +geom_bar(stat = 'identity', width = 0.8,color="black",fill="red") +scale_fill_manual(guide = FALSE)+geom_text(aes(y = B.r+0.005, label = ifelse(B.p==TRUE,"","*")),size = 5, fontface = "bold") +xlab("Environmental factor")+ylab(expression(r^"2"))+scale_y_continuous(expand = c(0,0))+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p2

#mantel检验
data <- as.data.frame(t(data))
species.distance<-vegdist(data,method = 'bray')
soil <- NULL
for (i in 2:ncol(env)) {dd <- mantel(species.distance, vegdist(env[,i], method = "euclidean"), method = "pearson", permutations = 9999, na.rm = TRUE)soil <- rbind(soil,c(colnames(env)[i],dd$statistic, dd$signif))
}
head(soil,n=3)
soil <- data.frame(aa=rownames(B.rda.env),M_r=soil[,2],M.p=soil[,3])
rownames(soil)=soil$aa
soil=arrange(soil,M_r)
soil[c(3)]=soil[c(3)]>0.05 # 将p<0.05标记为FALSE，p>0.05标记为TRUE，同上
soil$aa = factor(soil$aa,order = T,levels = row.names(soil))
soil$M_r=round(as.numeric(soil$M_r),4)
head(soil,n=3)
#mantel检验可视化
p3 <- ggplot(soil, aes(x =aa, y = M_r),size=2) +geom_bar(stat = 'identity', width = 0.8,color="black",fill="red") +scale_fill_manual(guide = FALSE)+geom_text(aes(y = M_r+0.005, label = ifelse(M.p==TRUE,"","*")),size = 5, fontface = "bold") +xlab("Environmental factor")+ylab(expression(r))+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p3

##################群落结构差异的统计检验(三种方法及可视化)
####Adonis
otu <- data.frame(data)#mantel检验时已经转置，勿要转置
head(otu)
#样本分组文件
group <- read.delim('metadata.txt', sep = '\t', stringsAsFactors = FALSE)
head(group)
#使用 Bray-Curtis 距离测度  unifrac_binary
adonis_result <- adonis(otu~Group, group, distance = 'Bray-Curtis', permutations = 999)
adonis_result$aov.tab
group_name <- unique(group$Group)adonis_result_two <- NULL
for (i in 1:(length(group_name) - 1)) {for (j in (i + 1):length(group_name)) {group_ij <- subset(group, Group %in% c(group_name[i], group_name[j]))otu_ij <- otu[group_ij$SampleID, ]adonis_result_otu_ij <- adonis(otu_ij~Group, group_ij, permutations = 999, distance = 'bray')adonis_result_two <- rbind(adonis_result_two, c(paste(group_name[i], group_name[j], sep = '_'),       'Bray-Curtis', unlist(data.frame(adonis_result_otu_ij$aov.tab, check.names = FALSE)[1, ])))}
}
adonis_result_two <- data.frame(adonis_result_two, stringsAsFactors = FALSE)
names(adonis_result_two) <- c('group', 'distance', 'Df', 'Sums of squares', 'Mean squares', 'F.Model', 'R2', 'P')
adonis_result_two$R2<- as.numeric(adonis_result_two$R2)
adonis_result_two$P <- as.numeric(adonis_result_two$P)
#p值Benjamini校正 作图时自己选择用哪一个P值(下同)
adonis_result_two$P_adj_BH <- p.adjust(adonis_result_two$'P', method = 'BH')
head(adonis_result_two)
adonis_result_two=arrange(adonis_result_two,R2)
adonis_result_two
adonis_result_two$P=adonis_result_two$P>0.05#将p<0.05标记为FALSE，p>0.05标记为TRUE，使用此数据绘制柱形图(下同)
adonis_result_two$tax = factor(adonis_result_two$group,order = T,levels = adonis_result_two$group)
adonis <- ggplot(adonis_result_two, aes(x =tax, y = R2),size=2) +geom_bar(stat = 'identity', width = 0.8,color="black",fill="red")+scale_fill_manual(guide = FALSE)+geom_text(aes(y = R2-0.02, label = ifelse(P==TRUE,"","*")),#选择作图的P值(P_adj_BH)(下同)size = 5, fontface = "bold") +xlab("")+ylab(expression(r^"2"))+scale_y_continuous(expand = c(0,0))+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 0))
adonis ##图片大小在保存时自己调合适

####Anosim
anosim_result <- anosim(otu, group$Group, distance = 'bray', permutations = 999)
anosim_result$signif    #p 值
anosim_result$statistic    #R 值
group_name <- unique(group$Group)
anosim_result_two <- NULL
for (i in 1:(length(group_name) - 1)) {for (j in (i + 1):length(group_name)) {group_ij <- subset(group, Group %in% c(group_name[i], group_name[j]))otu_ij <- otu[group_ij$SampleID,]anosim_result_otu_ij <- anosim(otu_ij, group_ij$Group, permutations = 999, distance = 'bray')anosim_result_two <- rbind(anosim_result_two, c(paste(group_name[i], group_name[j], sep =        '_'), 'Bray-Curtis', anosim_result_otu_ij$statistic, anosim_result_otu_ij$signif))}
}
anosim_result_two <- data.frame(anosim_result_two, stringsAsFactors = FALSE)
names(anosim_result_two) <- c('group', 'distance', 'R', 'P')
anosim_result_two$R<- as.numeric(anosim_result_two$R)
anosim_result_two$P <- as.numeric(anosim_result_two$P)
anosim_result_two$P_adj_BH <- p.adjust(anosim_result_two$P, method = 'BH')
head(anosim_result_two)anosim_result_two=arrange(anosim_result_two,R)
anosim_result_two
anosim_result_two$P=anosim_result_two$P>0.05
anosim_result_two$tax = factor(anosim_result_two$group,order = T,levels = anosim_result_two$group)
anosim <- ggplot(anosim_result_two, aes(x =tax, y = R),size=2) +geom_bar(stat = 'identity', width = 0.8,color="black",fill="red")+scale_fill_manual(guide = FALSE)+geom_text(aes(y =R-0.02, label = ifelse(P==TRUE,"","*")),size = 5, fontface = "bold") +xlab("")+ylab(expression(r))+scale_y_continuous(expand = c(0,0))+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 0))
anosim

##MRPP
mrpp_result <- mrpp(otu, group$Group, distance = 'bray', permutations = 999)
mrpp_result$Pvalue    # p 值
mrpp_result$A  #相当于Anosim检验的R值
roup_name <- unique(group$Group)
mrpp_result_two <- NULL
for (i in 1:(length(group_name) - 1)) {for (j in (i + 1):length(group_name)) {group_ij <- subset(group, Group %in% c(group_name[i], group_name[j]))otu_ij <- otu[group_ij$SampleID,]mrpp_result_otu_ij <- mrpp(otu_ij, group_ij$Group, permutations = 999, distance = 'bray')    mrpp_result_two <- rbind(mrpp_result_two, c(paste(group_name[i], group_name[j], sep = '_'),      'Bray-Curtis', mrpp_result_otu_ij$A, mrpp_result_otu_ij$delta, mrpp_result_otu_ij$E.delta,       mrpp_result_otu_ij$Pvalue))}
}
mrpp_result_two <- data.frame(mrpp_result_two, stringsAsFactors = FALSE)
names(mrpp_result_two) <- c('group', 'distance', 'A', 'Observe_delta', 'Expect_delta', 'P')
mrpp_result_two$A<- as.numeric(mrpp_result_two$A)
mrpp_result_two$P <- as.numeric(mrpp_result_two$P)
mrpp_result_two$P_adj_BH <- p.adjust(mrpp_result_two$P, method = 'BH')
head(mrpp_result_two)
mrpp_result_two=arrange(mrpp_result_two,A)
mrpp_result_two
mrpp_result_two$P=mrpp_result_two$P>0.05
mrpp_result_two$tax = factor(mrpp_result_two$group,order = T,levels = mrpp_result_two$group)
MRPP <- ggplot(mrpp_result_two, aes(x =tax, y = A),size=2) +geom_bar(stat = 'identity', width = 0.8,color="black",fill="red")+scale_fill_manual(guide = FALSE)+geom_text(aes(y =A-0.02, label = ifelse(P==TRUE,"","*")),size = 5, fontface = "bold") +xlab("")+ylab(expression(r))+scale_y_continuous(expand = c(0,0))+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 0))
MRPP # 其余数据集重复绘制

#Output figure width and height
# Letter纸图片尺寸为单栏89 mm，双栏183 mm，页面最宽为247 mm 推荐比例16：10，
# 即半版89 mm x 56 mm; 183 mm x 114 mm
#
##################保存ggsave("./p1.pdf", p1, width = 183, height = 114, units = "mm")ggsave("./p2.pdf", p2, width = 183, height = 114, units = "mm")ggsave("./p3.pdf", p3, width = 183, height = 114, units = "mm")ggsave("./p4.pdf", adonis, width = 183, height = 114, units = "mm")ggsave("./p5.pdf", anosim, width = 183, height = 114, units = "mm")ggsave("./p6.pdf", MRPP, width = 183, height = 114, units = "mm")

参考资料

R绘图-RDA排序分析

R包vegan的置换多元方差分析（PERMANOVA）判断群落结构差异

R包vegan的相似性分析（ANOSIM）判断群落结构差异

R包vegan的MRPP分析判断群落结构差异

作者：黄静、旭日阳光

责编：马腾飞 南京农业大学

审核：刘永鑫 中科院遗传发育所

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发NatureCell专刊肠道指挥大脑

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命生命大跃进  细胞暗战 人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读

RDA_环境因子_群落结构_统计检验_可视化相关推荐

  1. C++_结构体中const使用场景_结构体_毕业设计案例_使用结构体数组_随机数种子---C++语言工作笔记027

    然后我们再看const符号,在结构体中的应用 首先我们新建一个结构体.student 然后定义一个结构体变量,并初始化 然后我们再写个方法,去打印这个结构体变量 可以看到我们用的是传值

  2. C++_选择结构_循环结构_for循环_敲桌子案例_嵌套循环_乘法口诀案例_跳转语句break---C++语言工作笔记018

    跟java ,一模一样 这样写也可以,带劲

  3. C++_选择结构_switch语句_循环结构while_while案例猜数字_do while循环_dowhile案例水仙花数_---C++语言工作笔记017

    跟java一样啊

  4. var _ 接口 = 结构体{}

    gin源码ps:把问题暴露在编译阶段,实例化Engine结构体,并立马丢掉,确保结构体实现了IRouter接口var _ 接口 = &结构体{}

  5. rda分析怎么做_群落分析的冗余分析（RDA）概述

    约束排序之冗余分析(RDA)概述 前篇先后简介了主成分分析(PCA).对应分析(CA).主坐标分析(PCoA)以及非度量多维尺度分析(NMDS).这些排序方法均属于非约束排序,只涉及一个数据矩阵,并在 ...

  6. C语言_函数结构体的调用

    C语言_函数结构体的调用 #include<stdio.h> //定义存储函数的结构体 struct map{//定义无参数类型返回void的函数指针void (*p)(); }; /** ...

  7. Java基础语法_循环结构【多测师_何sir】

    Java基础语法_循环结构 for 循环 while 循环 do-while 循环 Java 增强 for 循环 break 关键字 continue 关键字 for 循环 语法结构: for(初始化 ...

  8. Algorithms_基础数据结构(04)_线性表之链表_单向循环链表约瑟夫环问题

    文章目录 大纲图 链表的经典面试题目 如何设计一个LRU缓存淘汰算法 约瑟夫问题 结构 分析 大纲图 链表的经典面试题目 如何设计一个LRU缓存淘汰算法 tip:单向链表 约瑟夫问题 N个人围成一圈, ...

  9. c++语言编程,一个电灯两个开关控制,[理学]四川大学计算机学院精品课程_面向对象程序设计C++课件_游洪越_第一章绪论.ppt...

    [理学]四川大学计算机学院精品课程_面向对象程序设计C课件_游洪越_第一章绪论 主讲教师: 游洪跃 个人主页: /~youhongyue 邮件地址: youhongyao@ 教材:<C++面向对 ...

最新文章

  1. 微信8.0内测更新！！！（附内测体验资格）
  2. mysql $的用法_MYSQL limit用法
  3. win32汇编获取当前进程ID和可执行文件名
  4. Linux netstat查看网络连接信息
  5. 龙蜥社区成立系统运维SIG，开源sysAK系统运维工具集
  6. MATLAB基础教程（6）——使用matlab求解线性方程组
  7. 从潘叔到潘子，潘长江走下「神坛」
  8. 利用Session实现一次性验证码（多学一招）
  9. PLTS中计算Skew(计算延时差：对内/对间)
  10. FZU Problem 2198 快来快来数一数(矩阵快速幂 卡常数 +优化)
  11. 如何快速取消大量的合并单元格并向下填充数据
  12. 支付宝相关服务申请入口
  13. UNIX Time Sharing System - UNIX分时系统翻译
  14. ubuntu 内存占用过高导致卡死 解决办法
  15. JAVA-----乱码的处理 乱码的解决方法总结
  16. Matlab里for循环详解
  17. 内控遇到的问题及解决方法
  18. 加密狗软加密方案离线绑定与解绑
  19. 关于联通手机营业厅中的在线销户，大家有什么想说的？
  20. stdcall和cdecl

热门文章

  1. 《大话设计模式》读书总结
  2. 【Android 】零基础到飞升 | AbsoluteLayout(绝对布局)
  3. 第二证券|小鹏持续萎靡，理想蔚来逆势反弹破月销记录
  4. options should NOT have additional properties
  5. ZX297520V3T：Codec NAU88C22驱动调试
  6. 出海的成本越来越高，奈何
  7. 字典型列表转化为csv文件格式
  8. cxf 实名认证---全国公民身份信息系统
  9. firefox 14 vim化——Pentadactyl
  10. asp.net2.0中的ValidationGroup (转)