#1）Calculation of alpha diversity:

install.packages("ggplot2")
install.packages("ggpubr")
install.packages("vegan")
library(ggplot2)
library(ggpubr)
library(vegan)
otu=read.table('C:/Rtest/Yang/2022B.txt',row.names = 1,header=T,sep='\t')
otu=t(otu)
shannon=diversity(otu,"shannon",base =exp(1))
simpson=diversity(otu,"simpson")
alpha=data.frame(shannon,simpson,check.names=T)

write.table(alpha,"C:/Rtest/Yang/2022d.xls",sep='\t',quote=F)
map<-read.table('D:/r-working/mapping_file.txt',row.names = 1,header=T,comment.char ='',sep='\t',check.names=F)
group<-map["Group1"]
alpha<-alpha[match(rownames(group),rownames(alpha)),]
data<-data.frame(group,alpha,check.rows=T)
p=ggplot(data=data,aes(x=Group1,y=shannon))+geom_boxplot(fill=rainbow(7)[2])

#2）Redundancy analysis (RDA)
spe = read.table("spe.xls", header=T, row.names=1, sep="\t", comment.char="",stringsAsFactors = F,quote = "")
dim(spe)
head(spe)
gt = read.table("gene.xls", header=T, row.names=1, sep="\t", comment.char="",stringsAsFactors = F,quote = "")
dim(gt)
head(gt)

ENV = read.table("env.xls", header=T, row.names=1, sep="\t", comment.char="",stringsAsFactors = F,quote = "")
dim(ENV)
head(ENV)
decorana(t(spe))
decorana(t(gt))

rda=rda(t(spe),ENV[4:14],scale=T)
rda=rda(t(gt),ENV[4:14],scale=T)

rda.data=data.frame(rda$CCA$u[,1:2],ENV[1:3])
head(rda.data)
rda.spe=data.frame(rda$CCA$v[,1:2])
head(rda.spe)
rda.env=data.frame(rda$CCA$ biplot[,1:2])
head(rda.env)
rda1 =round(rda$CCA$eig[1]/sum(rda$CCA$eig)*100,2)
rda2 =round(rda$CCA$eig[2]/sum(rda$CCA$eig)*100,2)


cols=c("CK"="grey","NT"="red","PT"="blue","WL"="green")
cols
rda.plot=ggplot(data=rda.data,aes(RDA1,RDA2))+
  geom_point(aes(color=tillage,fill=tillage,shape=depth),size=3)+
  scale_color_manual(values=cols)+scale_fill_manual(values = cols)+
  labs(title="RDA plot",x=paste("RDA1",rda1," %"),y=paste("RDA2",rda2," %"))+
  theme_bw()+
  theme(axis.title = element_text(family = "serif", face = "bold", size = 18,colour = "black"))+
  theme(axis.text = element_text(family = "serif", face = "bold", size = 16,color="black"))+
  theme(panel.grid=element_blank())+
  geom_hline(yintercept=0)+geom_vline(xintercept=0)+
  geom_segment(data=rda.env,aes(x=0,y=0,xend=rda.env[,1],yend=rda.env[,2]),size=1,colour="black",arrow=arrow(angle = 35,length=unit(0.3,"cm")))+
  geom_text(data=rda.env,aes(x=rda.env[,1],y=rda.env[,2],label=rownames(rda.env)),size=3.5,      colour="black", hjust=(1-sign(rda.env[,1]))/2,angle=(180/pi)*atan(rda.env[,2]/rda.env[,1]))+
  theme(legend.position="top")
rda.plot


rda.plot+geom_mark_ellipse(aes(fill=tillage),alpha=0.1,tol=0.6,expand = unit(0, "mm")) + coord_cartesian(clip = "off")

ggsave(filename = "RDA.plot.pdf",plot = rda.plot,device="pdf",width = 8,height = 8)

#3) Principal component analysis （PCA）：

data<-iris
head(data)

dt<-as.matrix(scale(data[,1:4]))
head(dt)

