數(shù)據(jù)導(dǎo)入
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的行，對列排序遵倦，轉(zhuǎn)置
otu.helli=decostand(otu,method = "hellinger")
標(biāo)準(zhǔn)化
decorana(otu.helli) #執(zhí)行DCA去趨勢分析
fit.0 <- rda (otu.helli ~ 1, data = env) #因變量矩陣和指示變量逐個進(jìn)行RDA分析
fit.all <- rda (otu.helli ~ ., data = env) #因變量矩陣與所有環(huán)境變量進(jìn)行RDA分析
??ordiR2step
fit <- ordiR2step (fit.0, scope = formula (fit.all)) #執(zhí)行前向選擇
vif.cca(fit) # 除去值大于10的環(huán)境因子，因其與其他變量存在高度線性相關(guān)
(total = summary(fit)) # 查看RDA分析結(jié)果
summary(eigenvals(fit)) # 查看每個排序軸的特征值和能夠解釋的變異
查看模型統(tǒng)計(jì)特征
anova.cca(fit) # 查看RDA模型顯著性
anova.cca(fit, by="term", parallel=2) # 查看解釋變量的顯著性
anova.cca(fit, by="axis") # 查看排序軸的顯著性
# 挑選出保留下來的環(huán)境因子
subenv = env[,c("AN", "AP", "TP", "AK", "TN")]
# 新建空矩陣，用來保存結(jié)果
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 # 當(dāng)前變量所能解釋的方差
for (f in names(subenv)){
of = setdiff(names(subenv),f)
temp <- rda(otu.helli, subenv[f], subenv[,of])
temp$CCA$eig -> res[,f] #各個環(huán)境因子單獨(dú)解釋的方差
}
res[,ncol(res)] = total$constr.chi - sum(res[,-c(1:3)]) # 多個環(huán)境因子聯(lián)合解釋的方差
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圖
#獲取樣點(diǎn)坐標(biāo)
sites = total$sites[,1:2] %>% data.frame() %>% merge(grp[,1:2],by = 'row.names')
#獲取指示變量坐標(biāo)
biplot =total$biplot[,1:2] %>% data.frame()
biplot$env = rownames(biplot)
#定義樣點(diǎn)組合在圖例中的出場順序
sites$treatment = factor(sites$treatment, levels = c('S','SH','H'))
#重命名绊诲，須注意與上面的對應(yīng)關(guān)系
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)))+ #增加圖例中點(diǎn)的大小
geom_hline(yintercept=0, linetype=2,color='grey') + #添加經(jīng)過坐標(biāo)原點(diǎn)的橫線
geom_vline(xintercept=0, linetype=2,color='grey') + #添加經(jīng)過坐標(biāo)原點(diǎn)的縱線
geom_point(size = 3) +
stat_ellipse(show.legend = F) + #添加置信區(qū)間
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
結(jié)果展示
image.png
今天跟朋友討論了一下午的R語言，覺著在科研中的作用越來越重要褪贵，但是仔細(xì)想來掂之，在里面投入的時間和精力也不少抗俄，完全失去了科研的本質(zhì)，得不償失世舰，但好歹每天都有一點(diǎn)進(jìn)步动雹。疫情當(dāng)頭，不是我們沖鋒陷陣的時刻跟压，但是也是我們偷偷學(xué)習(xí)胰蝠，努力拔尖的高光時刻。中國加油裆馒，石家莊加油姊氓，自己也要加油。對于未來仍舊很渺茫喷好，仍不想放棄一次發(fā)光的機(jī)會翔横。
今天還是沒能完成規(guī)定的任務(wù)，看來走著走著還是要停下來思考一下梗搅，是不是方向定錯了禾唁，研究的內(nèi)容出了問題或者研究方法落后啥的！
