####Abundance and species diversity indices with depth------------

spe <- read.csv("spe.csv",row.names = 1,header = T) 
spe <- spe[,1:88]
env <- read.csv("env.csv",row.names = 1,header = T)

nonzero.spe <- which(rowSums(spe)>0)
spe <- spe[nonzero.spe,]
env <- env[nonzero.spe,]

ord.depth <- order(env$depth)

spe <- spe[ord.depth,]
env <- env[ord.depth,]


#Species richness

spe.count <- apply(spe>0, 1, sum)
spe.rich.depth <- data.frame(spe.count,env$depth)
spe.rich.depth$group <- cut(spe.rich.depth$env.depth,breaks = seq(500,4100,100),
                            labels = seq(500,4000,by=100),right = FALSE)
names(spe.rich.depth) <- c("spe.rich","depth","group")
spe.rich.depth <- as.data.frame(spe.rich.depth)

lable <- c("","","","","","","","","","","","","","","","","","","","","","","","","","",
           "","","","","","","","","","")
lable2 <- c("500","","","","","1000","","","","","1500","","","","","2000","","","","","2500","","","","","3000",
            "","","","","3500","","","","","4000")


p1 <- ggplot(spe.rich.depth,  aes(x=group, y = spe.rich))+
  geom_boxplot(outlier.size = 1,color="steelblue" ,fill = 'white')+
  scale_x_discrete(labels =lable)+
  theme(
    panel.grid = element_blank(), 
    panel.background = element_rect(fill = 'transparent', color = 'black'),
    axis.text.y = element_text(size=15,color="black"),
    axis.text.x = element_text(size=14,color="black"),
    axis.title = element_text(size=17)
  ) +
  labs(x = '', y = 'Richness')

p1


#Species abundance

spe.sum <- apply(spe, 1, sum)

spe.abun.depth <- data.frame(spe.sum,env$depth)
spe.abun.depth$group <- cut(spe.abun.depth$env.depth,breaks = seq(500,4100,100),
                            labels = seq(500,4000,by=100),right = FALSE)
names(spe.abun.depth) <- c("spe.abun","depth","group")


p2 <- ggplot(spe.abun.depth,  aes(x=group, y = spe.abun))+
  geom_boxplot(outlier.size = 1,color="steelblue" ,fill = 'white')+
  scale_x_discrete(labels =lable)+
  theme(
    panel.grid = element_blank(), 
    panel.background = element_rect(fill = 'transparent', color = 'black'),
    axis.text.y = element_text(size=15,color="black"),
    axis.text.x = element_text(size=14,color="black"),
    axis.title = element_text(size=17)
  ) +
  labs(x = '', y = 'Abundance/ind.')


p2

#shannon diversity

library(vegan)
library(picante)

#Shannon index
shannon_index <- diversity(spe, index = 'shannon', base = exp(1))

#Shannon diversity
shannon_diversity <- exp(1)^shannon_index

spe.shan.depth <- data.frame(shannon_diversity,env$depth)
spe.shan.depth$group <- cut(spe.shan.depth$env.depth,breaks = seq(500,4100,100),
                            labels = seq(500,4000,by=100),right = FALSE)
names(spe.shan.depth) <- c("spe.shan","depth","group")


p3 <- ggplot(spe.shan.depth,  aes(x=group, y = spe.shan))+
  geom_boxplot(outlier.size = 1, color = 'steelblue',fill="white")+
  scale_x_discrete(labels =lable2)+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text.y = element_text(size=15,color="black"),
        axis.text.x = element_text(size=14,color="black"),
        axis.title = element_text(size=17)
  ) +
  labs(x = 'Depth', y = 'Shannon’s index')


p3


######Pielou evenness index

richness <- rowSums(spe > 0)
pielou <- shannon_index / log(richness, exp(1))

spe.pielou.depth <- data.frame(pielou,env$depth)
spe.pielou.depth$group <- cut(spe.pielou.depth$env.depth,breaks = seq(500,4100,100),
                              labels = seq(500,4000,by=100),right = FALSE)
names(spe.pielou.depth) <- c("spe.pielou","depth","group")


p4 <- ggplot(spe.pielou.depth,  aes(x=group, y = spe.pielou))+
  geom_boxplot(outlier.size = 1, color = 'steelblue',fill="white")+
  scale_x_discrete(labels =lable2)+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text.y = element_text(size=15,color="black"),
        axis.text.x = element_text(size=14,color="black"),
        axis.title = element_text(size=17)
  ) +
  labs(x = 'Depth', y = 'Pielou’s evenness')

p4

library(cowplot)
plot6 <- plot_grid(p2,p1,p3,p4)
plot6

#存为pdf
pdf("Abundance and species diversity indices with depth.pdf")

plot6

garbage <- dev.off()

####Species cumulative curves-----

library(vegan)

spe2 <- read.csv("spe.csv",header=T,row.names = 1)

dev.new(title = "accumulation curve", noRStudioGD = TRUE)
sp1 <- specaccum(spe2,method = "random")
plot(sp1,ci.type = "poly",col = "blue",lwd=2,ci.lty = 0,ci.col = "lightblue",
     xlab = "Number of samples",ylab = "Number of morphospecies")

#存为pdf
pdf("Species cumulative curves.pdf")

plot(sp1,ci.type = "poly",col = "blue",lwd=2,ci.lty = 0,ci.col = "lightblue",
     cex.axis=1.3,cex.lab=1.4,
     xlab = "Number of samples",ylab = "Number of morphospecies")

garbage <- dev.off()


#####environmental variables

library(ggplot2)
library(cowplot)


env <- read.csv("env.csv",row.names = 1,header = T)

label.pairs2 <- c("","","","","","","")
label.pairs1 <- c("C","A","B","F","G","E","D")
label.pairs <- c("JL177","JL169","JL175","JL201","JL202","JL199","JL178")

p1 <- ggplot(env,  aes(x=reorder(group,depth), y = depth))+
  geom_boxplot(outlier.size = 1, fill = "steelblue")+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = '', y = 'Depth(m)')+
  scale_x_discrete(labels =label.pairs2 )
p2 <- ggplot(env,  aes(x=reorder(group,depth), y = slope))+
  geom_boxplot(outlier.size = 1, fill = 'steelblue')+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = '', y = 'Slope(°)')+
  scale_x_discrete(labels =label.pairs2 )
p3 <- ggplot(env,  aes(x=reorder(group,depth), y = curvature))+
  geom_boxplot(outlier.size = 1, fill = 'steelblue')+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = '', y = 'Curvature')+
  scale_x_discrete(labels =label.pairs2 )
p4 <- ggplot(env,  aes(x=reorder(group,depth), y = roughness))+
  geom_boxplot(outlier.size = 1, fill = 'steelblue')+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = '', y = 'Roughness')+
  scale_x_discrete(labels =label.pairs2 )
