# wilcox.test -------------------------------------------------------------
library(multtest)
##soil characterstics
#pH
#wilcox.test
pH_rural<-c(7.808,8.708,8.351,7.978)
pH_urban<-c(8.17,8.356,8.023,8.141,8.092,8.135,7.863)
pH_suburban<-c(8.48,8.249,8.237,8.06,8.316,8.312,8.301,8.038,7.999,8.166,8.118,8.391,8.044)
wilcox_pH_su <- wilcox.test(pH_suburban,pH_urban,exact=FALSE)
wilcox_pH_su
wilcox_pH_ru <- wilcox.test(pH_rural,pH_urban,exact=FALSE)
wilcox_pH_ru
wilcox_pH_sr <- wilcox.test(pH_suburban,pH_rural,exact=FALSE)
wilcox_pH_sr

#SHN
#wilcox.test
SHN_rural<-c(125.624,119.61,83.632,121.568)
SHN_urban<-c(72.083,168.268,104.636,91.217,109.327,76.144,109.325)
SHN_suburban<-c(481,133.61,137.085,363.427,98.912,52.902,117.348,65.974,130.855,121.561,56.39,53.721,71.353)
wilcox_SHN_su <- wilcox.test(SHN_suburban,SHN_urban,exact=FALSE)
wilcox_SHN_su
wilcox_SHN_ru <- wilcox.test(SHN_rural,SHN_urban,exact=FALSE)
wilcox_SHN_ru
wilcox_SHN_sr <- wilcox.test(SHN_suburban,SHN_rural,exact=FALSE)
wilcox_SHN_sr

#SAP
#wilcox.test
SAP_rural<-c(8.679,10.652,11.802,7.904)
SAP_urban<-c(15.733,9.172,10.048,8.943,10.889,8.567,7.836)
SAP_suburban<-c(10.1,9.214,10.617,5.18,8.828,8.041,10.547,9.841,9.46,8.142,8.978,7.647,8.924)
wilcox_SAP_su <- wilcox.test(SAP_suburban,SAP_urban,exact=FALSE)
wilcox_SAP_su
wilcox_SAP_ru <- wilcox.test(SAP_rural,SAP_urban,exact=FALSE)
wilcox_SAP_ru
wilcox_SAP_sr <- wilcox.test(SAP_suburban,SAP_rural,exact=FALSE)
wilcox_SAP_sr


##seed density
#spring
#wilcox.test
spring_rural<-c(533,1693,280,913)
spring_urban<-c(407,933,1140,213,487)
spring_suburban<-c(2380,1033,867,1213,607,527,1113,733,367,880,433)
wilcox_spring_su <- wilcox.test(spring_suburban,spring_urban,exact=FALSE)
wilcox_spring_su
wilcox_spring_ru <- wilcox.test(spring_rural,spring_urban,exact=FALSE)
wilcox_spring_ru
wilcox_spring_sr <- wilcox.test(spring_suburban,spring_rural,exact=FALSE)
wilcox_spring_sr


#autumn
#wilcox.test
autumn_rural<-c(2275,4651,1494,1528)
autumn_urban<-c(1154,951,1460,1664,13513,2852,2139)
autumn_suburban<-c(951,1935,4244,2580,951,1086,577,5500,1120,4108,2309)
wilcox_autumn_su <- wilcox.test(autumn_suburban,autumn_urban,exact=FALSE)
wilcox_autumn_su
wilcox_autumn_ru <- wilcox.test(autumn_rural,autumn_urban,exact=FALSE)
wilcox_autumn_ru
wilcox_autumn_sr <- wilcox.test(autumn_suburban,autumn_rural,exact=FALSE)
wilcox_autumn_sr


# NMDS --------------------------------------------------------------------
library(vegan)
library(ggplot2)

#autumn
dat1 <- read.csv("Supplemental Data S2.csv",row.names = 1)
dat1 <- dat1[, -(23:26)]
dat2 <- data.frame(t(dat1))
dat3<-dat2[,apply(dat2>0,2,sum)>2]
dat <-log(dat3+1)
rich<- read.csv("Supplemental Data S3.csv",row.names = 1)

#NMDS
nmds_dat <- metaMDS(dat, distance = 'euclidean', k = 2)
nmds_dat$stress

nmds_dat_site <- data.frame(nmds_dat$points)

nmds_dat_species <- data.frame(nmds_dat$species)

