library(dplyr)
library(dunn.test)
library(tibble)
library(reshape2)
library(ggplot2)

corr_multi<-function(table, genes1, genes2, font_size=5, ruta="./", exportar=F)
{
  cor.list<-list()
  cont<-1
  for (gene1 in genes1){
    for (gene2 in genes2){
      print(gene1)
      print(gene2)
      cor_temp<-cor.test(table[,gene1], table[,gene2])
      if (!is.na(cor_temp$p.value)){
        pval_plot<-paste("cor =",format(cor_temp$estimate, digits=2),if(cor_temp$p.value < 2.2e-16){" , p < 2.2e-16"}else{paste(", p=",format(cor_temp$p.value,digits=3))})
      }
      cor.list[[cont]]<-ggplotGrob(ggplot(data=table, aes(x=table[,gene1], y=table[,gene2]))+ 
                                     geom_point()+
                                     geom_smooth(method="lm")+
                                     geom_text(aes(x=min(table[,gene1], na.rm=T), y=(max(table[,gene2],na.rm = T)+max(table[,gene2],na.rm = T)*0.2), label=pval_plot), hjust="inward", family="serif", size=font_size)+
                                     labs(x=gene1, y=gene2)+
                                     theme_bw())
      if (exportar == T){
        png(paste0(ruta,gene1,"vs",gene2,".png"), res=900, width=3000, height=3000)
        plot(cor.list[[cont]])
        dev.off()
      }
      cont<-cont+1
    }
  }
  return(cor.list)
}

multi_stats<-function(table, value.var, x, group, stat.test, adjust="default", paired=F){
  ## Requires dplyr and tibble packets
  defaults=c("dunn"="none", "ttest"="holm", "wilcox"="holm")
  funs<-c("ttest"="pairwise.t.test", "wilcox"="pairwise.wilcox.test")
  if (adjust == "default"){adjust=defaults[stat.test]}
  stat.def<-as.data.frame(matrix(nrow=0, ncol=5))
  colnames(stat.def)<-c(x, "group1", "group2", "p.adj", "p.signif")
  for (point in unique(table[,x])){
    condition<-all(table %>% filter(table[,x] == point) %>% pull(value.var) == 0) == F
    len_group<-length(unique(table %>% filter(table[,x] == point) %>% pull(group)))
    if (condition == T & !is.na(condition) & len_group > 1){
      if(stat.test == "dunn"){
        test<-dunn.test(table %>% filter(table[,x] == point) %>% pull(value.var), table %>% filter(table[,x] == point) %>% pull(group), method=adjust)
        comp<-strsplit(test$comparisons, " - ")
        stat.temp<-data.frame(matrix(unlist(comp), nrow=length(comp), byrow=T), "p.adj"=test$P.adjusted, stringsAsFactors = F, check.names = F)
        colnames(stat.temp)[1:2]<-c("group1", "group2")
      }else if (stat.test %in% names(funs)){
        test<-get(funs[stat.test])(table %>% filter(table[,x] == point) %>% pull(value.var), table %>% filter(table[,x] == point) %>% pull(group), method=adjust, paired=paired)$p.value
        stat.temp<-melt(test)
        colnames(stat.temp)<-c("group1", "group2","p.adj")
      }
      
      stat.temp["p.signif"]<-case_when(
        stat.temp$p.adj >= 0.05 ~ "ns",
        stat.temp$p.adj < 0.0001 ~ "****",
        stat.temp$p.adj < 0.001 ~ "***",
        stat.temp$p.adj < 0.01 ~ "**",
        stat.temp$p.adj < 0.05 ~ "*"
      )
      stat.temp<-stat.temp %>% add_column(x=point, .before=T)
      colnames(stat.temp)[1]<-x
      stat.def<-rbind(stat.def, stat.temp)
    }
  }
  stat.def["Method"]<-stat.test
  return(stat.def)
}

