微生物环境因子分析(RDA/db-RDA)-ggvegan包
前言
在进行微生物多样性分析时,大家一定会做α,β多样性分析。通俗来讲,α多样性就是样本内的物种多样性。β多样性是指在地区尺度上,物种组成沿着某个梯度方向从一个群落到另一个群落的变化率。即沿着某一环境梯度,物种替代的速率、物种周转率等。
排序的过程是将样品或微生物物种排列在一定的空间, 使得排序轴能够反映一定的生态梯度 这些排序方法又可以分成间接梯度排序(indirect gradient analysis)和直接梯度排序(direct gradient analysis)。间接梯度排序又叫非约束性排序;寻求潜在的或在间接的环境梯度来解释物种数据的变化包括PCA,PCoA,NDMS,直接排序又叫约束性排序;它是指在特定的梯度上(环境轴) 上探讨物种的变化情况;方法包括 RDA, CCA, db-RDA。排序分析(Ordination analysis)。排序(ordination)的过程就是在一个可视化的低维空间或平面重新排列这些样本,使得样本之间的距离最大程度地反映出平面散点图内样本之间的关系信息。
db-RDA 介绍
distance-based redundancy analysis (db-RDA) 是目前在微生物领域应用的最为广泛的环境因子分析,该分析方法内置在R中的vegan包中。相信大家一定都知道vegan包,该R包是进行生态学(包括微生物多样性分析)研究的必备神器!vegan包中提供了所以基本排序分析的方法,可以说是一包在手搞定所有!关于vegan包的详细介绍,请大家查看vegan包的官方文档。本文还会给大家介绍另一款vegan的开挂版本,ggvegan的介绍与使用,ggvegan相当于在vegan软件包中内置了ggplot2, 绘制的图片比用vegan直接绘图更好看!
微生物环境因子分析
要进行微生物环境因子分析,我们需要两个文件,一个是微生物多样性的OTU 表格,另一个就是你所有样品的环境因子数据。比如,你进行土壤微生物研究,这时候你就需要知道你所测土壤的C,N,P,K等化学元素含量以及不同样地的气候信息等等,总之,在分析之前可以多准备些环境因子数据,后期我们还可以对这些环境因子进行共线性,以及环境因子与数据拟合优良性判断。
数据均一化
首先看看我们准备的OTU表格以及环境因子数据结构
(图1):OTU表格
(图2):环境因子数据
读取完数据之后,我们要把OTU的横轴和纵轴调换位置,然后把OTU表格也要进行hellinger转化,使数据均一性更好。并把环境因子进行log转化,以减少同一种环境因子之间本身数值大小造成的影响。
RDA和CCA模型筛选
数据都进行均一化之后,我们要进行RDA和CCA的模型筛选。先用species-sample资料做DCA分析看分析结果中Lengths of gradient的第一轴中的Axis lengths大小,如果大于4.0,就应该选CCA,如果在3.0-4.0之间,选RDA和CCA均可,如果小于3.0,RDA的结果要好于CCA。本例中,我们数据的Axis lengths 大小为0.63859,所有我们应该选择RDA进行分析!
(图3)
差膨胀因子分析
在筛选完RDA和CCA分析后,我们需要利用方差膨胀因子分析,对所有环境因子进行共线性分析。我们要依次删掉最大的变量,也就是删除掉共线性的环境因子,直到所有的变量都小于10。
检测最低AIC值
最后我们要用step模型检测最低AIC值,在这一步中该模型会自动筛选出最优的环境因子。当“none”位于最顶端时意味着改模型筛选结束,位于none值上方的环境因子即为与OTU拟合最好的环境因子。本例中,只有Mg这一个环境因子与我们的OTU拟合的最好。
ANOVA 显著性分析并出图
在进行完以上的数据筛选之后,我们可以用筛选的结果重新进行一次环境因子与OTU的线性回归分析,这样我们就拿到了最终的计算结果,并且用ANOVA进行显著性检验,并且通过该分析我们还可以看到所筛选的环境因子的整体贡献率,以及每个环境因子的单独贡献率。
本例中我们使用了内置ggplot2的vegan—ggvegan进行的分析。ggvegan的出图结果可以用内置的ggplot2进行优化,使你的图更为美观,其具体用法与ggplot2的图层叠加方式类似。详情大家可以参考ggvegan的官网
文中所有测试数据都已放在百度云盘中,请后台回复:“db-RDA”获取。
# 首先要安装devtools包,仅需安装一次 install.packages("devtools") # 加载devtools包 library(devtools) # 下载ggvegan包 devtools::install_github("gavinsimpson/ggvegan")library(ggvegan) otu.tab <- read.csv("otutab.txt", row.names = 1, header=T, sep="\t") env.data <- read.csv("new_meta.txt", row.names = 1, fill = T, header=T, sep="\t") #transform data otu <- t(otu.tab) #data normolization (Legendre and Gallagher,2001) ##by log env.data.log <- log1p(env.data)## ##delete NA env <- na.omit(env.data.log)###hellinger transform otu.hell <- decostand(otu, "hellinger")#DCA analysis sel <- decorana(otu.hell) selotu.tab.0 <- rda(otu.hell ~ 1, env) #no variables #Axis 第一项大于四应该用CCA分析 otu.tab.1<- rda(otu.hell ~ ., env) #我们在筛选完RDA和CCA分析后,我们需要对所有环境因子进行共线性分析,利用方差膨胀因子分析 vif.cca(otu.tab.1) #删除掉共线性的环境因子,删掉最大的变量,直到所有的变量都小于10 otu.tab.1 <- rda(otu.hell ~ N+P+K+Ca+Mg+pH+Al+Fe+Mn+Zn+Mo, env.data.log)vif.cca(otu.tab.1) #进一步筛选 otu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Al+Fe+Mn+Zn+Mo, env.data.log) vif.cca(otu.tab.1) #test again otu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env.data.log)#方差膨胀因子分析,目前所有变量都已经小于10 vif.cca(otu.tab.1) ##用step模型检测最低AIC值 mod.u <- step(otu.tab.0, scope = formula(otu.tab.1), test = "perm")# "perm"增加P值等参数 mod.d <- step(otu.tab.0, scope = (list(lower = formula(otu.tab.0), upper = formula(otu.tab.1)))) mod.d ##本处筛选的结果,找到一个Mg环境因子适合模型构建,为了下一步画图,我们 #保留所有非共线性的环境因子 #choose variables for best model and rda analysis again# (otu.rda.f <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env))anova(otu.rda.f) anova(otu.rda.f, by = "term") anova(otu.rda.f, by = "axis") #计算db-rda otu.tab.bray <- vegdist(otu.hell, "bray")#dbrda() for dbRDA # cca() for CCA otu.tab.b<- capscale(otu.tab.bray ~ ., env) ##绘制db-RDA图 plot(otu.tab.b ,display=c("si","bp","sp"))## 用ggvegan绘制RDA图 p<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), geom = "text", layers = c( "species","sites", "biplot", "centroids"), legend.position = "right", title = "db-RDA") ## 添加图层 p + theme_bw()+theme(panel.grid=element_blank())
Reference
Oksanen, J., Kindt, R., Legendre, P., O’Hara, B., Stevens, M. H. H., Oksanen, M. J., & Suggests, M. A. S. S. (2007). The vegan package. Community ecology package, 10, 631-637.
McArdle, B. H., & Anderson, M. J. (2001). Fitting multivariate models to community data: a comment on distance‐based redundancy analysis. Ecology, 82(1), 290-297.
猜你喜欢
10000+:菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑
系列教程:微生物组入门 Biostar 微生物组 宏基因组
专业技能:学术图表 高分文章 生信宝典 不可或缺的人
一文读懂:宏基因组 寄生虫益处 进化树
必备技能:提问 搜索 Endnote
文献阅读 热心肠 SemanticScholar Geenmedical
扩增子分析:图表解读 分析流程 统计绘图
16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
在线工具:16S预测培养基 生信绘图
科研经验:云笔记 云协作 公众号
编程模板: Shell R Perl
生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘
写在后面
为鼓励读者交流、快速解决科研困难,我们建立了“宏基因组”专业讨论群,目前己有国内外5000+ 一线科研人员加入。参与讨论,获得专业解答,欢迎分享此文至朋友圈,并扫码加主编好友带你入群,务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份,另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍末解决群内讨论,问题不私聊,帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战,关注“宏基因组”
点击阅读原文,跳转最新文章目录阅读
微生物环境因子分析(RDA/db-RDA)-ggvegan包相关推荐
- 用RDA进行微生物环境因子分析
本文首先发布于"宏基因组"公众号原创. 作者:舟行天下 编辑:metagenome 前言 在进行微生物多样性分析时,大家一定会做α,Β多样性分析.α多样品通俗来讲就是样本内的物种多 ...
- 用db-RDA进行微生物环境因子分析-“ggvegan“介绍
前言 在进行微生物多样性分析时,大家一定会做α,β多样性分析.通俗来讲,α多样性就是样本内的物种多样性.β多样性是指在地区尺度上,物种组成沿着某个梯度方向从一个群落到另一个群落的变化率.即沿着某一环境 ...
- oracle+11g+rda,Oracle RDA 4.20 初体验
RDA 全名RemoteDiagnostic Agent,是Oracle用来收集.分析数据库的工具,但统计信息远远大于只是数据库的,也可以说是现在一个Oracle dba 角色需要掌握的Oracle ...
- MySQL rpm包 二进制区别_Linux环境下安装mysql5.6(二进制包不是rpm格式)
一.准备: 1.CentOS release 6.8 2.mysql-5.6.31-linux-glibc2.5-x86_64.tar.gz 3.Linux下MySQL5.6与MySQL5.7安装方法 ...
- Win10 + VS2017 15.5.6 环境下解决 Python 3.6 环境无法刷新DB的问题
Win10 + VS2017 15.5.6 环境下解决 Python 3.6 环境无法刷新DB的问题 参考文章: (1)Win10 + VS2017 15.5.6 环境下解决 Python 3.6 环 ...
- linux打包java jar_在linux环境下修改可运行jar包配置并重新打包
在linux环境下修改可运行jar包配置并重新打包步骤: 1)mkdir xxx 2)mv XXX.jar XXX 3)jar xvf XXX.jar 4)mv XXX.jar ../ 5)vi XX ...
- 在内网环境使用pip离线安装python包
在公司的开发过程中,开发机器或生产机器或许并没有连接外网.这时python的pip和conda等安装方式就废掉了. 我们可以从外网提前下载第三方包,拷贝到内网机器中.而第三方包需要区分不同的运行环境, ...
- 条理清晰的搭建SSH环境之添加所需jar包
一.首先介绍要添加框架环境: JUnit Struts2 Hibernate Spring (1)配置JUnit /**-------------------------添加JUnit-------- ...
- Python环境(基于Pycharm和官方python包)搭建顺序
1.下载安装包 python官网下载3.7.2 Pycharm 社区版下载安装 2.直接使用Pycharm自带virtualEnv File - Settings - Project - Projec ...
最新文章
- linux mipi驱动分析_嵌入式技术在血液分析仪中的应用方案
- windows mysql dump_mysql在Windows下使用mysqldump命令手动备份数据库和自动备份数据库...
- 修改webpack配置,在react中使用less
- Python编程的Turtle 库画出“精美碎花小清新风格树”,速取代码!
- 在 Oracle Enterprise Linux 和 iSCSI 上构建您自己的 Oracle RAC 11g 集群 (2)
- 如何成为一名更出色的开发者?
- [轉]function, new function, new Function
- Linux环境下配置JDK,java环境
- 《东周列国志》第五回 宠虢公周郑交质 助卫逆鲁宋兴兵
- codecademy
- MAC下微信双开(一键命令)
- 试着用人话说说 使命 愿景 价值观,以及人的三观
- 阀门定位器调试顺序详解
- Python制作卡点视频
- 日本知名动画公司东映动画加入 The Sandbox 元宇宙
- 微信小程序---霍兰德职业兴趣测试、心里测评、性格测评
- ggplot2——柱状图
- NKOJ-Unknow 不死的 LYM
- Xanadu:逐梦一个关于光量子计算的“新上都计划”
- 美军派网络部队前往韩国 防“萨德”泄密
热门文章
- 年薪不到 25.2 万免费学,廖雪峰的“大数据高级开发”课程第5期开始招生
- 如何为MNIST手写数字分类开发CNN
- 微软某员工后悔跳槽阿里:工资才多20万不到,天天加班快崩溃!
- 书中自有BAT Offer!
- nodejs 开发企业微信第三方应用入门教程
- 推荐7款实用强大的国产windows软件,你值得拥有!
- 安利7款珍藏已久的windows软件,每一个都非常强大
- 谷歌内部考核制度OKR是怎么样的?你会用OKR吗?
- 高性能存储之--快速理解redis(简版)
- 如何用法向量求点到平面距离_高中数学丨2020新标课本,空间向量与二面角所有知识点,一张表搞定...