#usr/bin/env Rscript
1. cca 使用函数预先选择变量
2. cca 首先要要看模型是否显著
#类似地，环境因子（如温度、pH值、湿度等）也可能具有不同的量纲和尺度，例如温度的范围是0到40°C，而湿度的范围是0到100%。
1.如果环境和代谢都已经标准化了，那就不需要在scale=F
## 环境因子矩阵
env <- read.table(in_env,header = T,row.names = 1,fill = TRUE,check.names = FALSE,na.strings = "")
## 分组文件
metadata <- read.table(in_group,sep="\t",header=T, row.names = 1,comment.char = "",skip =0,stringsAsFactors = FALSE, check.names = FALSE)
#metadata <- read.table(in_group,header = T,row.names = 1,comment.char = "",skip =0)
metadata <- metadata[-1,]
#筛选环境因子进行分析,后续再进行筛选
#使用BioENV筛选出与群落数据最大相关的环境因子组合[这里可以设置为选择开关,默认使用所有环境因子作图
bioENV_results <- bioenv(species,env,method="spearman",use="pairwise.complete.obs")
##最佳环境因子组合中可能存在obs_missing现象，需要删除对应样本
env_factors_chr <- summary(bioENV_results)$variables[which(abs(summary(bioENV_results)$correlation) %in% max(abs(summary(bioENV_results)$correlation)))]
env_factors_list <- unlist(strsplit(env_factors_chr," ")) #最佳环境因子组合list
if (length(env_factors_list)==1){
stop("bioENV seleted just one env factor, please manually select more than one env factors for this analysis.")
}else{
sub_env <- env[,which(names(env) %in% env_factors_list)]
obs_missing <- rownames(which(is.na(sub_env),arr.ind = TRUE))
}
library(vegan)
library(ggplot2)
data(varespec)
data(varechem)
#属于预处理，样本的行要一样
{
group <- read_group(group)
colnames(group) <- c("Sample", "Group")
rownames(group) <- group$Sample
sample_id3 <- rownames(group)
sample_id <- intersect(sample_id1, sample_id2) #取交集
sample_id <- intersect(sample_id, sample_id3) #取交集
data1 <- data1[, sample_id] #匹配样本,防止输入的样本数据的顺序不一致
data2 <- data2[, sample_id]
data1 <- t(data1)
data2 <- t(data2)
group <- group[sample_id, ]}
#选择模型 同时说明非线性是否强
dca<- decorana(varespec, iweigh = 0, iresc=4) #列为物种
max_axis = max(dca$rproj[, 1])
#数据标准化
varespec <- decostand(varespec, method = "hellinger")
#run cca
cca_res <- cca(X=varespec, Y=varechem, scale = FALSE) #cca(Y, X, W) Y因变量和X 解释变量
cca_res #快速查看解释率
summary(cca_res)
df_arrows2 <- vegan:::ordiArrowMul(scores(cca_res)$biplot) *scores(cca_res)$biplot %>% as.data.frame()
##c 提取 CCA 典范系数 加权系数得到典范变量
cca_coef <- coef(cca_res)
#r2解释度矫正 能解释多少比例的变异
r2 <-RsquareAdj(cca_res)
cca_noadj <- r2$r.squared #原始 R2
cca_adj <- r2$adj.r.squared#校正后的 R2
cca_exp_adj <- cca_adj * cca_res$CCA$eig/sum(cca_res$CCA$eig) #占比
#otu_cca_eig_adj <- cca_exp_adj * var_cca$tot.chi #变异的多少
#全局检验 所有约束轴的置换检验
otu_cca_test <- vegan::anova.cca(cca_res, permutations = 999)
#逐轴检验，看每条轴是否显著，原则上不显著的轴要舍去
cca_test_axis <- vegan::anova.cca(cca_res, by = 'axis', permutations = 999) #于环境因子的顺序，排在前面的环境因子可能解释度最高且显著
cca_test_axis$`Pr(>F)` <- p.adjust(cca_test_axis$`Pr(>F)`, method = 'bonferroni') #置换检验的P值矫正
fitcca <- envfit(cca_res, varechem, permu=999) %>% head()
df <- as.data.frame(cca_res$CCA$biplot)
df1 <- scores(cca_res)$sites[,1]
df2 <- scores(cca_res)$sites[,2]
df_plot<- data.frame(rownames(scores(cca_res)$sites), df1 , df2)
colnames(df_plot) <- c("sample", "CCA1", "CCA2")
#计算方差膨胀因子
#vif.cca(cca_res)
#前向选择p值
{
#前向选择，这里以 ordiR2step() 的方法为例，基于 999 次置换检验，详情 ?ordiR2step
otu_cca_forward_pr <- ordiR2step(cca(varespec~1, varechem),
scope = formula(cca_res), R2scope = TRUE,
direction = 'forward', permutations = 999)
exp_adj <- RsquareAdj(otu_cca_forward_pr)$adj.r.squared * otu_cca_forward_pr$CCA$eig/sum(otu_cca_forward_pr$CCA$eig)
cca1_exp <- paste('CCA1:', round(exp_adj[1]*100, 2), '%')
cca2_exp <- paste('CCA2:', round(exp_adj[2]*100, 2), '%')
}
site <- scores(cca_res)$sites #样方投影到约束空间+其他空间
spec <- scores(cca_res)$species #物种坐标
envc<- scores(cca_res)$biplot %>% head()
rownames(varechem)
group <- data.frame(group=c(rep("A",6),rep("B",6),
rep("C",6),rep("D",6)))###对样方分组
colnames(grp)="group"###重命名列名
ggplot() +
geom_point(data = site, aes(CCA1, CCA2, color = group$group), size = 4) +
stat_ellipse(data = site, geom = "polygon", level = 0.9, linetype = 1, linewidth = 0.1,
aes(fill = group$group, x = CCA1, y = CCA2),alpha = 0.1, show.legend = T)+
geom_segment(data = spec, aes(x = 0, y = 0, xend = CCA1, yend = CCA2),
arrow =arrow(angle = 22.5, length = unit(0.35, "cm"), type = "closed"), linetype = 1, size = 0.6, colour = "#972b1d") +
geom_segment(data = envc, aes(x = 0, y = 0, xend = CCA1, yend = CCA2),
arrow = arrow(angle = 22.5, length = unit(0.35, "cm"),
type = "closed"),
linetype = 1, size = 0.6, colour = "#e6a84b")
#cca_res$CCA$biplot#输出所有环境因子得分
cca_res$CCA$u %>% head() #样本投影到约束空间
cca_res$CCA$v %>% head()
#otu_cca_forward_pr.scaling1 <- summary(otu_cca_forward_pr, scaling = 1) 标尺 I II型
var_cca_score <- scores(var_cca)
var_cca_score$sites#样方二维坐标值
var_cca$CCA$biplot#特征值
var_cca_score$species#物种坐标值
CCAE <- as.data.frame(var_cca$CCA$biplot[,1:2])#提取特征值前两轴，一般都是二维图形
CCAS1 <- var_cca_score$sites[,1]*0.3#提取样方的X坐标说白了缩放样方坐标距离，防止彼此距离太远没法不方便放在一个平面观察
CCAS2 <- var_cca_score$sites[,2]*0.3#提取样方的Y坐标
plotdata <- data.frame(rownames(var_cca_score$sites), CCAS1, CCAS2)#汇总成为一个数据集
colnames(plotdata) <- c("sample","CCAS1","CCAS2")#更改每列的名字
# 提取CCA得分（样本坐标和环境因子坐标）
cca_scores_sites <- scores(cca_result, display = "sites", scaling = 2)
cca_scores_env <- scores(cca_result, display = "bp", scaling = 2)
# 转换为数据框
cca_scores_sites_df <- as.data.frame(cca_scores_sites)
cca_scores_sites_df$Group <- group_data$Group # 添加分组信息
cca_scores_env_df <- as.data.frame(cca_scores_env) # 环境因子坐标