常用的排序方法如下：
間接排序法包括：
* PCA（principal components analysis涤浇，主成分分析）
* CA（correspondence analysis厕倍，對應(yīng)分析）
* DCA（Detrended correspondenceanalysis, 去趨勢對應(yīng)分析)
* PCoA（principal coordinate analysis，主坐標(biāo)分析）
* NMDS（non-metric multi-dimensional scaling，非度量多維尺度分析）悍赢；
直接排序法包括：
* RDA（Redundancy analysis，冗余分析）
* CCA（canonical correspondence analysis认境，典范對應(yīng)分析）
* dbRDA（distance based redundancy analysis硕淑，基于距離的冗余分析）
* CAP（canonical analysis of principal coordinates课竣，主要坐標(biāo)的典型分析）
其中PCA和RDA是基于線性模型（linear model）的嘉赎，
而CA、DCA于樟、CCA公条、DCCA是基于單峰（unimodal）模型
選擇單峰模型還是線性模型？
用DCA（vegan::decorana()）先對數(shù)據(jù)（site-species）進(jìn)行分析迂曲；
查看結(jié)果中的“Axis lengths”的第一軸DCA1的值,
根據(jù)該值判斷該采用線性模型還是單峰模型：
如果大于4.0赃份，就應(yīng)該選單峰模型；
如果3.0-4.0之間奢米，選線性模型或者單峰模型均可；
如果小于3.0, 線性模型的結(jié)果要好于單峰模型
rm(list = ls())
pacman::p_load(tidyverse,vegan,ggrepel)
species <- read.delim("phylum.top10.xls",header = T,
row.names = 1,check.names = F,sep="\t") %>%
t() %>% decostand(method = "hellinger")
envi <- read.delim("env1.log10.txt",
header = T,row.names = 1,check.names = F,sep="\t")
grp <- read.delim("group.xls",header = T,check.names = F,sep="\t")
decorana(species)
#根據(jù)看分析結(jié)果中Axis Lengths的第一軸的大小
#如果大于4.0,就應(yīng)選CCA（基于單峰模型纠永，典范對應(yīng)分析）
#如果在3.0-4.0之間鬓长，選RDA和CCA均可
#如果小于3.0, RDA的結(jié)果會(huì)更合理（基于線性模型，冗余分析）
p <- rda(species~.,envi) %>% summary()
sp <- as.data.frame(p$species[,1:2])*2
st <- as.data.frame(p$sites[,1:2])
yz <- as.data.frame(p$biplot[,1:2])
ggplot() + geom_point(data=st,aes(RDA1,RDA2,
shape=grp$group,fill=grp$group),size=3)+
scale_shape_manual(values =c(18,19,21))+
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 = "#E69F00")+
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=45,length = unit(0.35,"cm"),
type = "closed"),linetype=1,
size=0.6,colour = "steelblue")+
geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+
labs(x=paste("RDA1(", format(100*p$cont[[1]][2,1],
digits=4),"%)",sep=""),
y=paste("RDA2 (", format(100*p$cont[[1]][2,2],
digits=4),"%)", sep=""))+
geom_hline(yintercept=0,linetype=3,size=1) +
geom_vline(xintercept=0,linetype=3,size=1)+
guides(shape=guide_legend(title=NULL,color="black"),
fill=guide_legend(title=NULL))+
theme_bw()+theme(panel.grid=element_blank())
rda.jpeg
rm(list=ls())
pacman::p_load(tidyverse,vegan,ggrepel)
sample <- read.table("otu_table.txt", header = T,
row.names=1,sep="\t") %>% t()
env <- read.table("environment.txt", header=TRUE, row.names=1,sep="\t")
group <- read.table("groups.txt", header=T) %>% as.list()
decorana(sample)
cca <- cca(sample, env, scale = TRUE)
p <- scores(cca)
CCAE <- as.data.frame(cca$CCA$biplot[,1:2])
plotdata <- data.frame(rownames(p$sites), p$sites[,1],p$sites[,2], group$group)
colnames(plotdata) <- c("sample","CCAS1","CCAS2","group")
#----------------------------------------------------------------------------------
ggplot(plotdata, aes(CCAS1, CCAS2)) +
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")) +
geom_point(aes(fill = group,color=group),size = 3) +
labs(title = "CCA Plot",
x = round(cca$CCA$eig[1]/sum(cca$CCA$eig)*100,2) %>%
paste0("CCA1 ( ",.,"%"," )"),
y = round(cca$CCA$eig[2]/sum(cca$CCA$eig)*100,2) %>%
paste0("CCA2 ( ",.,"%"," )")) +
geom_segment(data = CCAE, aes(x = 0, y = 0, xend = CCAE[,1], yend = CCAE[,2]),
colour = "black", size = 1.2,
arrow = arrow(angle = 30, length = unit(0.4, "cm"))) +
geom_label_repel(data = CCAE, fill = "grey90",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 = 12),
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 = 12),
axis.title.y = element_text(colour="black", size = 12),
axis.text = element_text(colour = "black", size = 12),
legend.title = element_blank(),
legend.text = element_text(size = 12), legend.key = element_blank(),
plot.title = element_text(size = 14, colour = "black",
face = "plain", hjust = 0.5))
cca.jpeg