p5 <- ggplot(env,  aes(x=reorder(group,depth), y = fineBPI))+
  geom_boxplot(outlier.size = 1, fill = 'steelblue')+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = '', y = 'FineBPI')+
  scale_x_discrete(labels =label.pairs2 )
p6 <- ggplot(env,  aes(x=reorder(group,depth), y = broadBPI))+
  geom_boxplot(outlier.size = 1, fill = 'steelblue')+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = '', y = 'BroadBPI')+
  scale_x_discrete(labels =label.pairs2 )
p7 <- ggplot(env,  aes(x=reorder(group,depth), y = backscatter))+
  geom_boxplot(outlier.size = 1, fill = 'steelblue')+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = 'Dives', y = 'Backscatter(db)')+
  scale_x_discrete(labels =label.pairs1 )
p8 <- ggplot(env,  aes(x=reorder(group,depth), y = meanVelocity))+
  geom_boxplot(outlier.size = 1, fill = 'steelblue')+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = 'Dives', y = 'MeanVelocity(m/s)')+
  scale_x_discrete(labels =label.pairs1 )
p9 <- ggplot(env,  aes(x=reorder(group,depth), y = substrate))+
  geom_boxplot(outlier.size = 1, fill = 'steelblue')+
  theme(axis.text = element_text(size=10,angle = 0,color="black"),
        axis.title = element_text(size=12),
        panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black')) +
  labs(x = 'Dives', y = 'Substrate')+
  scale_x_discrete(labels =label.pairs1 )

plot6 <- plot_grid(p1,p2,p3,p4,p5,p6,p8,p7,p9,nrow=3,align = 'v',labels = "AUTO")

####存为pdf
pdf("environmental variables.pdf")

plot6

garbage <- dev.off()



########Cluster=====================================


library(ade4)
library(adespatial)
library(vegan)
library(gclus)
library(cluster)
library(pvclust)
library(RColorBrewer)
library(labdsv)
library(rioja)
library(indicspecies)
library(mvpart)
library(MVPARTwrap)
library(dendextend)
library(vegclust)
library(colorspace)
library(agricolae)
library(picante)


source("drawmap.R")
source("drawmap3.R")
source("hcoplot.R")
source("test.a.R")
source("coldiss.R")
source("bartlett.perm.R")
source("boxplerk.R")
source("boxplert.R")


grpdist <- function(X)
{
  require(cluster)
  gr <- as.data.frame(as.factor(X))
  distgr <- daisy(gr, "gower")
  distgr
}


#读入文件

spe <- read.csv("spe.csv",row.names = 1,header = T) 
env <- read.csv("env.csv",row.names = 1,header = T)


spe <- subset(spe, select = colSums(spe) != 0)
spe <- subset(spe, select = colSums(spe) != 1)
spe <- spe[,1:72]

nonzero.spe <- which(rowSums(spe)>0)

spe <- spe[nonzero.spe,]
env <- env[nonzero.spe,]

ord.depth <- order(env$depth)

spe <- spe[ord.depth,]
env <- env[ord.depth,]

colnames(spe)

# Fuzzy c-means clustering 

k <- 6		 
spe.fuz <- fanny(spe.ch, k = k, memb.exp = 1.2)
summary(spe.fuz)


spe.fuz$membership

spe.fuz$clustering
spefuz.g <- spe.fuz$clustering


spe.norm <- decostand(spe, "normalize")
spe.ch <- vegdist(spe.norm, "euc")

spe.ch <- dist.ldc(spe, "chord")


dev.new(title = "Fuzzy clustering of fish data - Silhouette plot",
        noRStudioGD = TRUE)
plot(
  silhouette(spe.fuz),
  main = "Silhouette plot - Fuzzy clustering",
  cex.names = 0.8,
  col = spe.fuz$silinfo$widths + 1
)


dc.pcoa <- cmdscale(spe.ch)
dc.scores <- scores(dc.pcoa, choices = c(1, 2))


# Noise clustering (vegclust) =====================================


k <- 6
## Create noise clustering with 6 clusters
spe.nc <- vegclust(
  spe.norm,
  mobileCenters = k,
  m = 1.5,
  dnoise = 0.75,
  method = "NC",
  nstart = 30
)
spe.nc

# Medoids of species
(medoids <- spe.nc$mobileCenters)

# Fuzzy membership matrix
spe.nc$memb

#  Cardinality of fuzzy clusters (i.e., the number of objects belonging to each cluster)
spe.nc$size

# Obtain hard membership vector, with 'N' for objects that are unclassified

spefuz.g <- defuzzify(spe.nc$memb)$cluster
clNum <- as.numeric(as.factor(spefuz.g))

ccccccluster <- as.data.frame(spefuz.g)

cluster <- as.data.frame(spefuz.g)

cluster <- read.csv("cluster.csv",header = T,row.names = 1)

#######plot 

a <- as.data.frame(dc.scores)
a <- data.frame(a,as.factor(cluster$cluster))
names(a) <- c("Dim1","Dim2","cluster")

dev.new(title = "Defuzzified site plot",
        noRStudioGD = TRUE) 


#存为pdf
pdf("Noise clustering .pdf")

p <- ggplot(data = a, aes(x = Dim1, y = Dim2)) +
  theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.key = element_rect(fill = 'transparent'),
        legend.position = c(0.078,0.85),
        legend.key.size = unit(15,"pt"),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13), 
        axis.text = element_text(size=14,color="black"),
        axis.title = element_text(size=15),
  ) +
  labs(x = 'MDS1', y = 'MDS2')

p + geom_point(aes(color=cluster,shape=cluster),size=3) +
  scale_color_manual(values = c(2:7,"grey"))+
  scale_shape_manual(values = c(1:7))


garbage <- dev.off()

#####Medoids of species


mediods.t <- t(medoids)

write.csv(mediods.t,file='output/mediods.t.6.csv')

##### environmental variables for each cluster.


env5.6 <- data.frame(env,cluster$cluster)
names(env5.6) <- c("depth","slope","aspect","curvature","roughness","broadBPI","fineBPI","backscatter","meanVelocity","substrate","cluster")


write.csv(env5.6,file='output/env5.6.csv')

###plot
library(ggplot2)
library(cowplot)

env5.6 <-  read.csv("env5.6.csv",header = T,row.names = 1)
env5.6$cluster <- cluster

label.pairs2 <- c(" "," "," "," "," "," "," ")

a <- ggplot(env5.6,aes(x=cluster,y=depth,fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3) +guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = '', y = 'Depth(m)')+
  scale_x_discrete(labels =label.pairs2 )
b <- ggplot(env5.6,aes(x=cluster,y=slope,fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3)+guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = '', y = 'Slope(°)')+
  scale_x_discrete(labels =label.pairs2 )