###ggplot2
#add group
nmds_dat_site$name <- rownames(nmds_dat_site)
nmds_dat_site$gradient <- c(rep('urban', 7), rep('suburban', 11), rep('rural', 4))
nmds_dat_site$richness <- rich$richness
nmds_dat_site$site <- rich$site


p <- ggplot(data = nmds_dat_site, aes(MDS1, MDS2)) +
  geom_point(aes(color= gradient,size=richness,shape=gradient,fill=gradient)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  scale_shape_manual(values = c(21, 21, 21)) +
  scale_fill_manual(values = c('black', 'grey', 'white')) +
  scale_size_area(max_size = 20) +
  theme(panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5,size=40, family="serif"), 
        legend.position = 'right',
        legend.text = element_text(size = 40, colour = "black",family="serif"),
        legend.title = element_text(size = 40, colour = "black",family="serif"),
        legend.key.size = unit(45, "pt"),
        legend.box.margin = margin(10,10,10,10, unit = "pt"),
        axis.text.x = element_text(size=40, family="serif", colour = "black"),
        axis.text.y = element_text(size=40, family="serif", colour = "black"),
        axis.line = element_line(linewidth = 1, colour = "black"),
        axis.ticks = element_line(linewidth = 1),
        axis.ticks.length.y = unit(.25, "cm"),
        axis.ticks.length.x = unit(.25, "cm"),
        axis.title.x = element_text(size=40, family="serif", colour = "black"),
        axis.title.y = element_text(size=40, family="serif", colour = "black")) +
  labs(x = 'NMDS1', y = 'NMDS2', title = 'autumn') +
  annotate('text', label = paste('Stress =',round(nmds_dat$stress, 4)), x = -3.5, y = 2.2, size = 10, colour = 'black', family="serif")	#标注应力函数值

p

#spring
library(vegan)
library(ggplot2)

dat1 <- read.csv("Supplemental Data S4.csv",row.names = 1)
dat1 <- dat1[, -(21:24)]
dat2 <- data.frame(t(dat1))
dat3<-dat2[,apply(dat2>0,2,sum)>2]
dat <-log(dat3+1)
rich<- read.csv("Supplemental Data S5.csv",row.names = 1)

#NMDS
nmds_dat <- metaMDS(dat, distance = 'euclidean', k = 2)
nmds_dat$stress


nmds_dat_site <- data.frame(nmds_dat$points)
nmds_dat_species <- data.frame(nmds_dat$species)

###ggplot2
nmds_dat_site$name <- rownames(nmds_dat_site)
nmds_dat_site$gradient <- c(rep('urban', 5), rep('suburban', 11), rep('rural', 4))
nmds_dat_site$richness <- rich$richness
nmds_dat_site$site <- rich$site


p <- ggplot(data = nmds_dat_site, aes(MDS1, MDS2)) +
  geom_point(aes(color= gradient,size=richness,shape=gradient,fill=gradient)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  scale_shape_manual(values = c(21, 21, 21)) +
  scale_fill_manual(values = c('black', 'grey', 'white')) +
  scale_size_area(max_size = 20) +
  theme(panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5,size=40, family="serif",), 
        legend.position = 'right',
        legend.text = element_text(size = 40, colour = "black",family="serif"),
        legend.title = element_text(size = 40, colour = "black",family="serif"),
        legend.key.size = unit(45, "pt"),
        axis.text.x = element_text(size=40, family="serif", colour = "black"),
        axis.text.y = element_text(size=40, family="serif", colour = "black"),
        axis.line = element_line(linewidth = 1, colour = "black"),
        axis.ticks = element_line(linewidth = 1),
        axis.ticks.length.y = unit(.25, "cm"),
        axis.ticks.length.x = unit(.25, "cm"),
        axis.title.x = element_text(size=40, family="serif", colour = "black"),
        axis.title.y = element_text(size=40, family="serif", colour = "black")) +
  labs(x = 'NMDS1', y = 'NMDS2', title = 'spring') +
  annotate('text', label = paste('Stress =',round(nmds_dat$stress, 4)), x = -4, y = 2.8, size = 10, colour = 'black', family="serif")	#标注应力函数值

p



# Boxplot -----------------------------------------------------------------

library(vegan)
library(ggplot2)
library(multtest)

#richness
#spring

