library(readr)
library(ggtree)
library(treeio)
library(DESeq2)
library(ggplot2)
library(stringr)
library(tidyverse)

df_handle<-function(
  otu_table,
  treefile,
  tiplabelfile,
  pattern){
  dat01<-read_tsv(otu_table)
  tiplabgroup<-read_tsv(tiplabelfile)
  tree<-read.tree(treefile)
  new.tree<-full_join(tree,tiplabgroup,by='label')
  ggtree(new.tree)+
    #geom_tiplab(align = TRUE)+
    geom_tippoint(aes(color=GROUP))+
    #xlim(NA,1.5)+
    theme(legend.position = "none")-> p1
  
  group.info.01<-data.frame(
    row.names = dat01 %>% select(matches(pattern)) %>% 
      colnames(),
    conditions = c(rep("trt",5),
                   rep("ctl",5))
  )
  group.info.01$conditions<-factor(group.info.01$conditions)
  dat01 %>% 
    select(matches(paste(pattern,"OTU ID",sep = "|"))) %>% 
    column_to_rownames('OTU ID') -> dat01.1
  dds <- DESeq2::DESeqDataSetFromMatrix(countData = dat01.1,
                                        colData=group.info.01,
                                        design = ~ conditions)
  dds_res <- DESeq2::DESeq(dds, sfType = "poscounts")
  
  res <- results(dds_res, 
                 tidy=T, 
                 format="DataFrame",
                 contrast = c("conditions","trt","ctl"))
  res %>% 
    mutate(str_extract(row,"f__\\S+;") %>% 
             str_split_fixed(pattern = ";",n=10) %>% 
             as.data.frame() %>% 
             select(1)) -> res.01
  
  
  res.01[which(res.01$V1%in%tree$tip.label),] %>% 
    filter(!is.na(log2FoldChange)) %>% 
    group_by(V1) %>% 
    summarise(log2FoldChange=mean(log2FoldChange)) -> res.02
  merge(res.02,tiplabgroup,by.x='V1',by.y="label") -> res.02
  res.02$V1<-factor(res.02$V1,
                    levels = p1$data %>% na.omit() %>% arrange(y) %>% pull(label))
  dat01.1 %>% rownames_to_column(var = "family") %>% 
    mutate(str_extract(family,"f__\\S+;") %>% 
             str_split_fixed(pattern = ";",n=10) %>% 
             as.data.frame() %>% 
             select(1)) -> dat01.2
  dat01.1 %>% rownames_to_column(var = "family") %>% 
    mutate(str_extract(family,"f__\\S+;") %>% 
             str_split_fixed(pattern = ";",n=10) %>% 
             as.data.frame() %>% 
             select(1)) %>% 
    merge(tiplabgroup,by.x="V1",by.y='label') %>% 
    select(-c('family','GROUP')) %>% 
    group_by(V1) %>% 
    summarise_all(mean) %>% 
    mutate(mrd=rowSums(.[2:11])/10) %>% 
    select(V1,mrd) -> dat01.2
  res.02$V1<-str_replace(res.02$V1,"f__","")
  res.02$GROUP<-str_replace(res.02$GROUP,"c__","")
  dat01.2$V1<-str_replace(dat01.2$V1,"f__","")
  new.tree@phylo$tip.label<-
    str_replace(new.tree@phylo$tip.label,
                "f__","")
  new.tree@data$GROUP<-
    str_replace(new.tree@data$GROUP,
                "c__","")
  return(list(bubble.1 = res.02,
              bubble.2 = dat01.2,
              tree.1 = new.tree))
}




bubble_figure<-function(dat){
  ggplot(data=dat,
         aes(x=log2FoldChange,y=V1))+
    geom_vline(xintercept = 0,lty="dashed",
               color="grey",size=1)+
    geom_point(aes(size=log2FoldChange,
                   fill=GROUP),
               shape=21)+
    guides(size='none',
           fill=guide_legend(ncol=6))+
    scale_fill_discrete(name=NULL)+
    #ggtitle("F_0_20VSW_0_20")+
    theme_minimal()+
    theme(legend.position = "bottom",
          legend.direction = "vertical",
          legend.justification = c(0,0),
          axis.ticks.y = element_blank(),
          axis.text.y = element_blank(),
          panel.grid = element_blank(),
          axis.line.x = element_line(),
          axis.ticks.x = element_line())+
    labs(y=NULL) -> p
  return(p)
}
#bubble_figure(plot_data[["bubble.1"]])

dat01<-read_tsv("data/20220414/data/otu_taxon_otu.full.xls")
dat01 %>% colnames()
dat01 %>% select(matches("[FW]\\d_0_20")) %>% colnames()
dat01 %>% select(matches("[FW]\\d_20_40")) %>% colnames()
dat01 %>% select(matches("[FW]\\d_40_60")) %>% colnames()

dat01 %>% select(matches("[HW]\\d_0_20")) %>% colnames()
dat01 %>% select(matches("[HW]\\d_20_40")) %>% colnames()
dat01 %>% select(matches("[HW]\\d_40_60")) %>% colnames()

dat01 %>% select(matches("[MW]\\d_0_20")) %>% colnames()
dat01 %>% select(matches("[MW]\\d_20_40")) %>% colnames()
dat01 %>% select(matches("[MW]\\d_40_60")) %>% colnames()

df_handle(otu_table = "data/20220414/data/otu_taxon_otu.full.xls",
          treefile = "data/20220414/data/phylo_tree .tre",
          tiplabelfile = "data/20220414/data/species_group .xls",
          pattern = "[FW]\\d_0_20") -> plot_data

trt_ctl<-c("[FW]\\d_0_20","[FW]\\d_20_40","[FW]\\d_40_60",
           "[HW]\\d_0_20","[HW]\\d_20_40","[HW]\\d_40_60",
           "[MW]\\d_0_20","[MW]\\d_20_40","[MW]\\d_40_60")

plot_list = list()
for (i in 1:9){
  df_handle(otu_table = "data/20220414/data/otu_taxon_otu.full.xls",
            treefile = "data/20220414/data/phylo_tree .tre",
            tiplabelfile = "data/20220414/data/species_group .xls",
            pattern = trt_ctl[i]) -> plot_data
  #print(plot_data)
  #print(bubble_figure(plot_data[["bubble.1"]]))
  bubble_figure(plot_data[["bubble.1"]]) -> plot_list[[i]]
}

p1<-ggtree(plot_data[["tree.1"]])+
  geom_tiplab(align = T,offset = 0.01)+
  geom_tippoint(aes(color=GROUP))+
  xlim(NA,2)+
  theme(legend.position = "none")
p1
p1+plot_spacer()+
  (plot_list[[1]]+
  plot_list[[2]]+
  plot_list[[3]]+
  plot_list[[4]]+
  plot_list[[5]]+
  plot_list[[6]]+
  plot_list[[7]]+
  plot_list[[8]]+
  plot_list[[9]]+
  plot_layout(guides = "collect",
              nrow=1)&theme(legend.position = "bottom"))+
  plot_layout(widths = c(1,1,9))

image.png