HarmonyOS鸿蒙最新微生物环境因子分析(RDA db-RDA)-ggvegan包(5)，面试须知

本文链接：https://blog.csdn.net/2401_84871031/article/details/138984350

深知大多数程序员，想要提升技能，往往是自己摸索成长，但自己不成体系的自学效果低效又漫长，而且极易碰到天花板技术停滞不前！

既有适合小白学习的零基础资料，也有适合3年以上经验的小伙伴深入学习提升的进阶课程，涵盖了95%以上鸿蒙开发知识点，真正体系化！

由于文件比较多，这里只是将部分目录截图出来，全套包含大厂面经、学习笔记、源码讲义、实战项目、大纲路线、讲解视频，并且后续会持续更新

读取完数据之后，我们要把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)
sel

otu.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 微生物组宏基因组

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