R语言之排序分析

2019-03-23 本文已影响178人 Oodelay

1. 数据准备

rm(list = ls()) #清空当前工作环境
library('vegan')
grp = read.table('grp.txt', row.names = 1, header = T, sep = '\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") #标准化

2. 模型选择

decorana(otu.helli) #执行DCA去趋势分析

根据看分析结果中Axis Lengths的第一轴的大小

如果大于4.0,就应选CCA（基于单峰模型，典范对应分析）
如果在3.0-4.0之间，选RDA和CCA均可
如果小于3.0, RDA的结果会更合理（基于线性模型，冗余分析）

3. 构建模型

fit.0 <- rda (otu.helli ~ 1, data = env) #因变量矩阵和指示变量逐个进行RDA分析
fit.all <- rda (otu.helli ~ ., data = env) #因变量矩阵与所有环境变量进行RDA分析
fit <- ordiR2step (fit.0, scope = formula (fit.all)) #执行前向选择

若是cca()，则使用ordistep()

4. 模型检验

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") # 查看排序轴的显著性

5. 求各个环境因子对整体变异的解释率（贡献率）

基本模型是 rda(x,y,z)或者rda(x~y+condition(z))

x为响应矩阵，通常为群落矩阵，必须有
y为限制性矩阵，通常为环境变量的矩阵，可无，则rda(x)相当于prcomp(x)执行主成分分析
z为条件矩阵，通常也是环境变量，意味着排除（partialed out）该条件矩阵的影响，可无，则为常规rda()或者cca()分析

# 挑选出保留下来的环境因子
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()

严格来说，其中存在包含关系constrained = AN + AP + AK + TN + TP + Joint.
这里使用饼图只是为了方便展示

5. 基本绘图

plot(fit, type="n") #绘制RDA图
text(fit, dis="cn") #绘制RDA图
points(fit, pch=21, col="blue", bg="blue", cex=2) #绘制RDA图
text(fit, "sites", col='red', cex=0.5,  adj = 0.5, pos=2) #绘制RDA图

可以在图窗口点击export将图复制为meta file到Power Point或者Visio中修改

6. 使用ggplot2绘图

#获取样点坐标
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())