dat1 <- read.csv("Supplemental Data S6.csv")
boxplot_spring <- boxplot(richness~gradient,dat1[dat1$season=="spring",])
boxplot_spring <- boxplot(sqrt(richness)~gradient,dat1[dat1$season=="spring",])
#wilcox.test
spring_rural<-c(10,17,9,11)
spring_urban<-c(7,14,14,6,12)
spring_suburban<-c(11,14,10,16,14,10,14,12,10,9,10)
wilcox_spring_su <- wilcox.test(spring_suburban,spring_urban,exact=FALSE)
wilcox_spring_su
wilcox_spring_ru <- wilcox.test(spring_rural,spring_urban,exact=FALSE)
wilcox_spring_ru
wilcox_spring_sr <- wilcox.test(spring_suburban,spring_rural,exact=FALSE)
wilcox_spring_sr

#autumn

boxplot_autumn <- boxplot(richness~gradient,dat1[dat1$season=="autumn",])
boxplot_autumn <- boxplot(sqrt(richness)~gradient,dat1[dat1$season=="autumn",])
#wilcox.test
autumn_rural<-c(10,14,8,9)
autumn_urban<-c(9,6,5,7,13,5,8)
autumn_suburban<-c(8,7,18,11,4,6,3,7,7,8,8)
wilcox_autumn_su <- wilcox.test(autumn_suburban,autumn_urban,exact=FALSE)
wilcox_autumn_su
wilcox_autumn_ru <- wilcox.test(autumn_rural,autumn_urban,exact=FALSE)
wilcox_autumn_ru
wilcox_autumn_sr <- wilcox.test(autumn_suburban,autumn_rural,exact=FALSE)
wilcox_autumn_sr


#plot
ggplot(dat1, aes(x=gradient, y=richness,fill=season)) +
  geom_boxplot(width=c(0.4),lwd=1.2,outlier.size=8,position = position_dodge(0.5)) +
  stat_boxplot(geom = "errorbar",width=0.2,lwd=1.2,position = position_dodge(0.5),
               aes(x=gradient, y=richness)) +
  guides(color="none") +
  labs(x="gradient",y="richness") +
  annotate("text",x=1.15:3.15, y=19.62, label="a", cex=20, family="serif")+
  annotate("text",x=0.9:3, y=19.62, label="A", cex=20, family="serif")+
  theme_classic () +
  theme(plot.title=element_text(hjust=0.5,size=45, family="serif", colour = "black"),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        legend.key.size = unit(50, "pt"),
        legend.text = element_text(size=45, family="serif", colour = "black"),
        legend.title = element_text(size=45, family="serif", colour = "black"),
        axis.text.x = element_text(size=45, family="serif", colour = "black"),
        axis.text.y = element_text(size=45, family="serif", colour = "black"),
        axis.line = element_line(linewidth = 1.3, colour = "black"),
        axis.ticks = element_line(linewidth = 1.3),
        axis.ticks.length.y = unit(.5, "cm"),
        axis.ticks.length.x = unit(.5, "cm"),
        axis.title.x = element_text(size=45, family="serif", colour = "black"),
        axis.title.y = element_text(size=45, family="serif", colour = "black"))




#shannon
dat1 <- read.csv("Supplemental Data S7.csv")

#spring
boxplot_spring <- boxplot(shannon~gradient,dat1[dat1$season=="spring",])
boxplot_spring <- boxplot(sqrt(shannon)~gradient,dat1[dat1$season=="spring",])
#wilcox.test
spring_rural<-c(2.15872938,2.36259219,1.99452149,1.69991082)
spring_urban<-c(1.7771392,2.26728855,2.42108702,1.50807607,2.28106461)
spring_suburban<-c(0.80060649,2.2201494,1.79457813,2.32478906,2.27108047,
                   1.96819122,2.02539472,1.90635838,2.17141748,1.79874667,
                   1.95731695)
wilcox_spring_su <- wilcox.test(spring_suburban,spring_urban,exact=FALSE)
wilcox_spring_su
wilcox_spring_ru <- wilcox.test(spring_rural,spring_urban,exact=FALSE)
wilcox_spring_ru
wilcox_spring_sr <- wilcox.test(spring_suburban,spring_rural,exact=FALSE)
wilcox_spring_sr

#autumn
boxplot_autumn <- boxplot(shannon~gradient,dat1[dat1$season=="autumn",])
boxplot_autumn <- boxplot(sqrt(shannon)~gradient,dat1[dat1$season=="autumn",])
#wilcox.test
autumn_rural<-c(1.78324273,2.42409947,1.71162407,1.9140931)
autumn_urban<-c(2.0220181,1.54885811,1.23324791,1.49561028,1.15445441,
                1.17085651,1.50603095)
