output$deseqContrastVarUI <- renderUI({
  validate(need(physeq(), ""))
  selectInput(
    "deseqContrastVar",
    label = "Experimental design : ",
    choices = c(sample_variables(physeq()))
  )
})

output$deseqContrastModUI <- renderUI({
  validate(need(physeq(), ""), need(input$deseqContrastVar, ""))
  checkboxGroupInput(
    "deseqContrastMod",
    label = "Contrast (exactly two required) : ",
    choices = NULL,
    inline = TRUE
  )
})

observe({
  validate(need(physeq(), ""), need(input$deseqContrastVar, ""))
  var <- levels(as.factor(get_variable(physeq(), input$deseqContrastVar)))
  updateCheckboxGroupInput(session,
                           inputId = "deseqContrastMod",
                           choices = var,
                           selected = var[c(1, 2)],
                           inline = TRUE
  )
})

output$deseqTitleUI <- renderUI({
  validate(need(physeq(), ""))
  textInput("deseqTitle",
            label = "Title : ",
            value = "Volcano Plot")
})

output$deseqPadjUI <- renderUI({
  validate(need(physeq(), ""))
  sliderInput("deseqPadj",
              label = "Adjusted p-value threshold :",
              min = 0,
              max = 1,
              value = 0.05)
})

output$deseqUI <- renderUI({
  validate(need(physeq(), ""))
  box(
    title = "Setting : " ,
    width = NULL,
    status = "primary",
    uiOutput("deseqContrastVarUI"),
    uiOutput("deseqContrastModUI"),
    uiOutput("deseqTitleUI"),
    uiOutput("deseqPadjUI")
  )
})

output$deseq <- metaRender2(renderPlot, {
  validate(
    need(physeq(), "Requires an abundance dataset"),
    need(length(input$deseqContrastMod) == 2, "Requires two conditions"))
  data <- physeq()
  
  metaExpr({
    deseq_data <- phyloseq_to_deseq2(data, ..(as.formula(paste("~", input$deseqContrastVar))))
    dds <- DESeq2::DESeq(deseq_data, sfType = "poscounts")
    
    results <- DESeq2::results(dds, 
                               tidy = TRUE, 
                               contrast = ..(c(input$deseqContrastVar, input$deseqContrastMod[1], input$deseqContrastMod[2]))) %>% 
      rename(OTU = row)
    
    da_volcano <- data.frame(
      #otu      = row.names(results),
      otu      = results$OTU,
      evidence = -log10(results$padj), 
      lfc      = results$log2FoldChange) %>% 
      na.omit()
    # add a threshol line
    y_axix_volcano_line <- -log10(..(input$deseqPadj))
    # Modify dataset to add new coloumn of colors
    da_volcano <- da_volcano %>%
      mutate(
        color = case_when(
          lfc > 0 & evidence > y_axix_volcano_line ~ "More", 
          lfc < 0 & evidence > y_axix_volcano_line ~ "Less", 
          TRUE                                     ~ "Equal"
        )
      )
    
    # Color corresponds to fold change directionality
    volcano_plot <- ggplot(da_volcano, 
                           aes(x = lfc, y = evidence, label = otu)) + 
      geom_point(aes(color = factor(color)), size = 1.75, alpha = 0.8, na.rm = T) + # add gene points
      geom_text() +
      theme_bw(base_size = 16) +                                                    # clean up theme
      theme(legend.position = "none") +                                             # remove legend 
      ggtitle(label = ..(input$deseqTitle)) +                                       # add title
      xlab(expression(log[2]("FoldChange"))) +                                      # x-axis label
      ylab(expression(-log[10]("adjusted p-value"))) +                              # y-axis label
      geom_vline(xintercept = 0, colour = "grey80", linetype = 2) +                                # add line at 0
      geom_hline(aes(yintercept = y_axix_volcano_line), yintercept = y_axix_volcano_line, colour = "grey80", linetype = 2) +
      annotate(geom = "text", 
               label = paste("padj =", ..(input$deseqPadj)), 
               x = min(da_volcano$lfc), 
               y = y_axix_volcano_line + 0.25, 
               size = 4,
               colour = "black",
               vjust = 0,
               hjust = 0) + # add pvalue threshold           
      scale_color_manual(values = c("More" = "red", "Less" = "chartreuse", "Equal" = "darkgray")) # change colors
    # Plot figure
    volcano_plot + scale_y_continuous(trans = "log1p")
  })
})

observeEvent(input$deseq_output_code,
             {
               displayCodeModal(
                 expandChain(
                   quote(library(phyloseq)),
                   quote(library(phyloseq.extended)),
                   quote(library(DESeq2)),
                   quote(library(ggplot2)),
                   quote(library(magrittr)),
                   quote(library(dplyr)),
                   
                   "# Replace `data` with you own data.",
                   output$deseq()
                 )
               )
             }
)