generate_labstats<-function(table_stat, table, value.var, x, group, y="max", bracket.offset=0.05, bracket.length=0.02){
  table[,group]<-as.factor(table[,group])
  table[,x]<-as.factor(table[,x])
  se<-function(x, na.rm=F) sd(x, na.rm = na.rm)/sqrt(length(x))
  if (y == "max"){ 
    formula<-as.formula(paste0(colnames(table_stat)[1], "~."))
    agg<-dcast(table, formula, value.var = value.var, fun.aggregate = max, na.rm=T)
  }else if (y == "mean"){ 
    formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
    agg<-dcast(table, formula, value.var = value.var, fun.aggregate = mean, na.rm=T)
    agg<- data.frame(x=agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
    colnames(agg)[1]<-x
  }else if (y == "mean+sd"){ 
    formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
    agg<- dcast(table, formula, value.var = value.var, fun.aggregate = function(x) mean(x,na.rm=T)+sd(x,na.rm=T))
    agg<- data.frame(x=agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
    colnames(agg)[1]<-x
  }else if (y == "mean+se"){ 
    formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
    agg<- dcast(table, formula, value.var = value.var, fun.aggregate = function(x) mean(x,na.rm=T)+se(x,na.rm=T))
    agg<- data.frame(agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
    colnames(agg)[1]<-x
  }
  
  t<-data.frame("y1"=merge(table_stat, agg ,sort=F)[,"."]+diff(range(table[value.var], na.rm = T))*bracket.offset,
                "y2"=merge(table_stat, agg ,sort=F)[,"."]+diff(range(table[value.var], na.rm = T))*bracket.offset,
                "x1"= match(table_stat[,x], unique(table[,x]))+
                  0.75*((match(table_stat$group1, levels(table[,group]))-0.5)/length(levels(table[,group]))-0.5),
                "x2"= match(table_stat[,x], unique(table[,x]))+
                  0.75*((match(table_stat$group2, levels(table[,group]))-0.5)/length(levels(table[,group]))-0.5)
  )
  for (dia in unique(table_stat[,1])){
    t[table_stat[,x] == dia,"y1"]<-seq(t[table_stat[,x] == dia,"y1"][1],
                                           t[table_stat[,x] == dia,"y1"][1]+diff(range(table[,value.var], na.rm = T))*0.05*(nrow(table_stat[table_stat[,x] == dia,])-1), 
                                           by=diff(range(table[,value.var], na.rm = T))*0.05)
    t[table_stat[,x] == dia,"y2"]<-t[table_stat[,x] == dia,"y1"]
  }
  
  t_def<-as.data.frame(matrix(ncol=4, nrow=0))
  for (row in 1:nrow(t)){
    t_def<-rbind(t_def, t[row,],
                 c(t[row,"y1"]-diff(range(table[,value.var], na.rm = T))*bracket.length, t[row,"y1"], t[row,"x1"], t[row,"x1"]),
                 c(t[row,"y1"]-diff(range(table[,value.var], na.rm = T))*bracket.length, t[row,"y1"], t[row,"x2"], t[row,"x2"]))
  }
  t_lab<-data.frame("x"=t$x1+(t$x2-t$x1)/2, "y"=t$y1+diff(range(table[,value.var], na.rm = T))*0.005, check.names = F)
  return(list("label"=t_lab, "brackets"=t_def))
}

secfile<-function(file){
  ext<-strsplit(file, ".", fixed = T)[[1]]
  ext<-ext[length(ext)]
  num<-1
  while(file.exists(file) == T){
    if(num == 1){
      file_tmp<-strsplit(file, ".", fixed=T)[[1]]
      file<-paste0(paste(file_tmp[1:(length(file_tmp)-1)], collapse = "."),"_",num,".",file_tmp[length(file_tmp)])
    }else{
      file_tmp<-paste(strsplit(file, "_", fixed=T)[[1]][-length(strsplit(file, "_", fixed=T)[[1]])], collapse = "_")
      file<-paste0(file_tmp, "_", num,".",ext)
    }
    num<-num+1
  }
  return(file)
}