c <- ggplot(env5.6,aes(x=cluster,y=log(curvature),fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3)+guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = '', y = 'Curvature')+
  scale_x_discrete(labels =label.pairs2 )
d <- ggplot(env5.6,aes(x=cluster,y=roughness,fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3) +guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = '', y = 'Roughness')+
  scale_x_discrete(labels =label.pairs2 )
e<- ggplot(env5.6,aes(x=cluster,y=fineBPI,fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3) +guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = '', y = 'FineBPI')+
  scale_x_discrete(labels =label.pairs2 )
f<- ggplot(env5.6,aes(x=cluster,y=broadBPI,fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3) +guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = '', y = 'BroadBPI')+
  scale_x_discrete(labels =label.pairs2 )
g<- ggplot(env5.6,aes(x=cluster,y=backscatter,fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3) +guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = 'Cluster', y = 'Backscatter(db)')
h<- ggplot(env5.6,aes(x=cluster,y=meanVelocity,fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3) +guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = 'Cluster', y = 'MeanVelocity(m/s)')
i<- ggplot(env5.6,aes(x=cluster,y=substrate,fill=cluster))+
  stat_boxplot(geom = "errorbar", width=0.3) +guides(fill=F)+geom_boxplot()+
  scale_fill_manual(values = c(2:7,"grey"))+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
  ) + labs(x = 'Cluster', y = 'Substrate')


plot6 <- plot_grid(a,b,c,d,e,f,g,i,h)

dev.new(title = "environmental variables for each cluster category",
        noRStudioGD = TRUE)
plot6


pdf("environmental variables for each cluster.pdf")

plot6

garbage <- dev.off()

####LCBD和SCBD=====================================

library(ade4)
library(adegraphics)
library(adespatial)
library(vegan)
library(vegetarian)
library(ggplot2)
library(FD)
library(taxize)

source("panelutils.R")
source("Rao.R")


#spe <- read.csv("SpeG.csv",row.names = 1,header = T)

# Computation using beta.div {adespatial} on Hellinger-transformed species data
spe.beta <- beta.div(spe, method = "hellinger", nperm = 9999)
summary(spe.beta)
spe.beta$beta  # SSTotal and BDTotal


# Computation using beta.div {adespatial} on ochiai-transformed species data
spe.beta <- beta.div(spe, method = "ochiai", nperm = 9999)
summary(spe.beta)
spe.beta$beta  # SSTotal和BDTotal

####SCBD

# Which species have a SCBD larger than the mean SCBD
spe.beta$SCBD[spe.beta$SCBD >= mean(spe.beta$SCBD)]

#the 20 species with the highest contribution to beta diversity

f5 <- data.frame(spe.beta$SCBD[spe.beta$SCBD >= mean(spe.beta$SCBD)])
f5 <- data.frame(row.names(f5),f5$spe.beta.SCBD.spe.beta.SCBD....mean.spe.beta.SCBD..)
abun <- c(28,602,2374,1892,962,102,21,219,21,60,30,81,6,236,51,21,79,92,176,227)
occur <- c(18,45,38,30,21,55,15,25,14,23,21,40,6,39,15,18,44,17,62,75)
f5 <- data.frame(f5, abun/10222,occur/130)
names(f5) <- c("species","SCBD","prop","freq")

f5 <- data.frame(f5[order(f5$SCBD,decreasing = T),])



####存为pdf
pdf("SCBD.pdf")

ggplot(data=f5, aes(x=reorder(species, -SCBD),y=SCBD))+
  geom_bar(stat="identity",width=0.8,fill="steelblue")+
  theme_classic()+theme(axis.text = element_text(size=13,color="black"),
                        axis.title = element_text(size=14),
                        axis.text.x = element_text(angle=45, hjust=1, vjust=1),legend.position = "top")+
  labs(x="Species",y="Value of contribution")+scale_y_continuous(sec.axis = sec_axis(~.*1000,name = 'Persent(100%)'))+
  scale_x_discrete(labels =names)+
  geom_line(data = f5,mapping = aes(x=reorder(species,-SCBD),y=freq*0.1,group=1),size=1,color="#1b7c3d")+
  geom_point(data = f5,mapping = aes(x=reorder(species,-SCBD),y=freq*0.1,group=1),size=3,stroke=1,shape=0,color="#1b7c3d")+
  geom_line(data = f5,mapping = aes(x=reorder(species,-SCBD),y=prop*0.1,group=2),size=1,color="#CB0000")+
  geom_point(data = f5,mapping = aes(x=reorder(species,-SCBD),y=prop*0.1,group=2),size=3,stroke=1,shape=2,color="#CB0000")+
  geom_rect(aes(xmin=10.5,xmax=12,ymin=0.07,ymax=0.073),fill='steelblue',color='steelblue')+
  annotate(geom='text',x=16.3,y=0.072,label='Contribution to beta diversity index',size=4)+
  annotate("segment",x = 10.5,xend = 12,y = 0.068,yend = 0.068,size=1,color="#CB0000")+
  geom_point(aes(x=11.25, y=0.068), size=3,stroke=1,shape=2,colour="#CB0000")+
  annotate(geom='text',x=15.2,y=0.068,label="Proportion of abundance",size=4)+
  annotate("segment",x = 10.5,xend = 12,y = 0.065,yend = 0.065,size=1,color="#1b7c3d")+
  geom_point(aes(x=11.25, y=0.065), size=3,stroke=1,shape=0,colour="#1b7c3d")+
  annotate(geom='text',x=15.2,y=0.065,label="Frequency of occurrence",size=4)

garbage <- dev.off()

####LCBD

spefuz.g <- read.csv("cluster.csv",header = T,row.names = 1)
f6 <- data.frame(spe.beta$LCBD)
f6 <- data.frame(row.names(f6),f6$spe.beta.LCBD)
c <- data.frame(env$depth,env$substrate,spefuz.g)
f6 <- data.frame(f6,c)
names(f6) <- c("samples","LCBD","depth","substrate","cluster")

substrate2 <- rep(".very_soft", nrow(env))
substrate2[env$substrate == 2] <- ".soft"
substrate2[env$substrate == 3] <- ".mixed"
substrate2[env$substrate == 4] <- ".hard"
substrate2[env$substrate == 5] <- ".very_hard"
substrate2 <- factor(substrate2, 
                     levels = c(".very_soft", ".soft", ".mixed", ".hard", ".very_hard"))
table(substrate2)
# Create an env3 data frame with slope as a qualitative variable
f7 <- f6
f7$substrate <- substrate2


###Contribution of samples to beta diversity for each taxonomic cluster


