利用R语言做冗余分析

在统计学中,冗余分析是通过原始变量与典型变量之间的相关性,分析引起原始变量变异的原因。以原始变量为因变量,典型变量为自变量,建立线性回归模型,则相应的确定系数等于因变量与典型变量间相关系数的平方。它描述了由于因变量和典型变量的线性关系引起的因变量变异在因变量总变异中的比例。 在生态学中应用比较多,也有用来做基因型与表型性状做关联分析的。利用R语言比较好实现,也方便做图。

主要需要用到vegan这个R包,其中比较关键的函数是rda这个函数,进行数据分析后可以用ggplot2进行绘图。具体代码如下:

library(vegan)
library(ggrepel)  ###这个包主要用来防止画图时标签重叠,不过好像用处也不是很大,不用的话将后面geom_text_repel改成geom_text就可以了,标签位置可以后期用AI进行调整
library(ggplot2)
scale_df <-function(xv)          ########数据标准化
{
  max_v <- max(xv)
  min_v <- min(xv)
  return(round((xv-min_v)/(max_v-min_v), 2))
}
df <- read.table("df.txt", sep = "\t", header = F)
df <- apply(df, 2, scale_df)
env <- as.data.frame(df[,-c(1:7)])    #######数据的第8列到最后为解释变量,主要是一些环境因子
phe <- as.data.frame(df[,c(1:7)])    ####数据的前7列为响应变量,主要是样本的表型数据
RDA <- rda(env~phe[,1] + phe[,2] + phe[,3] + phe[,4] + phe[,5] + phe[,6] + phe[,7], phe)  ######进行RDA分析
######################################################
ii <- summary(RDA)
sp <- as.data.frame(ii$species[,1:2]) * 2      ####为了避免向量的标签重叠将数据进行了变换,不会影响向量的形状,对结果不会有影响
st <- as.data.frame(ii$sites[,1:2])
yz <- as.data.frame(ii$biplot[,1:2]) * 3
row.names(sp) <- c("MWD", "WHC", "BD", "SP", "CEC", "pH", "S1", "S2", "S3")   ####对解释变量进行重命名
row.names(yz) <- c("SOC", "TN", "SOCS", "TNS", "POC", "AN", "POC/MOC")  ####对响应变量进行重命名
p <- ggplot() + 
  geom_point(data = st, aes(RDA1, RDA2), col = "gray") +  ###各样本的分布图,也可以去掉
  geom_segment(data = sp, aes(x=0, y=0, xend = RDA1, yend = RDA2),
               arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                             type = "closed"),linetype=1, size=0.6,colour = "red") +  ####解释变量的向量图
  geom_text_repel(data = sp, aes(RDA1, RDA2, label=row.names(sp))) +
  geom_segment(data = yz, aes(x=0, y=0, xend = RDA1, yend = RDA2),
               arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                             type = "closed"),linetype=1, size=0.6,colour = "blue") +    ####响应变量的向量图
  geom_text_repel(data = yz, aes(RDA1, RDA2, label=row.names(yz))) +
  labs(x="RDA axis1(76.89%)", y="RDA axis2(14.43%)") +
  geom_hline(yintercept=0,linetype=3, size=1) +   ####分割线
  geom_vline(xintercept=0,linetype=3, size=1)+
  theme_bw()+theme(panel.grid=element_blank())
pdf("RDA.pdf")
print(p)
dev.off()

主要参考资料:

冗余分析(RDA)在R语言中的实现
R语言做冗余分析(RDA)的一个简单小例子
R语言实现冗余分析(RDA)完整代码
第五章 非约束排序3 主成分分析(PCA)-基于Hellinger转化的数据以及其他内容