数据导入
library('vegan')
grp = read.table('grp.txt', row.names = 1, header = T)
env = grp[,-c(1:2)]
otu=read.table("OTU.txt",header=T,row.names = 1, sep = '\t') #读取因变量矩阵
library('dplyr')
otu = otu[,rownames(grp)] %>% .[rowSums(.) !=0,] %>% t() #去除均为0的行，对列排序，转置
otu.helli=decostand(otu,method = "hellinger")
标准化
decorana(otu.helli) #执行DCA去趋势分析
fit.0 <- rda (otu.helli ~ 1, data = env) #因变量矩阵和指示变量逐个进行RDA分析
fit.all <- rda (otu.helli ~ ., data = env) #因变量矩阵与所有环境变量进行RDA分析
??ordiR2step
fit <- ordiR2step (fit.0, scope = formula (fit.all)) #执行前向选择
vif.cca(fit) # 除去值大于10的环境因子，因其与其他变量存在高度线性相关
(total = summary(fit)) # 查看RDA分析结果
summary(eigenvals(fit)) # 查看每个排序轴的特征值和能够解释的变异
查看模型统计特征
anova.cca(fit) # 查看RDA模型显著性
anova.cca(fit, by="term", parallel=2) # 查看解释变量的显著性
anova.cca(fit, by="axis") # 查看排序轴的显著性
# 挑选出保留下来的环境因子
subenv = env[,c("AN", "AP", "TP", "AK", "TN")]
# 新建空矩阵，用来保存结果
res = matrix(0,nrow = 1, ncol = ncol(subenv)+4,
dimnames = list(c('Inertia'),c('total','Unconstrained','constrained',names(subenv),'joint')))
res[,1] = total$tot.chi #总体方差
res[,2] = total$unconst.chi # 未解释方差
res[,3] = total$constr.chi # 当前变量所能解释的方差
for (f in names(subenv)){
of = setdiff(names(subenv),f)
temp <- rda(otu.helli, subenv[f], subenv[,of])
temp$CCA$eig -> res[,f] #各个环境因子单独解释的方差
}
res[,ncol(res)] = total$constr.chi - sum(res[,-c(1:3)]) # 多个环境因子联合解释的方差
res = res/res[,1]
饼图展示效果
# jpeg('pie.jpg',width = 2000, height = 2000, res = 300)
pie(res[,-1],init.angle = 60,radius = 0.5,
labels = paste(colnames(res)[-1]," ", round((as.numeric(res[,-1])),3)*100,'%'),
main = "Variance explained by each factor",
col=rainbow(8))
# dev.off()
plot(fit, type="n") #绘制RDA图
text(fit, dis="cn") #绘制RDA图
points(fit, pch=21, col="red", bg="pink", cex=1) #绘制RDA图
text(fit, "sites", col='red', cex=0.5, adj = 0.5, pos=2) #绘制RDA图
#获取样点坐标
sites = total$sites[,1:2] %>% data.frame() %>% merge(grp[,1:2],by = 'row.names')
#获取指示变量坐标
biplot =total$biplot[,1:2] %>% data.frame()
biplot$env = rownames(biplot)
#定义样点组合在图例中的出场顺序
sites$treatment = factor(sites$treatment, levels = c('S','SH','H'))
#重命名，须注意与上面的对应关系
levels(sites$treatment) = c('水田','水旱轮作','旱田')
library('ggplot2')
library('ggrepel')
p <- ggplot(sites, aes(x = RDA1, y =RDA2, color = treatment)) +
theme_bw() +
labs( x= "RDA1 (23.42%)", y = "RDA2 (11.84%)", color = '') +
guides(color = guide_legend(override.aes = list(size=5)))+ #增加图例中点的大小
geom_hline(yintercept=0, linetype=2,color='grey') + #添加经过坐标原点的横线
geom_vline(xintercept=0, linetype=2,color='grey') + #添加经过坐标原点的纵线
geom_point(size = 3) +
stat_ellipse(show.legend = F) + #添加置信区间
geom_segment(data = biplot,
aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(length = unit(1/2, 'picas')), lwd = 1,
colour = "blue") + #添加箭头
geom_text_repel(data = biplot,
aes(x=RDA1,y=RDA2,label=env),
size= 5, fontface='bold',color='black')+ #添加指示变量文本
theme(legend.position = c(0.9,0.2),
legend.background = element_blank(),
legend.text = element_text(face = 'bold',color='black',size=12),
axis.title = element_text(face = 'bold',color='black',size=14),
axis.text = element_text(face = 'bold',color='black',size=12),
panel.grid = element_blank())
p
结果展示
image.png
今天跟朋友讨论了一下午的R语言，觉着在科研中的作用越来越重要，但是仔细想来，在里面投入的时间和精力也不少，完全失去了科研的本质，得不偿失，但好歹每天都有一点进步。疫情当头，不是我们冲锋陷阵的时刻，但是也是我们偷偷学习，努力拔尖的高光时刻。中国加油，石家庄加油，自己也要加油。对于未来仍旧很渺茫，仍不想放弃一次发光的机会。
今天还是没能完成规定的任务，看来走着走着还是要停下来思考一下，是不是方向定错了，研究的内容出了问题或者研究方法落后啥的！