p <- ggplot(data = f7, aes(x = depth, y = LCBD)) +
  theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.key = element_rect(fill = 'transparent'),
        axis.text = element_text(size=14,color="black"),
        axis.title = element_text(size=15),
        legend.position = c(0.17,0.9),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13), 
        legend.direction = "horizontal"
  ) +
  guides(
    color=guide_legend(title = "Cluster",
                       ncol = 2,
                       byrow = F,
                       reverse = F
    ))+
  labs(x = 'Depth', y = 'Value of contribution')

p + geom_point(aes(color=cluster),size=3) +
  scale_color_manual(values = c(2:7,"grey"))


pdf("LCBD.pdf")

p + geom_point(aes(color=cluster),size=3) +
  scale_color_manual(values = c(2:7,"grey"))

garbage <- dev.off()

####NMDS------

#rm(list = ls())
library(vegan)
library(metR)
library(ggrepel)
library(ggplot2)


spdt <- data.frame(env$depth,env$group,spe)
colnames(spdt)[colnames(spdt) == ' group '] <- 'Dives'
head(spdt)
#Taxa_classfi <- read.csv("Taxa_classfi.csv",row.names = 1)
#head(Taxa_classfi)
spmds <- metaMDS(spdt[,3:ncol(spdt)],distance =  "bray" )
spst<- as.data.frame(scores(spmds,"sites"))
spst <- data.frame(spst,tr=spdt$Depth,Time=spdt$Dives)
sp<- as.data.frame(scores(spmds, "species"))

grdi <- spdt$Depth
ordi <- ordisurf(spmds, as.numeric(grdi), plot = FALSE, bs = 'ds')
ordi.grid <- ordi$grid
ordi.sp <- expand.grid(x = ordi.grid$x, y = ordi.grid$y)
ordi.sp$z <- as.vector(ordi.grid$z)
ordi.sp.na <- data.frame(na.omit(ordi.sp))

NMDS <- data.frame(x = spmds$point[ ,1], y = spmds$point[ ,2], 
                   Type = spdt$Dives,tr=spdt$Depth)
NMDS$name <- rownames(NMDS)
NMDS.sp <- data.frame(MDS1=spmds$spec[ ,1],MDS2=spmds$spec[ ,2])
NMDS.sp$name <- rownames(NMDS.sp)
#Taxa_classfi <- Taxa_classfi[rownames(NMDS.sp),]
#NMDS.sp <- data.frame(NMDS.sp,Taxa_classfi)

mytheme <- theme(panel.grid = element_blank(), 
                 legend.key = element_blank(), 
                 panel.background = element_rect(color = 'black', fill = 'transparent'),
                 axis.title = element_text(size = 16,colour = "black"),
                 axis.text = element_text(size = 16,colour = "black"),
                 axis.text.x = element_text(colour = "black",size = 15),
                 axis.text.y = element_text(colour = "black",size = 15),
                 strip.text = element_blank(),
                 text = element_text(family = "sans",size = 16))



a <- ggplot()+
  stat_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..)) +
  geom_point(data = NMDS, aes(x, y, fill = factor(Type)), pch = 21, size = 5) +
  scale_fill_brewer(palette = "Set1") +
  geom_text_contour(data = ordi.sp.na, aes(x, y, z = z, color = ..level..), size = 6) +
  scale_colour_continuous(high = 'red', low = 'black') +
  geom_hline(yintercept=0,linetype=3,size=0.01)+ 
  geom_vline(xintercept=0,linetype=3,size=0.01)+
  labs(x = 'NMDS1', y = 'NMDS2', color = 'Depth(m)') +
  mytheme+
  labs(fill="Dives")+
  annotate("text",x=-1.5,y=0.5,hjust=0,vjust=0,
           label=paste("Stress:",round(spmds$stress,digits = 2)),size=6)+
  theme(legend.position = c(0.14,0.14),
        legend.direction = "vertical", legend.box = "horizontal",
        legend.spacing.x=unit(0.1, 'mm'),
        legend.key.size=unit(7, 'pt'))

a

####存为pdf
pdf("NMDS.pdf")

a

garbage <- dev.off()


#βdiversity=====================================

 
library(ade4)
library(adegraphics)
library(adespatial)
library(vegan)
library(vegetarian)
library(ggplot2)
library(FD)
library(taxize)
library(ggpubr)
library(ggpmisc)
library(reshape) 

source("panelutils.R")
source("Rao.R")

#读入文件

env <- env[,1:10]


# Jaccard-based Podani indices (Abundance data)
#fish.pod.j <- beta.div.comp(spe, coef = "J", quant = FALSE) 
fish.pod.j <- beta.div.comp(spe, coef = "J", quant = TRUE)  
# What is in the output object?
summary(fish.pod.j)
# Display summary statistics:
fish.pod.j$part


fish.rich <- as.matrix(fish.pod.j$rich)
fish.repl <- as.matrix(fish.pod.j$repl)
fish.jac <- as.matrix(fish.pod.j$D)


# Plot of the Jaccard, replacement and richness difference indices 
# between nearest neighbours
# First, extract the subdiagonals of the square dissimilarity 
# matrices
fish.repl.neigh <- diag(fish.repl[-1, ]) # Replacement
fish.rich.neigh <- diag(fish.rich[-1, ]) # Richness difference
fish.jac.neigh <- diag(fish.jac[-1, ]) # Jaccard Dj index


#####（Jaccard）

jac <- data.frame(row.names(env),env$depth,env$meanVelocity)
jac <- jac[-1,]
jac <- data.frame(jac,fish.jac.neigh)
names(jac) <- c("samples","depth","meanVelocity","jaccard")

abun <- data.frame(row.names(env),env$depth)
abun <- abun[-1,]
abun <- data.frame(abun,fish.rich.neigh)
names(abun) <- c("samples","depth","abun")

turn <- data.frame(row.names(env),env$depth)
turn <- turn[-1,]
turn <- data.frame(turn,fish.repl.neigh)
names(turn) <- c("samples","depth","turn")

all <- data.frame(turn,abun$abun,jac$jaccard)

names(all) <- c("samples","depth","turn","abun","jaccard")

m <- melt(all, id = c("samples", "depth")) 


m$variable <-
  factor(m$variable,
         levels = c("jaccard", "abun", "turn"))

####GAM

library(nlme)
library(mgcv)

a2 <- ggplot(data=m,aes(x=depth,y=value,color=variable))+
  geom_point(data=m,aes(x=depth,y=value,color=variable),size=3,alpha=0.3)+
  geom_smooth(method='gam',formula=y~s(x),se=TRUE,size=2,alpha=0.3,fill="#999999")+
  theme(panel.background=element_blank(),                      
        panel.border = element_rect(colour="black",fill=NA),
        axis.text = element_text(size=14,color="black"),
        axis.title = element_text(size=15),
        legend.position = "top",
        legend.key = element_blank(),
        legend.key.size = unit(13,"pt"),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13), 
  )+
  ylab("s(depth)")+
  xlab("Depth")+
  scale_color_discrete(
    name = "", 
    labels = c("Jaccard", "Abundance difference", "Turnover")
  )