autumn_suburban<-c(1.82723376,1.62605712,2.50096193,2.09133084,1.03004876,
                   1.41201134,0.97184877,1.15499414,1.61402675,1.59794497,
                   1.79501279)
wilcox_autumn_su <- wilcox.test(autumn_suburban,autumn_urban,exact=FALSE)
wilcox_autumn_su
wilcox_autumn_ru <- wilcox.test(autumn_rural,autumn_urban,exact=FALSE)
wilcox_autumn_ru
wilcox_autumn_sr <- wilcox.test(autumn_suburban,autumn_rural,exact=FALSE)
wilcox_autumn_sr

#plot
ggplot(dat1, aes(x=gradient, y=shannon,fill=season)) +
  geom_boxplot(width=c(0.4),lwd=1.2,outlier.size=8,position = position_dodge(0.5)) +
  stat_boxplot(geom = "errorbar",width=0.2,lwd=1.2,position = position_dodge(0.5),
               aes(x=gradient, y=shannon)) +
  guides(color="none") +
  labs(x="gradient",y="Shannon") +
  annotate("text",x=1.15:3.15, y=2.7, label="a", cex=20, family="serif")+
  annotate("text",x=0.9, y=2.7, label="A",cex=20, family="serif")+
  annotate("text",x=1.88, y=2.7, label="A",cex=20, family="serif")+
  annotate("text",x=2.88, y=2.7, label="B",cex=20, family="serif")+
  theme_classic () +
  theme(plot.title=element_text(hjust=0.5,size=45, family="serif", colour = "black"),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        legend.key.size = unit(50, "pt"),
        legend.text = element_text(size=45, family="serif", colour = "black"),
        legend.title = element_text(size=45, family="serif", colour = "black"),
        axis.text.x = element_text(size=45, family="serif", colour = "black"),
        axis.text.y = element_text(size=45, family="serif", colour = "black"),
        axis.line = element_line(linewidth = 1.3, colour = "black"),
        axis.ticks = element_line(linewidth = 1.3),
        axis.ticks.length.y = unit(.5, "cm"),
        axis.ticks.length.x = unit(.5, "cm"),
        axis.title.x = element_text(size=48, family="serif", colour = "black"),
        axis.title.y = element_text(size=48, family="serif", colour = "black"))



##Heatmap####
##autumn
library(ComplexHeatmap)
library(ggplot2)
library(colorRamp2)

install.packages("plotrix")
library(plotrix)

dat1 <- read.csv("Supplemental Data S8.csv")
dim(dat1)
rownames(dat1)=dat1$Species

mat<-dat1
mat <- dat1[,2:43]
mat <- as.matrix(mat)

png("1128.png",units="in", width=15, height=15,res=1000)

col_fun <- colorRamp2(c(0, 1,1000,2000, 10000),
                      c("#EFEEEE","#ffffff","#f239AE", "#f31BA4", "#f311A0"))

Heatmap(mat,name = "Density(Seed/m2)",col=col_fun,
        #km = 4, 
        #top_annotation = ha,
        cluster_columns = F, 
        cluster_rows = F,
        show_row_names = T,
        row_names_side = "left",
        row_names_gp = gpar(fontsize = 10, col = "black", 
                            fontface = "italic",family="serif"),
        show_column_names = T,
        show_row_dend = F,
        #column_km =3,
        show_column_dend = F) +
  Heatmap(dat1$Sourcetype,
          name = "Source type", col= c("#F2ABAB","#ADCFFA"),
          width = unit(5, "mm")) +
  Heatmap(dat1$Lifeform, 
          name = "Life form",
          width = unit(5, "mm"),
          col = c( "#C4ADFB","#78C857","#8D6F64"))+
  Heatmap(dat1$Dispersalmodes,
          name = "Dispersal model", 
          width = unit(5, "mm"),
          col = c( "#EFC3CA","#FBD53B","#9DC6FA","#ADABAB"))+
  Heatmap(dat1$Growthtype,
          name = "Growth type",
          col= c("#CC6CE7","#348BF8","#FE9900","#f38181",
                 "#a5dee5","#e0f9b5","#FDDA54","#ffcfdf"),
          width = unit(5, "mm")) 


dev.off()