rm1<-cor(dt)
rs1<-eigen(rm1)
val <- rs1$values

(Standard_deviation <- sqrt(val))

(Proportion_of_Variance <- val/sum(val))
(Cumulative_Proportion <- cumsum(Proportion_of_Variance))

par(mar=c(6,6,2,2))
plot(rs1$values,type="b",
     cex=2,
     cex.lab=2,
     cex.axis=2,
     lty=2,
     lwd=2,
     xlab = "PCA",
     ylab="Eigenvalue",)

(U<-as.matrix(rs1$vectors))

PC <-dt %*% U
colnames(PC) <- c("PC1","PC2","PC3","PC4")
head(PC)

df<-data.frame(PC,iris$Species)
head(df)

library(ggplot2)

xlab<-paste0("PC1(",round(Proportion_of_Variance[1]*100,2),"%)")
ylab<-paste0("PC2(",round(Proportion_of_Variance[2]*100,2),"%)")

p1<-ggplot(data = df,aes(x=PC1,y=PC2,color=iris.Species))+
  stat_ellipse(aes(fill=iris.Species),
               type ="norm", geom ="polygon",alpha=0.2,color=NA)+
  geom_point()+labs(x=xlab,y=ylab,color="")+
  guides(fill=F)

library(scatterplot3d)
color = c(rep('purple',50),rep('orange',50),rep('blue',50))
scatterplot3d(df2[,1:3],color=color,
              pch = 16,angle=30,
              box=T,type="p",
              lty.hide=2,lty.grid = 2)
legend("topleft",c('Setosa','Versicolor','Virginica'),
       fill=c('purple','orange','blue'),box.col=NA)

4）Correlation analysis 
install.packages("corrplot")
install.packages("ggplot2")
install.packages("ggpubr")
library(corrplot)
library(ggplot2)
library(ggpubr)
td <-read.table("testdata.txt", header=T)
library(corrplot)
cor (td, method="pearson")
cor(x, y = NULL, use = "everything",
method = c("pearson", "kendall", "spearman"))
tdc <- cor (td, method="pearson")
corrplot(tdc)
corrplot(tdc)
corrplot(tdc, method = "ellipse", 
         type = "upper",
         tl.col = "black", tl.cex = 1.2, tl.srt = 45
         )

method = c("circle", "square", "ellipse", "number", "shade", "color", "pie"),
type = c("full", "lower", "upper"),
corrplot(tdc, method = "ellipse", type = "upper",
         tl.col = "black", tl.cex = 0.8, tl.srt = 45,tl.pos = "lt")
corrplot(tdc, method = "number", type = "lower",
         tl.col = "n", tl.cex = 0.8, tl.pos = "n",
         add = T)

#5) Radar map：
# install.packages("devtools")
devtools::install_github("ricardo-bion/ggradar")
library(ggradar)
  df<-read.csv("C:/Rtest/Yang/leida.csv", header = T)
df
df[, 1] <- paste0("G", 1:6)
colnames(df) <- c("treatment", paste0("B", 1:9))

ggradar(plot.data =df, 
        background.circle.colour="white", 
        #axis.labels=c("R","a","b","c","d","e","f","g","h","i")
        gridline.min.colour = "blue", #The color of the circle with the smallest axis value
        gridline.mid.colour="red", #The color of the axis value median circle
        gridline.max.colour="black", #The color of the maximum axis value circle
        gridline.min.linetype = 1, #Line type for minimum circle with axis value
        gridline.mid.linetype = 2, #The line type of the axis value median circle
        gridline.max.linetype = 3, #Line type for maximum axis value circle
        group.line.width=3  #Radar chart line width
)
ggradar(df,
        background.circle.colour = "#D7D6D1",
        gridline.min.colour = "grey",
        gridline.mid.colour="grey",
        gridline.max.colour="black",
        centre.y = 0,
        grid.min = 5,
        grid.mid = 20,
        grid.max = 35,

        values.radar = c(5, 20, 35),
        group.line.width=1.5,
        group.point.size=4)