a2


#####β diversity(0-1)=====================================



# Jaccard-based Podani indices (Abundance data)
fish.pod.j <- beta.div.comp(spe, coef = "J", quant = FALSE) 
#fish.pod.j <- beta.div.comp(spe, coef = "J", quant = TRUE)  
# What is in the output object?
summary(fish.pod.j)
# Display summary statistics:
fish.pod.j$part


fish.rich <- as.matrix(fish.pod.j$rich)
fish.repl <- as.matrix(fish.pod.j$repl)
fish.jac <- as.matrix(fish.pod.j$D)


# Plot of the Jaccard, replacement and richness difference indices 
# between nearest neighbours
# First, extract the subdiagonals of the square dissimilarity 
# matrices
fish.repl.neigh <- diag(fish.repl[-1, ]) # Replacement
fish.rich.neigh <- diag(fish.rich[-1, ]) # Richness difference
fish.jac.neigh <- diag(fish.jac[-1, ]) # Jaccard Dj index

jac <- data.frame(row.names(env),env$depth)
jac <- jac[-1,]
jac <- data.frame(jac,fish.jac.neigh)
names(jac) <- c("samples","depth","jaccard")

rich <- data.frame(row.names(env),env$depth)
rich <- rich[-1,]
rich <- data.frame(rich,fish.rich.neigh)
names(rich) <- c("samples","depth","rich")

turn <- data.frame(row.names(env),env$depth)
turn <- turn[-1,]
turn <- data.frame(turn,fish.repl.neigh)
names(turn) <- c("samples","depth","turn")

all <- data.frame(turn,rich$rich,jac$jaccard)

names(all) <- c("samples","depth","turn","rich","jaccard")

m <- melt(all, id = c("samples", "depth")) 


m$variable <-
  factor(m$variable,
         levels = c("jaccard", "rich", "turn"))

####GAM

library(nlme)
library(mgcv)

a3 <- ggplot(data=m,aes(x=depth,y=value,color=variable))+
  geom_point(data=m,aes(x=depth,y=value,color=variable),size=3,alpha=0.3)+
  geom_smooth(method='gam',formula=y~s(x),se=TRUE,size=2,alpha=0.3,fill="#999999")+
  theme(panel.background=element_blank(),                      
        panel.border = element_rect(colour="black",fill=NA),
        axis.text = element_text(size=14,color="black"),
        axis.title = element_text(size=15),
        legend.position = "top",
        legend.key = element_blank(),
        legend.key.size = unit(13,"pt"),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13), 
  )+
  ylab("s(depth)")+
  xlab("Depth")+
  scale_color_discrete(
    name = "", 
    labels = c("Jaccard", "Richness difference", "Turnover")
  )

a3



library(ggpubr)

a <- ggarrange( a2, a3, 
                labels = c("A", "B"),
                ncol = 1)
a



pdf("beta diversity.pdf")

a

garbage <- dev.off()


#GAM to analyze the beta diversity pattern of the extension gradient 
#and test for differences in the beta diversity pattern among the three parallel survey lines along the depth gradient

library(ade4)
library(adegraphics)
library(adespatial)
library(vegan)
library(vegetarian)
library(ggplot2)
library(FD)
library(taxize)
library(reshape)
library(dplyr)


source("panelutils.R")
source("Rao.R")


spe <- read.csv("spe.csv",row.names = 1,header = T) 
env <- read.csv("env.csv",row.names = 1,header = T)


spe <- subset(spe, select = colSums(spe) != 0)
spe <- subset(spe, select = colSums(spe) != 1)
spe <- spe[,1:72]

nonzero.spe <- which(rowSums(spe)>0)

spe <- spe[nonzero.spe,]
env <- env[nonzero.spe,]

ord.depth <- order(env$depth)

spe <- spe[ord.depth,]
env <- env[ord.depth,]

colnames(spe)
colnames(env)



A.spe <- which(env$group=="A")
speA <- spe[A.spe,]
envA <- env[A.spe,]

spe <- speA
env <- envA

ord.depth <- order(envA$depth)
#ord.depth <- order(env$substrate)
spe2 <- spe[ord.depth,]
env2 <- env[ord.depth,]
dim(spe2)
# Replacement, richness and nestedness indices of the data ===

# Jaccard-based Podani indices (presence-absence data)
#fish.pod.j <- beta.div.comp(spe, coef = "J", quant = FALSE) 
fish.pod.j <- beta.div.comp(spe2, coef = "J", quant = TRUE)  
# What is in the output object?
summary(fish.pod.j)
# Display summary statistics:
fish.pod.j$part

fish.rich <- as.matrix(fish.pod.j$rich)
fish.repl <- as.matrix(fish.pod.j$repl)
fish.jac <- as.matrix(fish.pod.j$D)


# Plot of the Jaccard, replacement and richness difference indices 
# between nearest neighbours
# First, extract the subdiagonals of the square dissimilarity 
# matrices
fish.repl.neigh <- diag(fish.repl[-1, ]) # Replacement
fish.rich.neigh <- diag(fish.rich[-1, ]) # Richness difference
fish.jac.neigh <- diag(fish.jac[-1, ]) # Jaccard Dj index

#write.csv(fish.repl.neigh,file='output/neighrepl.csv')
#write.csv(fish.rich.neigh,file='output/neighrich.csv')
#write.csv(fish.jac.neigh,file='output/neighjac.csv')

########Abandance-DiveA


jacA <- data.frame(row.names(env2),env2$depth)
jacA <- jacA[-1,]
a <- c(rep("A", 15))
jacA <- data.frame(jacA,fish.jac.neigh,a)
names(jacA) <- c("samples","depth","jaccard","dives")

abunA <- data.frame(row.names(env2),env2$depth)
abunA <- abunA[-1,]
abunA <- data.frame(abunA,fish.rich.neigh,a)
names(abunA) <- c("samples","depth","abun","dives")

turnA <- data.frame(row.names(env2),env2$depth)
turnA <- turnA[-1,]
turnA <- data.frame(turnA,fish.repl.neigh,a)
names(turnA) <- c("samples","depth","turn","dives")

########Abandance-DiveB

spe <- read.csv("spe.csv",row.names = 1,header = T) 
env <- read.csv("env.csv",row.names = 1,header = T)


spe <- subset(spe, select = colSums(spe) != 0)
spe <- subset(spe, select = colSums(spe) != 1)
spe <- spe[,1:72]

nonzero.spe <- which(rowSums(spe)>0)

spe <- spe[nonzero.spe,]
env <- env[nonzero.spe,]

