环境因子分析
---
title: "环境因子分析"
author: "Omicsgene"
date: "2020年11月17日"
output: html_document
---
 | |
---|---|---
问答社区：http://www.omicsclass.com/ | 组学大讲堂公众号|
---
||
---|---|---
课程推荐1：[R语言入门与基础绘图](https://bdtcd.xetslk.com/s/2G8tHr)|课程推荐2：[R语言绘图（ggplot）](https://bdtcd.xetslk.com/s/2G8tHr)|所有生信课程:[点击](https://study.omicsclass.com/index)
---
## 1.加载需要的包
```{r warning=FALSE}
# 2.1 安装CRAN来源常用包
#设置镜像，
local({r <- getOption("repos")
r["CRAN"] <- "http://mirrors.tuna.tsinghua.edu.cn/CRAN/"
options(repos=r)})
# 依赖包列表：
package_list <- c("grid","ggplot2","gridExtra","vegan","reshape2","RColorBrewer","ggrepel","corrplot","tidyverse")
# 判断R包加载是否成功如果加载不成功自动安装
for(p in package_list){
if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){
install.packages(p, warn.conflicts = FALSE)
suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))
}
}
# 2.2 安装bioconductor常用包
# options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
# package_list <- c("digest","AnnotationDbi", "impute", "GO.db", "preprocessCore","WGCNA","multtest")
# for(p in package_list){
# if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){
# if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
# BiocManager::install(p)
# suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))
# }
# }
```
## 2.读入数据并做处理
```{r }
#设置工作目录
#setwd("D:\\网络培训\\微生物扩增子分析实操\\ENV环境因子分析")
#读取数据
env=read.table("data/env.txt",header=T,sep="\t",row.names = 1,comment.char = "",stringsAsFactors = F)
otu=read.table("data/feature_table.txt",header=T,sep="\t",row.names = 1,comment.char = "")
fastmap=read.table("data/fastmap.txt",header=T,sep="\t",row.names = 1,comment.char = "")
#去除含有缺失数据的样本
env = na.omit(env)
#物种丰度过滤,物种在所有样本中的丰度值要大于样本数(可选过滤)
#otu = otu[,colSums(otu) > ncol(otu)]
#转置数据：vegan 要求行为样本
otu=t(otu)
```
数据抽平标准化 vegan包中：rrarefy，按照最小值进行抽平
数据的前期处理已经抽平，这一步选做
```{r}
range(rowSums(otu))
otu = as.data.frame(rrarefy(otu, range(rowSums(otu))[1])) #按最小值抽平
range(rowSums(otu))
```
筛选样本保持一致:fastmap otu env 中样本统一
```{r}
inter=intersect(rownames(env),rownames(otu))
env=env[inter,]
otu=otu[inter,]
fastmap=fastmap[inter,]
head(fastmap)
head(env)
# 保存筛选样本的抽平后的OTU表，后续绘图会用到计算richness
otu_count = otu
```
## 2.1 数据转换选择
otu table Hellinger 转换选择：如果转化，后面分析RDA时即为tb-RDA分析：
Hellinger预转化（处理包含很多0值的群落物种数据时，推荐使用）
```{r}
otu = decostand(otu, "hellinger")
# 将环境因子进行log+1自然对数转化，使得环境因子数据更均一
env = log1p(env)
```
## 3 环境因子之间相关性分析(pearson)
```{r}
#相关性分析
env_corr=cor(env,method="pearson") #c("pearson", "kendall", "spearman")
#设置色板
palette_1 <- RColorBrewer::brewer.pal(n=11, name = "RdYlBu")
palette_2 <- rev(palette_1) # 色板反向
#绘图
corrplot(env_corr, method = "square", type = "lower", order="AOE",
col = palette_2, tl.pos = "tp", tl.col = "blue", cl.pos = "r",
title = "env corr", mar = c(2,4,4,1))
corrplot(env_corr, method = "number", type = "upper",order="AOE",
col = palette_2,add = TRUE, diag = FALSE, tl.pos = "n", cl.pos = "n",
number.cex=0.5,mar = c(2,4,4,1))
```
加入参数order=”AOE”，使相关性高的一类因子放在一起。另外”FPC”，”hclust”也有类似的效果。
Hb, hemoglobin; ALB, albumin; BUN, blood urea nitrogen; SCr, serum creatinine; eGFR, estimated glomerular filtration rate; P, phosphate.
Hb，血红蛋白； ALB，白蛋白； BUN，血液尿素氮； SCr，血清肌酐； eGFR，估计的肾小球滤过率； P，磷酸盐。
### 相关性图输出
```{r}
pdf(file="env_cor.pdf",w=10,h=10)
corrplot(env_corr, method = "square", type = "lower", order="AOE",
col = palette_2, tl.pos = "tp", tl.col = "blue", cl.pos = "r",
title = "env corr", mar = c(2,4,4,1))
corrplot(env_corr, method = "number", type = "upper",order="AOE",
col = palette_2,add = TRUE, diag = FALSE, tl.pos = "n", cl.pos = "n",
number.cex=0.5,mar = c(2,4,4,1))
dev.off()
```
## 4 环境因子与群落相关性分析
### 4.1 曼特尔检验Mantel test
Mantel test由Nathan Mantel在1976年提出，在微生物研究上有着广泛的应用，是对两个矩阵相关关系的非参数统计方法。
通常将微生物群落作为一个距离矩阵（如bray distance matrix），环境变量作为另一个距离矩阵（如pH、温度、地理、体重、BMI值、血糖、免疫因子等），再检验两个矩阵之间的相关性。
另外还有：Partial Mantel test 偏曼特尔检验分析，可以将某些环境因子作为控制矩阵，计算微生物群落矩阵与剩余目标环境因子矩阵的相关性。即排除不关注的环境因子，研究关注的环境因子与微生物之间的关系。
```{r}
#otutable转换成bray 矩阵
otu_bray<-vegdist(otu, method="bray")
#eGFR 转换矩阵
env_eGFR<-vegdist(env$eGFR, method="euclidean")
##BUN 转换矩阵
env_BUN<-vegdist(env$BUN, method="euclidean")
##Mantel test相关性分析，不考虑交互作用
mantel(otu_bray, env_eGFR, permutations=999)
mantel(otu_bray, env_BUN, permutations=999)
##相关性分析，考虑交互作用
mantel(otu_bray, env_eGFR+env_BUN, permutations=999)
##mantel.partial 当这些环境变量间存在自相关关系时，需要控制其他变量
mantel.partial(otu_bray, env_eGFR, env_BUN, permutations=999)
##批量mantel分析函数
mt<-function(otu,env){
vars <- colnames(env)
models<-list()
for (i in seq_along(vars)){
otu_bray<-vegdist(otu,method = "bray")
env_dis<-vegdist(env[vars[i]],method = "euclidean")
model <- mantel(otu_bray,env_dis, permutations=999)
name <- vars[i]
statistic <- model$statistic
signif <- model$signif
models[[i]] <- data.frame(name = name, statistic = statistic, signif = signif, row.names = NULL)
}
models %>% bind_rows()
}
mantel_result=mt(otu,env)
mantel_result
```
- 第1列为环境因子名称。
- 第2列statistic(r值)统计量，statistic大于0表示正相关，statistic小于0表示负相关，statistic值的绝对值越大则越相关。
- 第3列signif表示显著性水平，小于0.05表示具有显著性相关。
## 4.2 环境因子与群落多样性指数相关性分析（拟合展示）
```{r}
richness = as.data.frame(specnumber(otu_count))
colnames(richness) = "richness"
#或者可计算shannon指数
#alpha_shannon <- as.data.frame(diversity(otu, "shannon") )
#colnames(alpha_shannon)="shannon"
#head(alpha_shannon)
richness.lm = lm(richness$richness ~ ., data = env)
summary(richness.lm)
######绘图展示################################
#添加患者stage分组
df = cbind(group=fastmap$group, richness, env)
head(df,5)
#数据转换成长形数据方便后续绘图
df_all = melt(df, id=1:2)
head(df_all)
# 绘制所有环境因子
# 保留前4列不变
p=ggplot(df_all,aes(value, richness)) +
geom_point(aes(colour=group)) +
geom_smooth(method="lm", size=1, se=T) +
facet_wrap( ~ variable , scales="free", ncol=4) +theme_bw()
p
```
## 5 环境因子解释度分析
找到影响微生物群落显著的环境因素。
### adonis 分析筛选
Adonis又称置换多元方差分析或非参数多元方差分析（PERMANOVA）。它利用半度量（如Bray-Curtis）或度量的距离矩阵（如Euclidean）对总方差进行分解，分析不同分组因素对样本差异的解释度，并使用置换检验对其统计学意义进行显著性分析。
也可以衡量 环境因子（连续变量）解释度
```{r}
# adonis统计otu与所有环境因子关联
# 本方法计算的结果与dbrda和capscale中的anova.cca方法得到的结果一样。
##对分组解释度检验
otu_group.adonis = adonis(vegdist(otu,method="bray") ~group, data=fastmap)
otu_group.adonis
#多个
otu_group.adonis = adonis(vegdist(otu,method="bray") ~age+group+gender, data=fastmap)
otu_group.adonis
#考虑交互
otu_group.adonis = adonis(vegdist(otu,method="bray") ~age+group*gender, data=fastmap)
otu_group.adonis
##对数量型环境因子解释度检验
otu_env.adonis = adonis(vegdist(otu,method="bray") ~., data=env)
otu_env.adonis
#想知道两种环境因子交互作用
otu_env.adonis = adonis(vegdist(otu,method="bray") ~eGFR*BUN, data=env)
otu_env.adonis
otu_env.adonis = adonis(vegdist(otu,method="bray") ~eGFR, data=env)
otu_env.adonis
```
Df，表示自由度；SumsOfSqs，总方差，又称离差平方和；Mean Sqs，表示均方（差），即SumsOfSqs/Df；F.Model，表示F检验值；R2，表示不同分组对样本差异的解释度，即分组方差与总方差的比值，R2越大表示分组对差异的解释度越高；Pr，表示P值，小于0.05说明本次检验的可信度高。
## 5 环境因子排序分析 RDA/CCA
CCA/RDA分析主要用来反映物种或功能与环境因子之间的关系。
冗余分析(redundancy analysis, RDA)或者典范对应分析(canonical correspondence analysis, CCA)是基于对应分析发展而来的一种排序方法，将对应分析与多元回归分析相结合，每一步计算均与环境因子进行回归，又称多元直接梯度分析。
### 5.1环境因子筛选
有时解释变量（环境因子）太多，我们需要想办法减少解释变量的数量。
减少解释变量的数量原因：
- 第一是寻求简约的模型，利于我们对模型的解读；
- 第二是当解释变量过多时会导致模型混乱，例如有些解释变量之间可能存在较强的线性相关，即共线性问题，可能会造成回归系数不稳定。
#### 方法1：bioenv筛选最佳因子组合
挑选合适的环境因子进行后续分析，当环境因子的个数较多时，可以使用bioenv 分析，挑选出相关性最高的环境因子的组合；
用于后续的分解分析
```{r}
res <- bioenv(otu,env,method="pearson")
res
summary(res)
res$whichbest
```
#### 方法2：方差膨胀因子分析:
VIF（variance inflation factor）方差膨胀因子分析，可分析环境因子之间是否存在相关性
删除具有共线性的因子：
计算每一个要进行筛选的环境因子对应的 VIF 值，通常的规则是 VIF 值大于20就被认定为该环境因子与其他环境因子有较强的相关性，应选择删掉其中的一个环境因子。进行多次筛选，直到选出所有的环境因子对应的 VIF 值全部小于20为止，更严格可筛选到10以下。
筛选原理：https://www.omicsclass.com/article/1359
```{r}
#所有的环境因子计算vif值
otu.tab.1<- rda(otu ~ ., env)
vif=vif.cca(otu.tab.1)
vif
#去掉vif最大的一个环境因子再次计算vif值，eGFR，直到所有的vif值小于20为止
otu.tab.1<- rda(otu ~ BMI+hypertension+WBC+RBC+Hb+PLT+X24hUTP+P+BUN+UA+Scr+T_chol+TG, env)
vif=vif.cca(otu.tab.1)
vif
#文章中的环境因子组合结果
otu.tab.1<- rda(otu ~Scr+BUN+ALB+P+Hb+eGFR, env)
vif=vif.cca(otu.tab.1)
vif
```
### 5.2 RDA/CCA分析
CCA的全称是典范对应分析 (Canonical correspondence analysis)，RDA的全称是冗余分析 (Redundancy analysis)。使用物种和环境因子进行的排序分析成为直接排序。
- CCA和RDA都属于排序分析，理论上其实和降维分析有点类似，都是尽量的将变量进行复杂的整合，使得排在前几位的变量组合能够尽可能的解释更多的信息。
- CCA和RDA都属于直接排序分析，区别是使用的排序方法不一样，CCA使用的是单峰模型，而RDA使用的是线性模型。
- CCA对应的间接排序方法是对应分析 (Correspondence analysis, CA)，RDA对应的间接排序方式是主成分分析 (Principal components analysis, PCA)。
- 在间接排序中的结果中，坐标轴是分析变量也就是物种的复杂函数组合，而在直接排序中结果中，坐标轴是环境因子的复杂函数组合。
我们在阅读文献的时候，会发现有的文章使用的是CCA，而有的文章使用的是RDA，那么对于初学者来说就会产生一种困惑，到底我是使用CCA还是RDA呢？
DCA分析可判别使用哪个方法更好: RDA还是CCA：
- 根据看分析结果中Axis Lengths的第一轴的大小
- 如果大于4.0,就应选CCA（基于单峰模型，典范对应分析）
- 如果在3.0-4.0之间，选RDA和CCA均可
- 如果小于3.0, RDA的结果会更合理（基于线性模型，冗余分析）
```{r}
#DCA 分析判断
dca = decorana(otu)
dca
```
这里可以看到：2.9257 应选择RDA分析
#### 5.2.1 RDA分析
```{r}
rda_result= rda(otu~Scr+BUN+ALB+P+Hb+eGFR,data=env,scale=TRUE)
rda_result
```
这里有6个环境因子，限制排序的轴为RDA开头，其他非限制轴为PC开头；
- 限制排序总比例解释比例为：Proportion Explained :0.05051
- 第一轴解释了：Proportion Explained 2.3010/163
- 第二轴：解释：1.7111/163
也可以总结一下 不用自己计算了
```{r}
summary(rda_result)
```
#### 5.2.2 CCA分析
这里选择文章中的环境因子，进行分析
```{r}
cca_result= cca(otu~Scr+BUN+ALB+P+Hb+eGFR,data=env,scale=TRUE)
cca_result
```
这里有6个环境因子，限制排序的轴为CCA开头，其他非限制轴为CA开头；
- 总比例解释比例为：Proportion Explained :
- 第一轴解释了：Proportion Explained
- 第二轴：解释：
```{r}
summary(cca_result)
```
### 5.2.3 环境因子解释是否显著
#### 以CCA分析为例
```{r}
#环境因子对微生物群落的解释是否显著？
#总体检验,所有约束轴的置换检验，以 999 次为例
cca_test <- anova (cca_result, permutations = 999)
cca_test
#各轴解释度检验显著性，环境因子对应的轴解释
cca_test_axis <- anova(cca_result, by = 'axis', permutations = 999)
cca_test_axis
#矫正p值，因为检验的环境因子越多，会得到越多显著的结果
cca_test_axis$`Pr(>F)` <- p.adjust(cca_test_axis$`Pr(>F)`, method = 'bonferroni')
cca_test_axis
#各环境因子显著性与解释变差置换检验 这里的anova.cca 也可以写成anova
cca_test_env <- anova.cca(cca_result, by = 'margin', permutations = 999)
cca_test_env
#p 值校正（Bonferroni 为例）
cca_test_env$`Pr(>F)` <- p.adjust(cca_test_env$`Pr(>F)`, method = 'bonferroni')
```
#### 以RDA分析为例
```{r}
#环境因子对微生物群落的解释是否显著？
#总体检验,所有约束轴的置换检验，以 999 次为例
rda_test <- anova(rda_result, permutations = 999)
rda_test
#各轴解释度检验显著性，环境因子对应的轴解释
rda_test_axis <- anova(rda_result, by = 'axis', permutations = 999)
rda_test_axis
#矫正p值，因为检验的环境因子越多，会得到越多显著的结果
rda_test_axis$`Pr(>F)` <- p.adjust(rda_test_axis$`Pr(>F)`, method = 'bonferroni')
rda_test_axis
#各环境因子显著性与解释变差置换检验 这里的anova.cca 也可以写成anova
rda_test_env <- anova.cca(rda_result, by = 'margin', permutations = 999)
rda_test_env
#p 值校正（Bonferroni 为例）
rda_test_env$`Pr(>F)` <- p.adjust(rda_test_env$`Pr(>F)`, method = 'bonferroni')
```
### 5.2.4 单独的某个环境因子解释多少？(选修)
解释方差的置换检验，以 eGFR 所能解释的全部方差为例；999 次置换，详情 ?anova.cca
例如，对eGFR所解释变差部分的置换检验如下所示。
```{r}
anova.cca(rda(otu, env['eGFR']), permutations = 999)
#若只考虑eGFR单独解释的方差部分，需将其它变量作为协变量如Hb
anova.cca(rda(otu, env['eGFR'], env[c("Hb")]), permutations = 999)
```
## 6 RDA CCA 分析结果排序图展示
I型标尺:如果分析兴趣在于解读样方之间的关系，设定scaling=1。
II型标尺:如果分析兴趣在于解读物种之间的关系，设定scaling=2。
display 参数解释：sp = "species", wa = "sites", bp = "biplot"
```{r}
#提取排序得分值
#scaling: "none", "sites", "species", or "symmetric", which correspond to the values 0, 1, 2, and 3 respectively
ccascore <- scores(cca_result,choices = 1:4,display =c("sp","wa","bp") ,scaling = 1)
ccascore
#输出结果保存
write.csv(ccascore$sites,file="cca.sample.csv") #样本排序坐标
write.csv(cca_result$CCA$biplot,file="cca.env.csv") #环境因子排序坐标
write.csv(ccascore$species,file="cca.species.csv") #物种排序坐标
##绘图开始准备
#提取环境因子坐标
CCAE <- as.data.frame(cca_result$CCA$biplot[,1:2])
#提取样方坐标
CCAS1 <- ccascore$sites[,1]
CCAS2 <- ccascore$sites[,2]
#创建绘图数据
plotdata <- data.frame(rownames(ccascore$sites), CCAS1, CCAS2, fastmap$group)
colnames(plotdata) <- c("sample","CCAS1","CCAS2","group")
head(plotdata)
#提取前两轴解释百分比 cca_result$CCA$tot.chi 表示Constrained，cca_result$CA$tot.chi Unconstrained序部分 rda也是一样
cca1 <- round(cca_result$CCA$eig[1]/cca_result$CCA$tot.chi*100,2)
cca2 <- round(cca_result$CCA$eig[2]/cca_result$CCA$tot.chi*100,2)
#模型总解释比例
cca <- round(sum(cca_result$CCA$eig)/cca_result$tot.chi*100,2)
##开始绘图
#颜色，形状准备
mycol=c(brewer.pal(9,"Set1"),brewer.pal(7,"Set2"),brewer.pal(12,"Set3"),brewer.pal(12,"Paired"),brewer.pal(8,"Dark2"),brewer.pal(7,"Accent"))
myshape <- c(15:20,1:14,21:25)
#CCA plot
plot_CCA<-ggplot(plotdata, aes(CCAS1, CCAS2)) +
geom_point(aes(color = group,shape=group),size = 4) +
scale_shape_manual(values = myshape)+
scale_color_manual(values = mycol)+
labs(title = paste("Var= ",cca,"%,"," Pvalue=",cca_test$`Pr(>F)`, sep = "") ) +
xlab(paste("CCA1 ( ",cca1,"%"," )", sep = "")) +
ylab(paste("CCA2 ( ",cca2,"%"," )", sep = "")) +
geom_segment(data = CCAE, aes(x = 0, y = 0, xend = CCAE[,1], yend = CCAE[,2]),
colour = "blue", size = 1.2,
arrow = arrow(angle = 30, length = unit(0.4, "cm"))) +
geom_label_repel(data = CCAE, fill = "white",segment.colour = "black", aes(x = CCAE[,1], y = CCAE[,2], label = rownames(CCAE))) +
geom_vline(aes(xintercept = 0), linetype = "dotted") +
geom_hline(aes(yintercept = 0), linetype = "dotted") +
theme(panel.background = element_rect(fill = "white", colour = "black"),
panel.grid = element_blank(),
axis.title = element_text(color = "black", size = 18),
axis.ticks.length = unit(0.4,"lines"),
axis.ticks = element_line(color = "black"),
axis.line = element_line(colour = "black"),
axis.title.x = element_text(colour = "black", size = 18),
axis.title.y = element_text(colour="black", size = 18),
axis.text = element_text(colour = "black", size = 18),
legend.title = element_blank(),
legend.text = element_text(size = 18), legend.key = element_blank(),
plot.title = element_text(size = 22, colour = "black",
face = "bold", hjust = 0.5))
#是否添加样本名称到排序图中，如果添加运行下面代码，
#plot_CCA=plot_CCA+geom_label_repel(aes(CCAS1, CCAS2, label = sample), fill = "white",
# color = "black", box.padding = unit(0.6,"lines"),
# segment.colour = "grey50",
# label.padding = unit(0.35,"lines"))
plot_CCA
```
CCA和RDA图查看说明：
- CCA和RDA的结果图中使用点代表不同的样本，从原点发出的箭头代表不同的环境因子。
- 箭头的长度代表该环境因子对群落变化影响的强度，箭头的长度越长，表示环境因子的影响越大。
- 箭头与坐标轴的夹角代表该环境因子与坐标轴的相关性，夹角越小，代表相关性越高。
- 样本点到环境因子箭头极其延长线的垂直距离表示环境因子对样本的影响强度，样本点与箭头距离越近，该环境因子对样本的作用越强。
- 样本位于箭头同方向，表示环境因子与样本物种群落的变化正相关，样本位于箭头的反方向，表示环境因子与样本物种群落的变化负相关。
- 图像中坐标轴标签中的数值，代表了坐标轴所代表的环境因子组合对物种群落变化的解释比例。
## 7 db-RDA 分析
采用的距离为欧式距离以外的距离作为输入进行rda分析,例如bray curtis RDAcca分析就变成了 db-RDA分析;普通的RDA 默认为欧式距离；
CAP（Canonical analysis of principal coordinates）分析就是db-RDA。
### 7.1 分步完成db-RDA 分析
```{r}
#计算bray curtis 距离矩阵
dis_bray <- vegdist(otu, method = 'bray')
#或者直接使用已准备好的距离矩阵
#dis_bray <- as.dist(read.delim('data/bray_distance.txt', row.names = 1, sep = '\t', check.names = FALSE))
#db-rda分析可以capscale 完成
db_rda <- capscale(dis_bray~Scr + BUN + ALB + P + Hb + eGFR, env, add = TRUE)
db_rda
```
### 7.2 一步到位直接出结果
capscale()是vegan中打包好的db-RDA执行函数，可以直接一步出结果。
```{r}
db_rda <- capscale(otu~Scr + BUN + ALB + P + Hb + eGFR, env, distance = 'bray', add = TRUE,sqrt.dist = TRUE)
db_rda
#检验一下分析结果是否有显著性
db_rda_test = anova(db_rda, permutations = 1000, parallel = 4)
db_rda_test
```
### 7.3 结果可视化
```{r}
#前4轴
db_rdascore <- scores(db_rda,choices = 1:4,display =c("sp","wa","bp") ,scaling = 1)
#输出结果
write.csv(db_rdascore$sites,file="db_rda.sample.csv")
write.csv(db_rda$CCA$biplot,file="db_rda.env.csv")
write.csv(db_rdascore$species,file="db_rda.species.csv")
#提取环境因子坐标
CAPE <- as.data.frame(db_rda$CCA$biplot[,1:2])
#提取样方坐标
CAPS1 <- db_rdascore$sites[,1]
CAPS2 <- db_rdascore$sites[,2]
#创建绘图数据
plotdata <- data.frame(rownames(ccascore$sites), CAPS1, CAPS2, fastmap$group)
colnames(plotdata) <- c("sample","CAPS1","CAPS2","group")
head(plotdata)
#提取前两轴解释百分比
cap1 <- round(db_rda$CCA$eig[1]/db_rda$CCA$tot.chi*100,2)
cap2 <- round(db_rda$CCA$eig[2]/db_rda$CCA$tot.chi*100,2)
#总解释比例
cap <- round(sum(db_rda$CCA$eig)/db_rda$tot.chi*100,2)
#开始绘图
mycol=c(brewer.pal(9,"Set1"),brewer.pal(7,"Set2"),brewer.pal(12,"Set3"),brewer.pal(12,"Paired"),brewer.pal(8,"Dark2"),brewer.pal(7,"Accent"))
myshape <- c(15:20,1:14,21:25)
#CAP plot
p<-ggplot(plotdata, aes(CAPS1, CAPS2)) +
geom_point(aes(color = group,shape=group),size = 4) +
scale_shape_manual(values = myshape)+
scale_color_manual(values = mycol)+
labs(title = paste("Var= ",cap,"%,"," Pvalue=",db_rda_test$`Pr(>F)`, sep = "") ) +
xlab(paste("CAP1 ( ",cap1,"%"," )", sep = "")) +
ylab(paste("CAP2 ( ",cap2,"%"," )", sep = "")) +
geom_segment(data = CAPE, aes(x = 0, y = 0, xend = CAPE[,1], yend = CAPE[,2]),
colour = "blue", size = 1.2,
arrow = arrow(angle = 30, length = unit(0.4, "cm"))) +
geom_label_repel(data = CAPE, fill = "white",segment.colour = "black", aes(x = CAPE[,1], y = CAPE[,2], label = rownames(CAPE))) +
geom_vline(aes(xintercept = 0), linetype = "dotted") +
geom_hline(aes(yintercept = 0), linetype = "dotted") +
theme(panel.background = element_rect(fill = "white", colour = "black"),
panel.grid = element_blank(),
axis.title = element_text(color = "black", size = 18),
axis.ticks.length = unit(0.4,"lines"),
axis.ticks = element_line(color = "black"),
axis.line = element_line(colour = "black"),
axis.title.x = element_text(colour = "black", size = 18),
axis.title.y = element_text(colour="black", size = 18),
axis.text = element_text(colour = "black", size = 18),
legend.title = element_blank(),
legend.text = element_text(size = 18), legend.key = element_blank(),
plot.title = element_text(size = 22, colour = "black",
face = "bold", hjust = 0.5))
#是否添加样本名称到排序图中，如果添加运行下面代码，
#p=p+geom_label_repel(aes(CAPS1, CAPS2, label = sample), fill = "white",
# color = "black", box.padding = unit(0.6,"lines"),
# segment.colour = "grey50",
# label.padding = unit(0.35,"lines"))
p
pdf(file="db-rda.pdf", height=7, width=7)
print(p)
dev.off()
png(filename="db-rda.png", height=7*300, width=7*300, res=300, units="px")
print(p)
dev.off()
```
- 发表于 2021-02-26 15:28
- 阅读 ( 7213 )
- 分类：宏基因组
你可能感兴趣的文章
- MENA如何结合OTU共现性网络与环境因子 11528 浏览