# 6)Nightingale Rose Chart
rm(list=ls())
setwd("C:/Rtest/Yang/Rose4")
install.packages("ggprism")
install.packages("ggthemes")
library(ggplot2)
library(ggprism)
library(ggthemes)

df <- read.table("wm.txt",header = T, check.names = F)
df
p<-ggplot(df, aes(x = sample, y = value, fill = sample)) +
  geom_bar(stat = "identity", color = "white",
           lwd = 1, show.legend = FALSE,width = 0.8)+
  scale_fill_manual(values=c("#E74A32","#4CBCD6","#00A188","#F39C80","#8592B5","#92D2C3","#D00000","#7F6146","#8592B5"))+
  ylim(0,35)+theme(axis.text.x = element_text(size = 10,face = "bold"))+
   theme_gdocs()
p+coord_polar()

df1 <- read.table("SpWM.txt",header = T, check.names = F)
df1
p1<-ggplot(df1, aes(x = sample, y = value, fill = sample)) +
  geom_bar(stat = "identity", color = "white",
           lwd = 1, show.legend = FALSE,width = 0.8)+
  scale_fill_manual(values=c("#E74A32","#4CBCD6","#00A188","#F39C80","#8592B5","#92D2C3","#D00000","#7F6146","#8592B5"))+
  ylim(0,35)+theme(axis.text.x = element_text(size = 10,face = "bold"))+
  theme_gdocs()
p1+coord_polar()

df2 <- read.table("PWM.txt",header = T, check.names = F)
df2
p2<-ggplot(df2, aes(x = sample, y = value, fill = sample)) +
  geom_bar(stat = "identity", color = "white",
           lwd = 1, show.legend = FALSE,width = 0.8)+
  scale_fill_manual(values=c("#E74A32","#4CBCD6","#00A188","#F39C80","#8592B5","#92D2C3","#D00000","#7F6146","#8592B5"))+
  ylim(0,35)+theme(axis.text.x = element_text(size = 10,face = "bold"))+
  theme_gdocs()
p2+coord_polar()



df3 <- read.table("SWM.txt",header = T, check.names = F)
df3
p3<-ggplot(df3, aes(x = sample, y = value, fill = sample)) +
  geom_bar(stat = "identity", color = "white",
           lwd = 1, show.legend = FALSE,width = 0.8)+
  scale_fill_manual(values=c("#E74A32","#4CBCD6","#00A188","#F39C80","#8592B5","#92D2C3","#D00000","#7F6146","#8592B5"))+
  ylim(0,35)+theme(axis.text.x = element_text(size = 10,face = "bold"))+
  theme_gdocs()
p3+coord_polar()


df4 <- read.table("SmWM.txt",header = T, check.names = F)
df4
p4<-ggplot(df4, aes(x = sample, y = value, fill = sample)) +
  geom_bar(stat = "identity", color = "white",
           lwd = 1, show.legend = FALSE,width = 0.8)+
  scale_fill_manual(values=c("#E74A32","#4CBCD6","#00A188","#F39C80","#8592B5","#92D2C3","#D00000","#7F6146","#8592B5"))+
  ylim(0,35)+theme(axis.text.x = element_text(size = 10,face = "bold"))+
  theme_gdocs()
p4+coord_polar()


df5 <- read.table("RSWM.txt",header = T, check.names = F)
df5
p5<-ggplot(df5, aes(x = sample, y = value, fill = sample)) +
  geom_bar(stat = "identity", color = "white",
           lwd = 1, show.legend = FALSE,width = 0.8)+
  scale_fill_manual(values=c("#E74A32","#4CBCD6","#00A188","#F39C80","#8592B5","#92D2C3","#D00000","#7F6146","#8592B5"))+
  ylim(0,35)+theme(axis.text.x = element_text(size = 10,face = "bold"))+
  theme_gdocs()
p5+coord_polar()