ord.depth <- order(env$depth)

spe <- spe[ord.depth,]
env <- env[ord.depth,]

colnames(spe)
colnames(env)


B.spe <- which(env$group=="B")
speB <- spe[B.spe,]
envB <- env[B.spe,]

spe <- speB
env <- envB


ord.depth <- order(envB$depth)
#ord.depth <- order(env$substrate)
spe2 <- spe[ord.depth,]
env2 <- env[ord.depth,]
dim(spe2)
# Replacement, richness and nestedness indices of the fish data ===
# Jaccard-based Podani indices (presence-absence data)
#fish.pod.j <- beta.div.comp(spe, coef = "J", quant = FALSE) 
fish.pod.j <- beta.div.comp(spe2, coef = "J", quant = TRUE)  
# What is in the output object?
summary(fish.pod.j)
# Display summary statistics:
fish.pod.j$part


fish.rich <- as.matrix(fish.pod.j$rich)
fish.repl <- as.matrix(fish.pod.j$repl)
fish.jac <- as.matrix(fish.pod.j$D)


# Plot of the Jaccard, replacement and richness difference indices 
# between nearest neighbours
# First, extract the subdiagonals of the square dissimilarity 
# matrices
fish.repl.neigh <- diag(fish.repl[-1, ]) # Replacement
fish.rich.neigh <- diag(fish.rich[-1, ]) # Richness difference
fish.jac.neigh <- diag(fish.jac[-1, ]) # Jaccard Dj index


jacB <- data.frame(row.names(env2),env2$depth)
jacB <- jacB[-1,]
b <- c(rep("B", 13))
jacB <- data.frame(jacB,fish.jac.neigh,b)
names(jacB) <- c("samples","depth","jaccard","dives")

abunB <- data.frame(row.names(env2),env2$depth)
abunB <- abunB[-1,]
abunB <- data.frame(abunB,fish.rich.neigh,b)
names(abunB) <- c("samples","depth","abun","dives")

turnB <- data.frame(row.names(env2),env2$depth)
turnB <- turnB[-1,]
turnB <- data.frame(turnB,fish.repl.neigh,b)
names(turnB) <- c("samples","depth","turn","dives")


########Abandance-DiveF

spe <- read.csv("spe.csv",row.names = 1,header = T) 
env <- read.csv("env.csv",row.names = 1,header = T)


spe <- subset(spe, select = colSums(spe) != 0)
spe <- subset(spe, select = colSums(spe) != 1)
spe <- spe[,1:72]

nonzero.spe <- which(rowSums(spe)>0)

spe <- spe[nonzero.spe,]
env <- env[nonzero.spe,]

ord.depth <- order(env$depth)

spe <- spe[ord.depth,]
env <- env[ord.depth,]

colnames(spe)
colnames(env)



F.spe <- which(env$group=="F")
speF <- spe[F.spe,]
envF <- env[F.spe,]

spe <- speF
env <- envF


ord.depth <- order(envF$depth)
#ord.depth <- order(env$substrate)
spe2 <- spe[ord.depth,]
env2 <- env[ord.depth,]
dim(spe2)
# Replacement, richness and nestedness indices of the fish data ===

# Jaccard-based Podani indices (presence-absence data)
#fish.pod.j <- beta.div.comp(spe, coef = "J", quant = FALSE) 
fish.pod.j <- beta.div.comp(spe2, coef = "J", quant = TRUE)  
# What is in the output object?
summary(fish.pod.j)
# Display summary statistics:
fish.pod.j$part


fish.rich <- as.matrix(fish.pod.j$rich)
fish.repl <- as.matrix(fish.pod.j$repl)
fish.jac <- as.matrix(fish.pod.j$D)


# Plot of the Jaccard, replacement and richness difference indices 
# between nearest neighbours
# First, extract the subdiagonals of the square dissimilarity 
# matrices
fish.repl.neigh <- diag(fish.repl[-1, ]) # Replacement
fish.rich.neigh <- diag(fish.rich[-1, ]) # Richness difference
fish.jac.neigh <- diag(fish.jac[-1, ]) # Jaccard Dj index


jacF <- data.frame(row.names(env2),env2$depth)
jacF <- jacF[-1,]
f <- c(rep("F", 25))
jacF <- data.frame(jacF,fish.jac.neigh,f)
names(jacF) <- c("samples","depth","jaccard","dives")

abunF <- data.frame(row.names(env2),env2$depth)
abunF <- abunF[-1,]
abunF <- data.frame(abunF,fish.rich.neigh,f)
names(abunF) <- c("samples","depth","abun","dives")

turnF <- data.frame(row.names(env2),env2$depth)
turnF <- turnF[-1,]
turnF <- data.frame(turnF,fish.repl.neigh,f)
names(turnF) <- c("samples","depth","turn","dives")


####Jaccard=====================================
library(dplyr)
m <- bind_rows(jacA, jacB,jacF)


m$dives <-
  factor(m$dives,
         levels = c("A", "B", "F"))



b1 <- ggplot(data=m,aes(x=depth,y=jaccard,color=dives))+
  geom_point(data=m,aes(x=depth,y=jaccard,color=dives),size=1,alpha=0.3)+
  geom_smooth(method='gam',formula=y~s(x),se=TRUE,size=1,alpha=0.3,fill="#999999")+
  theme(panel.background=element_blank(),                      
        panel.border = element_rect(colour="black",fill=NA),
        axis.text = element_text(size=10,color="black"),
        axis.title = element_text(size=10),
        axis.title.x = element_text(angle = 180),
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5),
        axis.text.y = element_text(angle = 90,hjust=0.5),
        legend.position = c(0.9,0.2),
        legend.key = element_blank(),
        legend.key.size = unit(5,"pt")
  )+
  ylab("Dissimilarity")+
  xlab("Depth")+
  coord_cartesian(xlim=c(500, 2000),ylim=c(-0.12, 1.1))+
  scale_x_continuous(breaks=seq(500, 2000, 200))+
  scale_y_continuous(breaks=seq(0, 1, 0.2))+
  scale_color_discrete(
    name = "", 
    labels = c("A", "B", "F"),
    guide = guide_legend(
      direction = "horizontal",
      title.position = "top",
      label.position = "top",
      label.hjust = 0.5,
      label.vjust = 1,
      label.theme = element_text(size=5,angle = 90)
    )
  )+
  annotate(geom = "text", x = 490, y =-0.1 , label = "A",size=6,
           angle = 90)

b1


ggsave(filename = "77-jaccard.png",
       plot = b1,
       width = 9,
       height = 6,
       units = "cm",
       dpi = 300)


####abun=====================================

n <- bind_rows(abunA, abunB,abunF)



n$dives <-
  factor(n$dives,
         levels = c("A", "B", "F"))



