# Huang WC, Liu Y, Zhang X, Zhang CJ, Zou D, Zheng S, Xu W, Luo Z, Liu F, Li M. 
#Comparative genomic analysis reveals metabolic flexibility of Woesearchaeota. 
#Nat Commun. 2021 Sep 6;12(1):5281. doi: 10.1038/s41467-021-25565-9
library(ggplot2)
library(ggmap)
library(sp)
library(maptools)
library(maps)
woese<-read.table('Data4.tsv',sep='\t',row.names=1,header=T,check.name=FALSE)
colors<-c('#d6274a','#4ea94a','#f6d534','#4862ab','#ea7a3a','#7e3a8e','#75c7d0','#b0569c','#f5bbbe','#b9cf39','#955f25')
#ncol<-c('#b0569c','#ea7a3a','#75c7d0','#b9cf39','#7e3a8e','#4ea94a','#f6d534','#4862ab','#955f25','#f5bbbe','#d6274a')
names(colors)<-unique(woese$biotope)
mp<-NULL
mapworld<-borders("world",colour = "gray50",fill="white") #
mp2<-mp+geom_jitter(data=woese,aes(x=longitude_deg,y=latitude_deg,size=woese$`relative abundance of Woesearchaeota (%)`,color=woese$biotope)) +
  scale_color_manual(values=colors) +
  xlab('longitude') +
  ylab('latitude')
mp2
mp2<-mp2+theme(legend.position='bottom')
mp2
library(Rtsne)
library(vegan)
wunifrac<-read.table('unweighted_unifrac_otu_woese_beta.tsv',sep='\t',header=T,row.names = 1)
wunifrac[,ncol(wunifrac)]
wunifrac.m<-as.matrix(wunifrac)
tsne_model<-Rtsne(wunifrac.m,num_threads = 38,check_duplicates = F)
df_tsne<-as.data.frame(tsne_model$Y)
rownames(df_tsne)<-rownames(wunifrac)
p<-ggplot(df_tsne,aes(x=V1,y=V2))+geom_point()
meta<-woese
p
salinity<-meta[which(rownames(meta)==rownames(df_tsne)),'empo_2']
salinity<-meta[rownames(df_tsne),'empo_2']
biotope<-meta[which(rownames(meta)==rownames(df_tsne)),'biotope']
biotope<-meta[rownames(df_tsne),'biotope']
df_tsne$saline<-salinity
df_tsne$biotope<-biotope

##
##PERMANOVA on unifrac using empo_2 and biotope

empo_2<-adonis(wunifrac.m~salinity)

biotope_test<-adonis(wunifrac.m~biotope)

library(grid)
grob2 <- grobTree(textGrob("p-value=0.001\nR^2=0.325", x=0.75,  y=0.95, hjust=0,
                           gp=gpar(col="black", fontface="italic")))
grob1<-grobTree(textGrob("p-value=0.001\nR^2=0.138", x=0.75,  y=0.95, hjust=0,
                         gp=gpar(col="black", fontface="italic")))

p<-ggplot(df_tsne,aes(x=V1,y=V2,color=saline))+geom_point()+xlab('tSNE 1')+ylab('tSNE 2')+annotation_custom(grob1)
#p<-p+theme(legend.position='none')
p1<-ggplot(df_tsne,aes(x=V1,y=V2,color=biotope))+geom_point()+scale_color_manual(values = colors)+xlab('tSNE 1')+ylab('tSNE 2')+annotation_custom(grob2)
#p1<-p1+theme(legend.position='none')

##alpha diversity
alpha_div<-read.table('alpha_div_pd_whole_tree.txt',sep='\t',header=T,row.names = 1)
alpha_div$biotope<-meta[rownames(alpha_div),'biotope']
p0<-ggplot(alpha_div,aes(x=biotope,y=PD_whole_tree,fill=biotope))+
  geom_boxplot()+
  scale_fill_manual(values = colors)+
  guides(fill=FALSE)+
  theme(axis.text.x=element_text(angle=45,hjust = 1,vjust=1))
## Show data points if n less than and equal 10
saline_env_alpha<-alpha_div[which(alpha_div$biotope=='saline environment'),]
marine_water_alpha<-alpha_div[which(alpha_div$biotope=='marine water'),]
saltmarshes_alpha<-alpha_div[which(alpha_div$biotope=='Saltmarshes'),]
p0<-p0+
  geom_point(data=saline_env_alpha,aes(x=biotope,y=PD_whole_tree))+
  geom_point(data=marine_water_alpha,aes(x=biotope,y=PD_whole_tree))+
  geom_point(data=saltmarshes_alpha,aes(x=biotope,y=PD_whole_tree))

library(gridExtra)
lay1 <- rbind(
  c(1,1),
  c(2,3))
plots<-list(p0,p,p1)
final_plot<-grid.arrange(grobs=plots,layout_matrix=lay1)
#ggsave('G:/woese_revision/fig1/fig1_revised.pdf',final_plot,height =10 ,width=14,dpi=300,useDingbats=FALSE)
#ggsave('G:/woese_revision/fig1/fig1c.pdf',p,height =5 ,width=5,dpi=300,useDingbats=FALSE)

###Post hoc test
a<-aov(PD_whole_tree~biotope,data=alpha_div)
tk<-TukeyHSD(a)
###Post hoc abudance test
abu_t<-aov(woese$`relative abundance of Woesearchaeota (%)`~biotope,data=woese)
tK_abu_t<-TukeyHSD(abu_t)

##Mantel 
simple_cor<-meta[8:15]
simple_cor$Woesearchaeota<-meta[17]
simple_cor<-simple_cor[rownames(alpha_div),]
mantel_f<-function(otu,env){
  library(vegan)
  library(dplyr)
  vars<-colnames(env)
  models<-list()
  for (i in seq_along(vars)){
    if (vars[i]!='Woesearchaeota'){
      env_sub<-na.omit(env[vars[i]])
      otu_sub<-otu[rownames(env_sub),rownames(env_sub)]
      otu_dist<-as.matrix(otu_sub)
      dis<-vegdist(env_sub,method='euclidean')
      model<-vegan::mantel(otu_dist,dis,permutations = 999)
      name<-vars[i]
      num<-length(rownames(env_sub))
      statistic<-model$statistic
      signif<-model$signif
      models[[i]]<-data.frame(name=name,statistic=statistic,signif=signif,number=num,row.names=NULL)
    }
  }
  models %>% bind_rows()
}

info<-mantel_f(wunifrac,simple_cor)

library(gridExtra)
lay1 <- rbind(
  c(1,1),
  c(2,2),
  c(3,4))
plots<-list(mp2,p0,p,p1)
final_plot<-grid.arrange(grobs=plots,layout_matrix=lay1)
ggsave('fig1.pdf',final_plot,height =18,width=13,units = "in")