b2 <- ggplot(data=n,aes(x=depth,y=abun,color=dives))+
  geom_point(data=n,aes(x=depth,y=abun,color=dives),size=1,alpha=0.3)+
  geom_smooth(method='gam',formula=y~s(x),se=TRUE,size=1,alpha=0.3,fill="#999999")+
  theme(panel.background=element_blank(),                      
        panel.border = element_rect(colour="black",fill=NA),
        axis.text = element_text(size=10,color="black"),
        axis.title = element_text(size=10),
        axis.title.x = element_text(angle = 180),
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5),
        axis.text.y = element_text(angle = 90,hjust=0.5),
        legend.position = c(0.9,0.2),
        legend.key = element_blank(),
        legend.key.size = unit(5,"pt")
  )+
  ylab("Abundance Difference")+
  xlab("Depth")+
  coord_cartesian(xlim=c(500, 2000),ylim=c(-0.12, 1.1))+
  scale_x_continuous(breaks=seq(500, 2000, 200))+
  scale_y_continuous(breaks=seq(0, 1, 0.2))+
  scale_color_discrete(
    name = "", 
    labels = c("A", "B", "F"),
    guide = guide_legend(
      direction = "horizontal",
      title.position = "top",
      label.position = "top",
      label.hjust = 0.5,
      label.vjust = 1,
      label.theme = element_text(size=5,angle = 90)
    )
  )+
  annotate(geom = "text", x = 490, y =-0.1 , label = "B",size=6,
           angle = 90)

b2


ggsave(filename = "77-abun.png",
       plot = b2,
       width = 9,
       height = 6,
       units = "cm",
       dpi = 300)

####turn=====================================

p <- bind_rows(turnA, turnB,turnF)



p$dives <-
  factor(p$dives,
         levels = c("A", "B", "F"))



b3 <- ggplot(data=p,aes(x=depth,y=turn,color=dives))+
  geom_point(data=p,aes(x=depth,y=turn,color=dives),size=1,alpha=0.3)+
  geom_smooth(method='gam',formula=y~s(x),se=TRUE,size=1,alpha=0.3,fill="#999999")+
  theme(panel.background=element_blank(),                      
        panel.border = element_rect(colour="black",fill=NA),
        axis.text = element_text(size=10,color="black"),
        axis.title = element_text(size=10),
        axis.title.x = element_text(angle = 180),
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5),
        axis.text.y = element_text(angle = 90,hjust=0.5),
        legend.position = c(0.9,0.2),
        legend.key = element_blank(),
        legend.key.size = unit(5,"pt")
  )+
  ylab("Turnover")+
  xlab("Depth")+
  coord_cartesian(xlim=c(500, 2000),ylim=c(-0.12, 1.1))+
  scale_x_continuous(breaks=seq(500, 2000, 200))+
  scale_y_continuous(breaks=seq(0, 1, 0.2))+
  scale_color_discrete(
    name = "", 
    labels = c("A", "B", "F"),
    guide = guide_legend(
      direction = "horizontal",
      title.position = "top",
      label.position = "top",
      label.hjust =0.5,
      label.vjust = 1,
      label.theme = element_text(size=5,angle = 90)
    )
  )+
  annotate(geom = "text", x = 490, y =-0.1 , label = "C",size=6,
           angle = 90)

b3

 
ggsave(filename = "77-turn.png",
       plot = b3,
       width = 9,
       height = 6,
       units = "cm",
       dpi = 1200)



library(ggplot2)
library(cowplot)

b4 <- plot_grid(b3,b2,b1,ncol = 1, align = 'v')

b4

ggsave(filename = "Abundance- parallel survey lines β diversity.png",
       plot = b4,
       width = 9,
       height = 17,
       units = "cm",
       dpi = 300)

#######0-1  is similar to abundance 
 

######dbRDA=====================================

# Load packages, functions and data 
library(ade4)
library(adegraphics)
library(adespatial)
library(cocorresp)
library(vegan)
library(vegan3d)
library(ape)   # For PCoA with Lingoes correction (not in the book)
library(MASS)
library(ellipse)
library(FactoMineR)
library(rrcov)
library(rdacca.hp)

# Source additional functions that will be used later in this
# Chapter. Our scripts assume that files to be read are in
# the working directory.
source("hcoplot.R")
source("triplot.rda.R")
source("plot.lda.R")
source("polyvars.R")
source("screestick.R")

#读入文件
spe <- read.csv("spe.csv",row.names = 1,header = T) 
env <- read.csv("env.csv",row.names = 1,header = T)


spe <- subset(spe, select = colSums(spe) != 0)
spe <- subset(spe, select = colSums(spe) != 1)
spe <- spe[,1:72]

nonzero.spe <- which(rowSums(spe)>0)

spe <- spe[nonzero.spe,]
env <- env[nonzero.spe,]

ord.depth <- order(env$depth)

spe <- spe[ord.depth,]
env <- env[ord.depth,]

colnames(spe)
colnames(env)
env <- env[,1:10]

cluster <- read.csv("cluster.csv",row.names = 1,header = T)

#spe <- read.csv("Spe6.csv",row.names = 1,header = T)
#env <- read.csv("env6.csv",row.names = 1,header = T)

#substrates
# Recode the substrate variable into a factor (qualitative) 
# variable to show how these are handled in the ordinations
substrate2 <- rep(".soft", nrow(env))
#substrate2[env$substrate == 1] <- ".soft"
substrate2[env$substrate == 2] <- ".mod_soft"
substrate2[env$substrate == 3] <- ".mixed"
substrate2[env$substrate == 4] <- ".mod_hard"
substrate2[env$substrate == 5] <- ".hard"

substrate2 <- factor(substrate2, 
                     levels = c(".soft", ".mod_soft", ".mixed", ".mod_hard", ".hard"))
table(substrate2)
# Create an env3 data frame with slope as a qualitative variable
env3 <- env
env3$substrate <- substrate2

#env4 <- data.frame(env3,cluster)

#尝试归一化
env.scale <- scale(env,center = T,scale = T)
env4 <- subset(env.scale,select = -c(substrate))
env4<- data.frame(env4,env3$substrate)

env3 <- env4

colnames(env3)[colnames(env3) == 'env3.substrate'] <- 'substrate'
#log(x+1)
# Forward selection using vegan's ordistep()
spe.rda.all <- 
  capscale(log1p(spe) ~., data=env3, 
           distance = "bray", 
           add = "lingoes", 
           comm = spe)

mod0 <- 
  capscale(log1p(spe) ~1, data=env3, 
           distance = "bray", 
           add = "lingoes", 
           comm = spe)

step.forward <- 
  ordistep(mod0, 
           scope = formula(spe.rda.all), 
           direction = "forward", 
           #trace=FALSE,
           permutations = how(nperm = 499)
  )
RsquareAdj(step.forward)

step.forward$anova

#selected env variables: all
spe.rda.pars <- 
  capscale(log1p(spe) ~ meanVelocity + broadBPI + depth + curvature + fineBPI + substrate + slope +backscatter,#手动加的backscatter
           data=env3, 
           distance = "bray", 
           method="db",
           add = "lingoes", 
           comm = spe)
summary(spe.rda.pars,scaling=1)


#perc <- round(100*(summary(spe.rda.pars)$cont$importance[2, 1:2]), 2)
#perc

#### CAP1:0.3603   CAP2:0.1873   CAP3:0.1050


sc_si <- scores(spe.rda.pars, display="sites", choices=c(1,2), scaling=1)
sc_sp <- scores(spe.rda.pars, display="species", choices=c(1,2), scaling=1)
sc_bp <- scores(spe.rda.pars, display="bp", choices=c(1, 2), scaling=1)
sc_centr <-read.csv("center.csv",row.names = 1,header = T)


sc_si2 <- data.frame(sc_si,cluster)


row.names(sc_sp)
sel.sp=c(41,15,50,42,29,4,7,9,10,13,33,63,69,70,49,26,25,20)
sc_sp[sel.sp,]

#select environmental variables
sel.bp <- c(1,2,3,4,5,10,11)
sc_bp[sel.bp,]

##### 1.2-----
library(vegan)
library(ggplot2)
library(ggrepel)
library(ggforce)

windowsFonts(A=windowsFont("Times New Roman"),
             B=windowsFont("Arial"))

p <- ggplot() +
  geom_point(data = sc_si2,aes(CAP1*15,CAP2*15,color=sc_si2$cluster),size=2)+
  scale_shape_manual(values = c(1:7))+
  scale_color_manual(name="Cluster",values=c(2:7,"grey"))+
  scale_fill_manual(values=c(2:7,"grey"))+
  theme(#text = element_text(family = "B"),
    panel.background=element_blank(),                      
    panel.border = element_rect(colour="black",fill=NA),
    axis.text = element_text(size=12,color="black"),
    axis.title = element_text(size=14),
    legend.position = c(0.1,0.16),
    legend.key = element_blank(),
    legend.key.size = unit(0.1,"cm"),
    legend.title = element_text(size=10),
    legend.text = element_text(size=10)
  )+
  geom_segment(data = as.data.frame(sc_sp[sel.sp,]),aes(x = 0, y = 0, xend = CAP1, yend = CAP2), 
               arrow = arrow(angle=22.5,length = unit(0.4,"cm")),
               linetype=1,size=0.4,colour = "red")+
  geom_text_repel(data = as.data.frame(sc_sp[sel.sp,]),aes(CAP1,CAP2,label=row.names(sc_sp[sel.sp,])),size=4,colour = "red")+
  geom_segment(data = as.data.frame(sc_bp[sel.bp,]),aes(x = 0, y = 0, xend = CAP1*25, yend = CAP2*25), 
               arrow = arrow(angle=22.5,length = unit(0.4,"cm")),
               linetype=1, size=0.4,colour = "blue")+
  geom_text_repel(data =as.data.frame(sc_bp[sel.bp,]),aes(CAP1*25,CAP2*25,label=row.names(sc_bp[sel.bp,])),size=4,colour = "blue")+
  labs(x="dbRDA1(36.03%)",y="dbRDA2(18.73%)")+
  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())+
  geom_point(data = sc_centr,aes(CAP1*15,CAP2*15),size=2,color="blue",shape=18)+
  geom_text_repel(data = as.data.frame(sc_centr),aes(CAP1*15,CAP2*15,label=row.names(sc_centr)),size=4,colour = "blue")

p

##
pdf("RDA13.pdf")

p

garbage <- dev.off()

#######2.3-------

sc_si23 <- scores(spe.rda.pars, display="sites", choices=c(2,3), scaling=1)
sc_sp23 <- scores(spe.rda.pars, display="species", choices=c(2,3), scaling=1)
sc_bp23 <- scores(spe.rda.pars, display="bp", choices=c(2, 3), scaling=1)
sc_centr <-read.csv("center.csv",row.names = 1,header = T)


sc_si2 <- data.frame(sc_si23,cluster)


row.names(sc_sp23)
sel.sp23=c(41,15,50,42,29,4,7,9,10,13,33,63,69,70,49,26,25)
sc_sp23[sel.sp23,]

#select environmental variables
sel.bp23 <- c(1,2,3,4,5,10,11)
sc_bp23[sel.bp23,]


library(vegan)
library(ggplot2)
library(ggrepel)
library(ggforce)

p23 <- ggplot() +
  geom_point(data = sc_si2,aes(CAP2*15,CAP3*15,color=sc_si2$cluster),size=2)+
  scale_shape_manual(values = c(1:7))+
  scale_color_manual(name="Cluster",values=c(2:7,"grey"))+
  scale_fill_manual(values=c(2:7,"grey"))+
  theme(panel.background=element_blank(),                      
        panel.border = element_rect(colour="black",fill=NA),
        axis.text = element_text(size=12,color="black"),
        axis.title = element_text(size=14),
        legend.position = c(0.1,0.16),
        legend.key = element_blank(),
        legend.key.size = unit(0.1,"cm"),
        legend.title = element_text(size=10),
        legend.text = element_text(size=10)
  )+
  geom_segment(data = as.data.frame(sc_sp23[sel.sp23,]),aes(x = 0, y = 0, xend = CAP2, yend = CAP3), 
               arrow = arrow(angle=22.5,length = unit(0.4,"cm")),
               linetype=1, size=0.6,colour = "red")+
  geom_text_repel(data = as.data.frame(sc_sp23[sel.sp23,]),aes(CAP2,CAP3,label=row.names(sc_sp23[sel.sp23,])),size=4,colour = "red")+
  geom_segment(data = as.data.frame(sc_bp23[sel.bp23,]),aes(x = 0, y = 0, xend = CAP2*65, yend = CAP3*65), 
               arrow = arrow(angle=22.5,length = unit(0.4,"cm")),
               linetype=1, size=0.4,colour = "blue")+
  geom_text_repel(data =as.data.frame(sc_bp23[sel.bp23,]),aes(CAP2*65,CAP3*65,label=row.names(sc_bp23[sel.bp23,])),size=4,colour = "blue")+
  labs(x="dbRDA2(18.73%)",y="dbRDA3(10.50%)")+
  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=)+
  geom_point(data = sc_centr,aes(CAP2*15,CAP3*15),size=2,color="blue",shape=18)+
  geom_text_repel(data = as.data.frame(sc_centr),aes(CAP2*15,CAP3*15,label=row.names(sc_centr)),size=4,colour = "blue")

p23


pdf("RDA23.pdf")

p23

garbage <- dev.off()