inst/shiny-app/server.r

In omicplotR: Visual Exploration of Omic Datasets Using a Shiny App

server <- function(input, output, session) {

################################################################################
    #downloads

    filtering_options <- reactive({
        #get the reactive inputs for downloadable file
        reactive_values <- paste("min.count <- ", input$mincounts, "\n",
        "min.prop <- ", input$minprop, "\n",
        "max.prop <- ", input$maxprop, "\n",
        "min.sum <- ", input$minsum, "\n",
        "min.reads <- ", input$minreads, "\n",
        "taxselect <- ", as.numeric(input$taxlevel), "\n",
        "taxoncheck <- ", as.numeric(input$taxlevel), "\n",
        "arrowcheck <- ", input$arrowcheckbox, "\n",
        "scale.slider <- ", input$scale, "\n",
        "removenames <- ", input$removesamplenames, "\n",
        "var.filt <- ", input$varslider, "\n",
        "abund <- ", input$abundcutoffbarplot, "\n",
        "dist <- ", input$dismethod, "\n",
        "clust <- ", input$clustermethod, "\n",
        "vals <- ", vals$data, "\n",
        "metaval <- ", metaval$data,
        sep = "")
    })

    #combine reactive values with script for PCA biplot

    metadata_script <- reactive({

        info <- sessionInfo()
        loaded <- loadedNamespaces()

        # from sessionInfo
        package <- grep("^package:", search(), value = TRUE)
        keep <- vapply(package, function(x) x == "package:base" ||
                    !is.null(attr(as.environment(x), "path")), NA)
        package <- .rmpkg(package[keep])

        # get rid of non-loaded packages (i.e., base)
        loaded <- loaded[!(loaded %in% package)]

        packages <- list()
        packages[loaded] <- loaded

        for (i in seq(loaded)) {

            packages[loaded[i]] <- paste(info$loadedOnly[[loaded[i]]]$Package, info$loadedOnly[[loaded[i]]]$Version, sep = "-")

        }

        packages.lines <- as.data.frame(unlist(packages))
        colnames(packages.lines) <- "loaded_packages"


        # get other packages
        other <- list()

        for (i in seq(info$otherPkgs)) {
            other[i] <- paste(info$otherPkgs[[i]]$Package, info$loadedOnly[[loaded[i]]]$Version, sep = "-")
        }

        other.packages <- as.data.frame(unlist(other))
        colnames(other.packages)[1] <- "other_packages"

        #return all the necessary info

        # generate output in the desired order
        output <- c("Session Info", "\n", info$R.version$version.string, info$platform, info$running,  "\n", "Loaded Packages", as.character(packages.lines$loaded_packages), "\n", "Other attached packages", as.character(other.packages$other_packages))

    })

    PCA_script <- reactive({
        y <- filtering_options()
        metadata <- metadata_script()
        PCA_script <- c(metadata, "\n", "Filtering options", y, readLines("./PCA_script.R"))
    })

    rab_script <- reactive({
        y <- filtering_options()
        metadata <- metadata_script()
        rab_script <- c(metadata, "\n", "Filtering options", y, readLines("./rab_script.R"))
    })

    effect_script <- reactive({
        y <- filtering_options()
        metadata <- metadata_script()
        effect_script <- c(metadata, "\n", "Filtering options", y, readLines("./effect_script.R"))
    })

    file_name <- reactive({
        inFile <- input$file1
        file_name <- inFile$name
    })

    #download script for PCA biplot
    output$PCA_download <- downloadHandler(
        filename=function(){paste(file_name(), "_PCA.r", sep = "")},
        content= function(file) {
            writeLines(PCA_script(), file)
        }
    )

    #download script for rab plots
    output$rab_download <- downloadHandler(
        filename=function(){paste(file_name(), "_rab.r", sep = "")},
        content= function(file) {
            writeLines(rab_script(), file)
        }
    )

    #download script for effect plots
    output$effect_download <- downloadHandler(
        filename=function(){paste(file_name(), "_effect.r", sep = "")},
        content= function(file) {
            writeLines(effect_script(), file)
        }
    )

    output$ebi_counts_download <- downloadHandler(
        filename = "ebi_counts.txt",
        content = function(file) {
            write.table(vals$ebi_counts, file, sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
        }
    )

################################################################################
  #observe events
  #quit button
  options(shiny.maxRequestSize=70*1024^2)
  observeEvent(input$stopApp, {
    stopApp(returnValue = invisible())
  })

  #reset filter
  observeEvent(input$reset, {
    vals$data <- NULL
    removeModal()
  })

  #pop up the modal
  observeEvent(input$choosemeta, {
    showModal(metaModal())
  })

  #convert input to reactive object
  observeEvent(input$update, {
    vals$data <-
    c(input$cnname, input$cnname2, input$cnname3, input$cnname4)
    metaval$data <- input$select_biplot_filter
    removeModal()
  })

  #pop up the modal for EBI input
    observeEvent(input$input_ebi_project, {
      showModal(ebiModal())
    })

    #convert input to reactive object
    observeEvent(input$update_EBI, {

        # make progress bar object
        progress <- shiny::Progress$new()
        # close on exit
        on.exit(progress$close())

        # set to zero
        progress$set(message = "Downloading dataset", value = 0)

        # get the MGnify project ID from the input
        sample_id <- input$ebi_id

        # base for MGnify api
        api_base = 'https://www.ebi.ac.uk/metagenomics/api/latest/'

        # define link to get samples
        samples_link <- paste0(api_base, "studies/", sample_id, "/samples")

        # extract number of pages to collect all samples by getting
        # number of pages from last page line
        last_link <- jsonlite::fromJSON(paste0(samples_link, "?page=1"))$links$last
        num_pages <- substr(last_link, nchar(last_link), nchar(last_link) + 1)

        # initiate all samples vector
        all_samples <- c()

        # update progress bar

        progress$inc(1/5, detail = "All samples found")

        # loop through all pages to get all sample names
        for (i in 1:num_pages) {
            current_samples <- jsonlite::fromJSON(paste0(samples_link, "?page=", i))$data$id
            all_samples <- c(all_samples, current_samples)
        }



        # for each sample name, get the list of runs
        runs_link <- paste0(api_base, "samples/")

        # loop through all samples to get all runs that are associated
        all_runs <- c()

        progress$inc(0, detail = "Downloading runs...")

        # this loops through the samples and gets all the run numbers. it takes a while for some reason (about 25 seconds for 70 samples)
        for (i in seq(all_samples)) {
            current_run <- jsonlite::fromJSON(paste0(runs_link, all_samples[i], "/runs"))$data$id
            all_runs <- c(all_runs, current_run)

            # print status message
        message(paste0("Finished getting run for sample ", all_samples[i]))
        progress$inc(0, detail = paste0("Downloading runs... ", i, "/", length(all_samples)))
        }

        progress$inc(1/5, detail = "All runs downloaded")

        # download the analysis accession for each of the runs
        go_slim_link <- paste0(api_base, "runs/")

        # get accession numbers for all runs
        all_acc <- c()

        progress$inc(0, detail = "Downloading accessions...")

        for (i in seq(all_runs)) {
            current_acc <- jsonlite::fromJSON(paste0(go_slim_link, all_runs[i], "/analyses"))$data$id
            all_acc <- c(all_acc, current_acc)
            # show message to give status
            message(paste0("Finished getting accession for run ", all_runs[i]))

            progress$inc(0, detail = paste0("Downloading accessions... ", i, "/", length(all_runs)))
        }

        progress$inc(1/5, detail = "All accessions found")

        # download counts for every accession
        analyses_link <- "https://www.ebi.ac.uk/metagenomics/api/v1/analyses/"

        # set up the counts table with the GO terms
        all_counts <- data.frame(GO = jsonlite::fromJSON(paste0(analyses_link, all_acc[i], "/go-slim"))$data$attributes$accession)

        all_counts$description <- jsonlite::fromJSON(paste0(analyses_link, all_acc[i], "/go-slim"))$data$attributes$description

        all_counts$category <- jsonlite::fromJSON(paste0(analyses_link, all_acc[i], "/go-slim"))$data$attributes$lineage

        # add columns for each run, name appropriately (+3 is because first three columns are descriptors
        all_counts[,4:(length(all_runs)+3)] <- "NA"

        # rename columns with run numbers
        colnames(all_counts)[4:length(all_counts)] <- all_runs

        progress$inc(0, detail = "Downloading count tables...")

        # add counts to the data table by column
        # this downloads go-slim counts for every sample.
        for (i in seq(all_acc)) {

            current_count <- jsonlite::fromJSON(paste0(analyses_link, all_acc[i], "/go-slim"))$data$attributes$count

            message(paste0("Finished accession ", all_acc[i]))

            progress$inc(0, detail = paste0("Downloading count tables... ", i, "/", length(all_acc)))

            # replace column column with counts for each accession, increments each iteration
            all_counts[,3+i] <- current_count
        }

        progress$inc(1/5, detail = "All count tables downloaded")

        vals$ebi_counts <- all_counts

        # count table generated
        progress$inc(0, detail = "Merged count table generated")

        # remove modal
        removeModal()

        # give pop up for saving file
        vals$download_finished <- TRUE
    })

    observeEvent(input$update_EBI, {
        if (vals$download_finished == TRUE) {
            showModal(ebi_download_modal())
        }
    })

    observeEvent(input$select_ef_columns, {
        showModal(select_columns())
    })

  #dendrogram brush ranges
  ranges <- reactiveValues(x = NULL, y = NULL)
  bp_ranges <- reactiveValues(x = NULL, y = NULL)

  ###########################################################################
  #custom UIs
  #make the pop up
  metaModal <- function (x) {
    modalDialog(
      renderUI({
        metadata <- metadata()
        #message if there is nothing uploaded
        if (is.null(metadata)) {
          selectInput(
            "select_biplot_filter",
            label = h3("Filter biplot by metadata column"),
            choices = list(
              "Input data" = 1
            ),
            selected = 1
          )
        } else {
          options <- colnames(metadata)
          selectInput(
            "select_biplot_filter",
            label = h3("Filter data by metadata column"),
            choices = options,
            selected = 1
          )
        }
      }),
      textInput("cnname", "Value 1", placeholder = "Enter value from metadata or leave blank"),
      textInput("cnname2", "Value 2", placeholder = "Enter value from metadata or leave blank"),
      textInput("cnname3", "Value 3", placeholder = "Enter value from metadata or leave blank"),
      textInput("cnname4", "Value 4", placeholder = "Enter value from metadata or leave blank"),
      title = "Choose metadata values to filter PCA biplot",
      "Enter metadata values of the metadata variable you have chosen (must be exact match).
      The biplot will be replotted with only the samples identified by
      the metadata values. Leave inputs blank if you don't need them.
      Samples that do not have metadata will be coloured black",
      footer = tagList(actionButton("update", "Update Filter"),
      actionButton("reset", "Reset Filter"),
      modalButton("Cancel")),
      easyClose = TRUE
    )
  }

  #pop up for inputting ebi project numberq
    ebiModal <- function (x) {
        modalDialog(
                textInput("ebi_id", "EBI Project ID", placeholder = "Project ID"),
                title = "Download EBI project to explore",
                "Enter the MGnify study ID (for example: MGYS00001110) to download the GO slim annotation for all samples. This will take a few minutes to gather all data depending on the size of the study. Internet access if required. Do not close window while download is in progress. File will be automatically inputted after downloading.",
              footer = tagList(actionButton("update_EBI", "Download and input"),
                modalButton("Cancel")),
                easyClose = TRUE
        )
    }

    ebi_download_modal <- function (x) {
        modalDialog(
            title = "Save counts table?", "Click yes to save the counts table locally. Input this file to avoid redownloading the data.",
            footer = tagList(h3(downloadButton("ebi_counts_download", "Yes")), modalButton("No"), modalButton("Close this window")), easyClose = TRUE
        )
    }

    select_columns <- function (x) {

        modalDialog(
            title = "Manually select conditions",
            "Enter the column numbers for each condition, separated by a comma. Condition 1 will be compared against condition 2. Scroll across table to see column names. Column numbers are shown below the column name.",
            DT::renderDataTable({
                data <- data.t()
                dataTable <- data.frame(matrix(0, ncol = length(colnames(data)), nrow = 1))
                colnames(dataTable) <- colnames(data)

                newrow <- 1:length(dataTable)
                dataTable <- rbind(newrow, dataTable)
                rownames(dataTable) <- c("Column number", "")
             dataTable

         }, options = list(pageLength = 5, scrollX = TRUE)),
            textInput("effect_cond_1", "Columns for condition 1", placeholder = "1, 2, 3, ..."),
            textInput("effect_cond_2", "Columns for condition 2", placeholder = "1, 2, 3, ..."),
            footer = tagList(actionButton("column_choice_ab", "Update conditions"), modalButton("Close window")), easyClose =TRUE
        )
    }

  output$taxchoice <- renderUI({
    if (input$taxoncheckbox) {
      radioButtons("taxlevel",
      "Choose taxonomic level to display",
      c("Kingdom" = 1,
      "Phylum" = 2,
      "Class" = 3,
      "Order" = 4,
      "Family" = 5,
      "Genus" = 6,
      "Species" = 7),
      selected = 4
    )
  }
})

#dynamic variance slider
output$varianceslider <- renderUI({
  data <- data()

  sliderInput(
    "varslider",
    "Variance cutoff",
    min = 0,
    max = 5,
    value = 0,
    step = 0.05
  )
})

#effect plot conditions
output$conditions<- renderUI({
  if (input$ep_chooseconds == 2) {
    meta <- metadata()
    options <- colnames(meta)
    cn <- input$colselect

    choice <- unique(meta[[cn]])
    tagList(
      selectInput("group1", "Choose group 1", choices = choice),
      selectInput("group2", "Choose group 2", choices = choice))
    }
  })

  #choose column from metadata
  output$colselectcond <- renderUI({
    meta <- metadata()
    options <- colnames(meta)
    selectInput("colselect", "Select column from metadata", choices = options)
  })

  ###########################################################################
  #data input
  #get data from uploaded file
  data <- reactive({
    if (input$exampledata) {
      read.table(
        "example_data.txt",
        header = TRUE,
        sep = "\t",
        stringsAsFactors = FALSE,
        quote = "",
        row.names = 1,
        check.names = FALSE,
        comment.char = "",
        na.strings = ""
      )
  } else if (input$exampledata2) {
        read.table(
          "selex.txt",
          header = TRUE,
          sep = "\t",
          stringsAsFactors = FALSE,
          quote = "",
          row.names = 1,
          check.names = FALSE,
          comment.char = "",
          na.strings = ""
      )} else if (!is.null(vals$ebi_counts)) {
            # load downloaded file as input
            data <- vals$ebi_counts

            # make rownames descriotion
            rownames(data) <- data$description

            # rearrange, put descriptors at end so it's easy to take off
            data <- cbind(data[,4:length(data)], GO = data[,1], category = data[,3])

        } else if (input$ebi_format == TRUE) {

            #reactive input file
            inFile <- input$file1

            #return NULL when no file is uploaded
            if (is.null(inFile))
            return(NULL)

            # input EBI formatted data
                data <- read.table(
                  inFile$datapath,
                  header = TRUE,
                  sep = "\t",
                  stringsAsFactors = FALSE,
                  quote = "",
                  row.names = NULL,
                  check.names = FALSE,
                  comment.char = "",
                  na.strings = "")

                  # make rownames descriotion
                  rownames(data) <- data$description

                  # rearrange, put descriptors at end so it's easy to take off
                  data <- cbind(data[,4:length(data)], GO = data[,1], category = data[,3])
        } else {
        #reactive input file
        inFile <- input$file1

        #return NULL when no file is uploaded
        if (is.null(inFile))
        return(NULL)

        #reads the file
        #if there's an error, it doesn't quit omicplotR
        data <- try(read.table(
          inFile$datapath,
          header = TRUE,
          sep = "\t",
          stringsAsFactors = FALSE,
          quote = "",
          row.names = 1,
          check.names = FALSE,
          comment.char = "",
          na.strings = ""), silent = TRUE)

        #if rownames error
        if(grepl("duplicate 'row.names'", data[1], fixed = TRUE)) {
          showModal(rownamesModal())
        } else {
          #if no error, then continue
          data <- data
        }
        }
}) #reactive

  #get metadata from uploaded file
  metadata <- reactive({

    if (input$exampledata) {
      read.table(
        "example_metadata.txt",
        header = TRUE,
        sep = "\t",
        stringsAsFactors = FALSE,
        quote = "",
        row.names = 1,
        check.names = FALSE,
        comment.char = ""
      )
    } else {
      #reactive input file
      inFile2 <- input$file2

      #return NULL when no file is uploaded
      if (is.null(inFile2))
      return(NULL)

      #reads the file
      read.table(
        inFile2$datapath,
        header = TRUE,
        sep = "\t",
        stringsAsFactors = FALSE,
        quote = "",
        row.names = 1,
        check.names = FALSE,
        comment.char = ""
      )
    }
  })

################################################################################
#tests to make sure input is correct format
formatModal <- function(failed = FALSE) {
    modalDialog(
        title="Incorrect format",
        "There are duplicate column names. Ensure all column names unique, then re-import your file.",
        easyClose = TRUE,
        footer = modalButton("Dismiss")
    )
}

  rownamesModal <- function(failed = FALSE) {
    modalDialog(
      title="Incorrect format",
      "There are duplicate row names. Ensure all row names unique, then re-import your file.",
      easyClose = TRUE,
      footer = modalButton("Dismiss")
    )
  }

#check data format
   observeEvent(input$showdata, {
     inFile <- data()

     #return NULL when no file is uploaded
     if (is.null(inFile)) {
         return(NULL)
     }

     #make frequency table of occurence of colnames
     c_occur <- data.frame(table(colnames(inFile)))

     #if frequency is more than one, show the "your format is wrong" modalDialog
     if (max(c_occur$Freq) > 1) {
         showModal(formatModal())
     } else {
         return(NULL)
     }
 })

#metadata check format
 observeEvent(input$showmetadata, {
   inFile2 <- metadata()

   #return NULL when no file is uploaded
   if (is.null(inFile2)) {
     return(NULL)
   }

   #make frequency table of occurence of colnames
   c_occur <- data.frame(table(colnames(inFile2)))

   #if frequency is more than one, show the "your format is wrong" modalDialog
   if (max(c_occur$Freq) > 1) {
     showModal(formatModal())
   } else {
     return(NULL)
   }
 })

################################################################################

  #input ALDEx2 table
  effect_input <- reactive({

    inFile3 <- input$effect_file

    #return NULL when no file is uploaded
    if (is.null(inFile3))
    return(NULL)

    read.table(
      inFile3$datapath,
      header = TRUE,
      sep = "\t",
      stringsAsFactors = FALSE,
      quote = "",
      row.names = 1,
      check.names = FALSE,
      comment.char = "",
      na.strings = ""
    )
  })

  #create reactive object for biplot filtering inputs
  vals <- reactiveValues(data = NULL)
  metaval <- reactiveValues(data = NULL)

  #get column selection for metadata
  column <- reactive({
    cn <- input$selectcolumn
  })

  ###########################################################################
  #data processing
  #generate taxonomy vector for biplots
  taxonomy <- reactive({
    data <- data.t()
    tax.result <- data.check()

    #if no taxonomy column, return NULL for later
    if (isTRUE(tax.result)) {
      tax <- NULL
    } else {
      #remove taxonomy column
      tax <- data$taxonomy
      data$taxonomy <- NULL

      #get genus names
      genus <- vapply(strsplit(as.character(tax), "[[:punct:]]"), "[", taxselected, FUN.VALUE=character(1))
    }
  })

  #if taxonomic column present, return boolean. no tax column returns TRUE
  data.check <- reactive ({
    data.r <- data()
    if (is.null(data.r$taxonomy)) {
      return(TRUE)
    } else {
      return(FALSE)
    }
  })

  #do transformations of data for biplot
  data.t <- reactive ({
    #get the data
    x <- data()
    meta <- metadata()
    cn <- column()


    #set to zero, but inputs from sliders/numeric inputs
    min.count <- input$mincounts
    min.prop <- input$minprop
    max.prop <- input$maxprop
    min.sum <- input$minsum
    min.reads <- input$minreads


    # validation statements
    validate((need((input$mincounts >=0), "Input minimum count per feature (0 or above)")))

   validate((need((input$minprop >=0), "Input minimum proportional abundance (0 or above)")))

   validate((need((input$maxprop <= 1), "Input maximum proportional abundance (1 or below)")))

   validate((need((input$minsum >= 0), "Input minimum count sum (0 or above)")))

   validate((need((input$minreads >= 0), "Input minimum count sum per samples (0 or above)")))

    #catch the errors
    validate((need(input$minprop, "Processing..")))

    #check for tax
    if (is.null(x$taxonomy)) {
      taxCheck <- TRUE

      if (input$ebi_format == TRUE) {
          # this removes the additional columns.
          order <- order(colnames(x[,1:(ncol(x) - 2) ]))
          test <- as.matrix(x[order])

      } else {
          test <- as.matrix(x[order(colnames(x))])
      }

     validate(need(min.reads < max(colSums(test)), "Minimum count sum per sample exceeds any sample. Reduce so not all samples are removed by filter."))

     # filter by min reads per sample
     test.0 <- test[, which(colSums(test) > min.reads)]

     validate(need(min.count < max(matrixStats::rowMaxs(test.0)), "Minimum count per feature exceeds maximum count. Reduce so not all features are removed by filter."))

     # filter by min count per feature
     test.1 <- test.0[which(matrixStats::rowMaxs(test.0) >= min.count), ]

     # filter by sum of count per feature
     test.2 <- test.1[which(rowSums(test.1) >= min.sum), ]

    } else {
      taxCheck <- FALSE

      test <- x[order(colnames(x[seq_len(length(x) - 1)]))]

      validate(need(min.reads < max(colSums(test)), "Minimum count sum per sample exceeds any sample. Reduce so not all samples are removed by filter."))

      # filter by min reads per sample
      test.0 <- as.matrix(test[, which(colSums(test[, seq_along(test)]) >=
            min.reads)])

      validate(need(length(which(matrixStats::rowMaxs(test.0) > min.count)) > 1, "Minimum count per feature exceeds maximum count. Reduce so not all features are removed by filter."))

      # filter by min count per feature
      test.1 <- test.0[which(matrixStats::rowMaxs(test.0) >= min.count),]

      # filter by sum of count per feature
      test.2 <- test.1[which(rowSums(test.1) >= min.sum), ]

    }

    #get filtered data if filtered
    if (is.null(vals$data)) {
      x <- x
    } else {

      x <- omicplotr.metadataFilter(x, meta, column = metaval$data, values = vals$data)
    }

    if (input$ebi_format == TRUE) {
        # remove last two columns (GO terms and descriptor columns)
        x.filt <- omicplotr.filter(x[,1:(ncol(x)-2)], min.reads = min.reads, min.count = min.count, min.prop = min.prop, max.prop = max.prop, min.sum = min.sum)

    } else {

    x.filt <- omicplotr.filter(x, min.reads = min.reads, min.count = min.count, min.prop = min.prop, max.prop = max.prop, min.sum = min.sum)
}
    return(x.filt)

  })

  #prcomp object
  data.prcomp <- reactive({
    #get data
    data.t <- data.t()

    # prevent further calculation if or less than 1 sample remains
    validate(need((dim(data.t)[1] > 0 & dim(data.t)[2] > 1), "Filtering has removed all features or samples. Reduce the stringency of the filters so you have features to plot."))

    var.filt <- input$varslider
    data <- data()

    validate(need(input$varslider, "Calculating..."))

    #check for tax
    if (is.null(data$taxonomy)) {
      taxCheck <- TRUE
    } else {
      taxCheck <- FALSE
    }

    #catch error if there is no data
    if (is.null(data)) {
      return(NULL)
    } else {
      if (isTRUE(taxCheck)) {
        data.t <- data.t
      } else {
        data.t$taxonomy <- NULL
      }


      # check fo variance filter
      if (any(data.t == 0)) {
        test.0 <- zCompositions::cmultRepl(t(data.t), label = 0, method = "CZM")
        }
      else {
        test.0 <- t(data.t)
        }
        test.clr <- as.matrix(log(test.0) - rowMeans(log(test.0)))

        # requires more than 1 feature to remain for plotting.
        validate(need(length(which(var.filt < matrixStats::colVars(test.clr))) > 1, "Variance filter exceeds maximum variance. Reduce the stringency of the filters so you have features to plot."))

      # set to CZM by default
      zeros <- "CZM"

      if (input$zero_replacement == 1) {
          zeros <- "CZM"
      }

      if (input$zero_replacement == 2) {
          zeros <- "pseudo"
      }

      data.pr <- omicplotr.clr(data.t, var.filt, zero.method=zeros)

    }
    return(data.pr)
  })

#these are for effect plot choices by metadata
cn <- ""
group1 <- ""
group2 <- ""
denom <- ""

 observeEvent(input$column_choice_ab, {
     effect_cond_2 <<- input$effect_cond_2
     effect_cond_1 <<- input$effect_cond_1

 })

  #computing aldex object
  d.clr <- reactive({

      # get reactive input from selecting conditions
      cn <- input$colselect
      group1 <- input$group1
      group2 <- input$group2
      denom <- input$denomchoice


    x <- data.t()

    meta <- metadata()

    #require user to click action button
    validate(need(input$effectplot_ab, ""))

    #get filtered metadata if it is filtered
    if (is.null(vals$data)) {
      x <- x
    } else {

      x <- omicplotr.metadataFilter(x, meta, column = metaval$data, values = vals$data)

      meta <- meta[which(rownames(meta) %in% colnames(x)),]
  }

    # Separate the taxonomy column from the counts
    if (is.null(x$taxonomy)) {
      d <- x
    } else {
      d <- x[, 0:(dim(x)[2] - 1)]
      taxon <- x[(dim(x)[2])]
    }

    if (input$ep_chooseconds == 1) {

        g1s <- effect_cond_1
        g2s <- effect_cond_2
        # rearrange the columns and make conditions
        # split on the comma and space to get all columns for each condtion
        group1 <- as.integer(unlist(strsplit(g1s, ", ")))
        group2 <- as.integer(unlist(strsplit(g2s, ", ")))

        # rearrange for aldex table
        d <- cbind(d[group1], d[group2])
        # make conditions length of each group

        # make conds a global variable that can be accessed by other functions
      conds <<- c(rep("group1", length(group1)),
      rep("group2", length(group2)))
    }

    if (input$ep_chooseconds ==2) {
      #filter the metadata and keep only the data which have been chosen
      group1 <- rownames(meta[which(meta[[cn]] == group1), ])
      group1.filt <- group1[group1 %in% (colnames(x))]
      data1.filt <- x[,which(colnames(x) %in% group1.filt)]
      one <- length(data1.filt)

      group2 <- rownames(meta[which(meta[[cn]] == group2), ])
      group2.filt <- group2[group2 %in% (colnames(x))]
      data2.filt <- x[,which(colnames(x) %in% group2.filt)]
      two <- length(data2.filt)

      # make conds a global variable that can be accessed by other functions
      conds <<- c(rep("Group 1", one),
      rep("Group 2", two))

      #combine
      d <- cbind(data1.filt, data2.filt)

    }

    # Generate ALDEx2 object and get a distribution of values that are centered
    # log ratio transformed at each feature*sample

    choices = c("all", "iqlr", "zero")

    method <- choices[as.numeric(denom)]
    withProgress(
      message = "Calculating",
      detail = "Calculating expected clr values", value = 1/2, {
        d.clr <- aldex.clr(d, mc.samples=128, conds = conds, denom = method, verbose=TRUE)
        incProgress(1/2, message = "clr values calculated")
      }
    )
    #output d.clr
    d.clr <- d.clr

  })

  aldex.obj <- reactive({
    x <- data.t()
    d.clr <- d.clr()
    meta <- metadata()

    #get filtered data if filtered
    if (is.null(vals$data)) {
      x <- x
    } else {

      x <- omicplotr.metadataFilter(x, meta, column = metaval$data, values = vals$data)

      meta <- meta[which(rownames(meta) %in% colnames(x)),]
  }

    if (is.null(x$taxonomy)) {
      d <- x
    } else {
      d <- x[, 0:(dim(x)[2] - 1)]
      taxon <- x[(dim(x)[2])]
    }

    withProgress(
      message = "Calculating t-test...",
      detail = "Calculating t test and effect size", value = 1/2, {
        # Run a t-test between the two groups specified by conds
        #TODO option to generate effect plots ttest to save time.
        # this will require changing aldex.plot in effect plot output
        x.tt <- aldex.ttest(d.clr, paired.test=FALSE)
        incProgress(1/4, message = "Calculating effect size...")
        # Generate effect sizes
        x.effect <- aldex.effect(d.clr, include.sample.summary=TRUE,
          verbose=TRUE)
          incProgress(1/4, message = "Efect sizes calculated")
        }

      )
      # Combine into one table
      x.all <- data.frame(x.tt, x.effect, stringsAsFactors=FALSE)
    })

    #stripchart data
    clr.strip <- reactive({
      x <- d.clr()
      codaSeq.repl <- function(x){
        # Initialize medians vector to be reshaped later to matrix
        MC.means <- numeric()


        # Iterating all samples
        for(i in seq_len(numConditions(x))) {
          sample.means <- numeric()

          # Get sample's MC Instances for features
          MC.matrix <- getMonteCarloInstance(x, i)

          # Iterating features
          for(j in seq_len(numFeatures(x))) {
            feature.means <- mean(MC.matrix[j, ])
            sample.means <- append(sample.means, feature.means)
          }

          # Add subject's MC instance medians to master vector
          MC.means <- append(MC.means, sample.means)
        }

        # Turn the vector of medians into a matrix and then data frame with
        result <- matrix(MC.means, nrow = numFeatures(x), ncol = numConditions(x))
        result <- data.frame(result, row.names = getFeatureNames(x))
        colnames(result) <- getSampleIDs(x)

        return(result)
      }

      clr.strip <- codaSeq.repl(x)
    })

    ###########################################################################
    # outputs
    #show removed OTUs/samples
    output$removedDT <- DT::renderDataTable({

      # catch error message
      validate(need(data() != "", ""),
          need(input$showremoved, ""))

          data.pr <- data.prcomp()
          data.in <- data()

      omicplotr.getRemovedSamples(data.in, data.pr)
    },
    #force datatable to size of window
    options = list(scrollX = TRUE)
  )

  output$removedDTotu <- DT::renderDataTable({


    validate(need(data() != "", ""),
              need(input$showremoved, "Click 'Show removed samples/OTUs' to view removed OTUs"))

    data.pr <- data.prcomp()
    data.in <- data()


    omicplotr.getRemovedFeatures(data.in, data.pr)
  },
  #force datatable to size of window
  options = list(scrollX = TRUE))

  #choose the metadata column you want to colour by
  output$choose_column <- renderUI({
    metadata <- metadata()

    #message if there is nothing uploaded
    if (is.null(metadata)) {
      selectInput(
        "selectcolumn",
        label = h3("Choose column to colour by"),
        choices = list(
          "No metadata file detected" = 1
        ),
        selected = 1
      )
    } else {
      options <- colnames(metadata)
      selectInput(
        "selectcolumn",
        label = h3("Choose column to colour by"),
        choices = options,
        selected = 1
      )
    }
  })

  output$textTitle <- renderText({
    "Choose filtering options"
  })

  output$nometadata <- renderText({
    metadata <- metadata()
    if (is.null(metadata)) {
      "Select metadata file"
    }
  })

  output$nodata <- renderText({
    metadata <- metadata()
    if (is.null(data)) {
      "Select data file"
    }
  })

  output$nostripchart <- renderText({
    validate(need(input$effectplot_ab, "Select conditions and click 'Generate effect plot'"))
  })

  output$test <- renderText({
    if (!is.null(vals$data)) {
      c("Data filtered by metadata column:", metaval$data,  "value:", vals$data)
    } else {
      return(NULL)
    }
  })

  output$filter_warning_effect <- renderText({
    if (!is.null(vals$data)) {
      c("Data filtered by metadata column:", metaval$data,  "value:", vals$data)
    } else {
      return(NULL)
    }
  })

    output$filter_warning_dendro <- renderText({
      if (!is.null(vals$data)) {
        c("Data filtered by metadata column:", metaval$data,  "value:", vals$data)
      } else {
        return(NULL)
      }
    })

  output$datatable <- DT::renderDataTable({
    validate((need(input$showdata, "Click 'Check data' to check format and display table")))

    if (input$showdata) {
      data <- data()

    }
  }, options = list(scrollX = TRUE)) #force data table to window size

  output$metadatatable <- DT::renderDataTable({

    validate((need(input$showmetadata, "Click 'Check metadata' to check format and display table")))

    if (input$showmetadata) {
      meta <- metadata()
    }
  }, options = list(scrollX = TRUE)) #force data table to window size

  output$contents <- renderTable(include.rownames = TRUE, {
    data()
  })

  ##############################################################################

  # PCA biplots

  #biplot
  output$biplot <- renderPlot({

      validate(need(data() !="", "Input a correctly formatted data file"))

    #get reactive objects
    data <- data.prcomp()
    tax <- data.t()
    removenames <- input$removesamplenames
    scale.slider <- input$scale
    arrowcheck <- input$arrowcheckbox
    taxoncheck <- input$taxoncheckbox
    taxselect <- as.numeric(input$taxlevel)
    #title <- input$biplot_title
    opacity <- input$opacity_samples_pca
    sample_size <- input$size_samples_pca

    x.var <- sum(data$sdev ^ 2)
    PC1 <- paste("PC 1 Variance: %", round(sum(data$sdev[1] ^ 2) / x.var * 100, 1))
    PC2 <- paste("PC 2 Variance: %", round(sum(data$sdev[2] ^ 2) / x.var*100, 1))

    #if no data don't display anything
    if (is.null(data))
    return(NULL)

    if (input$taxoncheckbox) {
      validate(need(input$taxlevel, "Processing..."))
    }

    taxonomy <- tax$taxonomy

    #if Show taxonomy checkbox is clicked, take chosen level
    if (isTRUE(taxoncheck)) {
        taxselect <- as.numeric(input$taxlevel)
    } else {
        taxselect <- 6
    }

    #get genus (or other level)
    genus <- vapply(strsplit(as.character(taxonomy), "[[:punct:]]"), "[", taxselect, FUN.VALUE=character(1))

    #biplot colouring options
    col = c(rgb(0,0,0,opacity), rgb(0, 0, 0, 0.2))

    #do checks for arrows
    if (isTRUE(arrowcheck)) {
      arrows = FALSE
    } else {
      arrows = TRUE
    }

    #if taxonomy is is null, use periods as points for features..
    #if taxoncheckbox is checked, show genus instead of points, otherwise
      if (isTRUE(taxoncheck)) {
        points <- genus
      } else {
        points <- c(rep(".", length(dimnames(data$rotation)[[1]])))
      }

    #remove sample names
    if (isTRUE(removenames)) {
      xlabs <- c(rep(".", length(dimnames(data$x)[[1]])))
      size = c(5, 0.8)
    } else {
      xlabs = unlist(dimnames(data$x)[1])
      size = c(sample_size, 0.8)
    }

    biplot(
      data,
      main = "Principal Component Analysis",
      cex.main = 1.5,
      cex = size,
      col = col,
      scale = scale.slider,
      var.axes = arrows,
      xlab = PC1,
      ylab = PC2,
      xlabs = xlabs,
      ylabs = points
    )
  })

  #coloured biplot
  output$coloredBiplot <- renderPlot({

      #file check
      validate(need(data() !="", "Input a correctly formatted data file"))

    #get reactive objects
    d <- data()
    data.prcomp <- data.prcomp()
    tax <- data.t()
    tax.check <- data.check()
    meta <- metadata()
    cn <- column()
    taxoncheck <- input$taxoncheckbox
    taxselect <- as.numeric(input$taxlevel)
    opacity <- input$opacity_samples_pca
    sample_size <- input$size_samples_pca

    validate(need(input$selectcolumn, "Click 'Colouring Options' to choose metadata for colouring biplot"))

    if (is.null(vals$data)) {
      #reactive metadata instead
      df <- as.data.frame(unlist(dimnames(data.prcomp$x)[1]))

    } else {
      #order the metadata
      meta <- meta[order(rownames(meta)), ]

      #get the input from the pop up
      select <- vals$data

      #get rid of blank values
      s <- select[select != ""]

      #filter metatable if it was filtered
      meta <- meta[which(meta[[metaval$data]] %in% s),]

      df <- unlist(dimnames(data.prcomp$x)[1])

      df <- as.data.frame(df[which(df %in% rownames(meta))])

    }

    if (input$taxoncheckbox) {
      validate(need(input$taxlevel, "Processing..."))
    }

    #if no data don't display anything
    if (is.null(meta))
    return(NULL)

    #get rid of taxonomy column if its there
    if (!isTRUE(tax.check)) {
      d$taxonomy <- NULL
    }

    colourvector <- omicplotr.colvec(data.prcomp, meta, opacity, cn, type = input$colouringtype)

    omicplotr.colouredPCA(data.prcomp, colourvector, scale = input$scale, arrows = input$arrowcheckbox, taxonomy = tax$taxonomy, show.taxonomy = taxoncheck, tax.level = taxselect, removenames = input$removesamplenames, names.cex=sample_size)
  })

  #histograms
  output$metahist <- renderPlot({

      # file check
      validate(need(data() !="", ""))

    #import reactive objects
    cn <- column()
    meta <- metadata()
    data <- data.prcomp()

    validate(need(input$selectcolumn, ""))

    if (is.null(meta))
    return(NULL)

    #get sample names
    df <- unlist(dimnames(data$x)[1])

    #filter rownames meta based on sample names in data
    meta <- meta[rownames(meta) %in% df,]

    #filter meta based on colnames data so it updates

    #colouring options
    if (input$colouringtype == 1) {

      validate(need(is.numeric(meta[[cn]]), "Nonnumeric data detected. Please click 'Nonnumeric metadata'"))

      #get length of unique values
      #dont laugh at my variables
      un <- unique(meta[[cn]])
      uniq <- un[!is.na(un)]
      uni <- length(uniq)

      #ramp it up
      colfunc <- colorRampPalette(c("red", "blue"))

      #make the colour
      col <- colfunc(uni)

      barplot(
        table(meta[[cn]]),
        space = 0,
        col = col,
        main = "Variable frequency",
        xlab = cn
      )
    } else if (input$colouringtype == 2) {

      validate(need(is.numeric(meta[[cn]]), "Nonnumeric data detected. Please click 'Nonnumeric metadata'"))

      #calculate quartile
      q <- quantile(meta[[cn]])

      #get freq/value
      t <- as.data.frame(table(meta[[cn]]))

      #convert factors to numeric
      t[, 1] <- as.numeric(as.character(t[, 1]))

      #add quartile rank column
      tn <-
      within(t, quartile <-
        as.integer(cut(
          t[, 1], quantile(meta[[cn]]), include.lowest = TRUE
        )))

        c <- colorRampPalette(c("red", "blue"))(4)

        #replace with colours
        tn$quartile[tn$quartile == 1] <- c[1]
        tn$quartile[tn$quartile == 2] <- c[2]
        tn$quartile[tn$quartile == 3] <- c[3]
        tn$quartile[tn$quartile == 4] <- c[4]

        barplot(
          table(meta[[cn]]),
          col = tn$quartile,
          space = 0,
          main = "Variable frequency",
          xlab = paste(cn, "by Quartiles")
        )

      } else if (input$colouringtype == 3) {
        #get frequency table
        t <- as.data.frame(table(meta[[cn]], exclude = NULL))

        #make colour column
        t$col <- NA

        #get unique features of columns in metadata
        u <- as.character(unique(meta[[cn]]))

        #colour palette
        colours <- c("indianred1", "steelblue3",  "skyblue1", "mediumorchid",
        "royalblue4", "olivedrab3", "pink", "#FFED6F", "mediumorchid3",
        "ivory2", "tan1", "aquamarine3", "#C0C0C0", "mediumvioletred",
        "#999933", "#666699", "#CC9933", "#006666", "#3399FF",
        "#993300", "#CCCC99", "#666666", "#FFCC66", "#6699CC",
        "#663366", "#9999CC", "#CCCCCC", "#669999", "#CCCC66",
        "#CC6600", "#9999FF", "#0066CC", "#99CCCC", "#999999",
        "#FFCC00", "#009999", "#FF9900", "#999966", "#66CCCC",
        "#339966", "#CCCC33", "#EDEDED"
      )

      #the tricky part here is to order the colouring based on what the coloredBiplot is doing,
      #so you need to order t based on the order of the unique values found in meta
      t <- t[match(u, t$Var1), ]

      #colour NAs black or follow color scheme
      if (NA %in% t$Var1) {
        y <- which(is.na(t$Var1))
        t$col[y] <- "black"
        for (i in (seq_len(length(t$Var1)))) {
          if (is.na(t$Var1[i])) {
            next
          } else {
            t$col[i] <- colours[i - 1]
          }
        }

        #colour NAs black
      } else {
        for (i in (seq_len(length(t$Var1)))) {
          t$col[i] <- colours[i]
        }
      }

      barplot(
        t$Freq,
        space = 0,
        col = t$col,
        main = "Number of occurences",
        names.arg = as.character(t$Var1),
        xlab = cn,
        ylab = "Frequency"
      )
    }
  })

  #generate screeplot
  output$screeplot <- renderPlot({

      # file check
      validate(need(data() !="", ""))

    data <- data.prcomp()

    #if no data don't display anything
    if (is.null(data.t()))
    return(NULL)

    #generate screeplot
    screeplot(data, main = "Screeplot of Variances", xlab = "Component")
  })

  #removed samples

  #colsums for samples
  output$colsums <- renderPlot({
    data <- data()

    if (is.null(data)) {
      return(NULL)
    }

    ab <- input$minreads

    if (is.null(data$taxonomy)){
      x <- colSums(data)

      plot(x, ylab = "Counts", xlab = "Sample Number", pch = 19, col = ifelse({x > ab}, "gray0", "red"), main = "Samples removed by filtering (count sum)")

      abline(h = ab, col = "red", lty = 2, lwd = 2)
      legend("topright", legend = c("Remaining", "Removed"), col = c("black", "red"), pch = 19)

    } else {
      x <- colSums(data[,(seq_along(data) - 1)])

      plot(x, ylab = "Counts", xlab = "Sample Number", pch = 19, col = ifelse({x > ab}, "gray0", "red"), main = "Samples removed by filtering")

      abline(h = ab, col = "red", lty = 2, lwd = 2)
      legend("topright", legend = c("Remaining", "Removed"), col = c("black", "red"), pch = 19)
    }
  })

  output$rowsums <- renderPlot({
    data <- data()

    if (is.null(data)) {
      return(NULL)
    }

    ab <- input$minsum

    if (is.null(data$taxonomy)){
      x <- rowSums(data)

      plot(x, ylab = "Counts", xlab = "Row Number", pch = 19, col = ifelse({x > ab}, "gray0", "grey"), main = "Rows removed by filtering (count sum)")

      abline(h = ab, col = "grey", lty = 2, lwd = 2)
      legend("topright", legend = c("Remaining", "Removed"), col = c("black", "grey"), pch = 19)

    } else {
      x <- rowSums(data[,(seq_along(data) - 1)])

      plot(x, ylab = "Counts", xlab = "Row Number", pch = 19, col = ifelse({x > ab}, "gray0", "grey"), main = "Rows removed by filtering")

      abline(h = ab, col = "grey", lty = 2, lwd = 2)
      legend("topright", legend = c("Remaining", "Removed"), col = c("black", "grey"), pch = 19)
    }
  })

##############################################################################

goslim_stripchart <- reactive({

        # sanity checks. user needs data and EBI formatted data.
        validate(
            need(data() != "", "You need data. Input your data."),
            need(input$ebi_format == TRUE, "This page requires a GO Slim formatted count table. Input the correctly formatted data, and check the GO Slim format box to generate a plot."),
            need(!(is.null(input$varslider)), "Principal component object needs to be calculated. Please view Colored PCA Biplot before GO slim chart.")
        )

         # get input date
         data <- data()
         data.prcomp <- data.prcomp()
         meta <- metadata()
         opacity <- input$opacity_samples_pca
         cn <- column()

        # rearrange the input to format for strip chart function
        data <- cbind(GO = data$GO, category = data$category, data[,1:(ncol(data) - 2)])

        # if metadata is absent, don't colour by sample.

        # metadata present

        if (!is.null(meta)) {

            validate(need(!is.null(cn), "Please select columns to colour by under the Options panel from Colored Biplot."))
            colourvector <- omicplotr.colvec(data.prcomp, meta, opacity, cn, type = input$colouringtype)

        } else {

            # no metadata
            # make colourvector for number of samples
            colourvector <- colorRampPalette(c("black", "white"))(length(data) - 2)
        }

         createStripcharts <- function(data) {

             #create dataframe with 0.5 added so there are no zeros for log function later

             data.n0 <- cbind(data[,1:2], data[,3:ncol(data)] + 0.5)

             #create dataframes with information needed (calculated percent reads)
             molFunSet <- subset(data.n0, category == "molecular_function")
             logMFSet <- cbind(molFunSet[,1:2], log(molFunSet[,3:ncol(molFunSet)]))
             finalMFSet <- (logMFSet[,3:ncol(logMFSet)]) - colMeans(logMFSet[,3:ncol(logMFSet)])
             #rownames(finalMFSet) <- molFunSet[,1]

             bioProcSet <- subset(data.n0, category == "biological_process")
             logBPSet <- cbind(bioProcSet[,1:2], log(bioProcSet[,3:ncol(bioProcSet)]))
             finalBPSet <- (logBPSet[,3:ncol(logBPSet)]) - colMeans(logBPSet[,3:ncol(logBPSet)])
             #rownames(finalBPSet) <- bioProcSet[,1]

             cellCompSet <- subset(data.n0, category == "cellular_component")
             logCCSet <- cbind(cellCompSet[,1:2], log(cellCompSet[,3:ncol(cellCompSet)]))
             finalCCSet <- (logCCSet[,3:ncol(logCCSet)]) - colMeans(logCCSet[,3:ncol(logCCSet)])
             #rownames(finalCCSet) <- cellCompSet[,1]

             #create room for description titles
             par(mar=c(5.1, 23.1, 4.1, 1.1))
             par(mfrow=c(1,3))

             #loop to create and colour stripcharts
             for (i in 1:length(finalMFSet)) {
                 if (i == 1) {
                     stripchart(finalMFSet[,i] ~ rownames(finalMFSet), method = "jitter", jitter = 0.2, pch = 19, las = 2, cex = 0.7, cex.axis = 0.8, group.names = rownames(finalMFSet), xlab = "Centered Log of Percent Reads", col = colourvector[i], main = "Molecular Function")
                 }
                 else {
                     stripchart(finalMFSet[,i] ~ rownames(finalMFSet), method = "jitter", jitter = 0.2, pch = 19, las = 2, cex = 0.7, add = TRUE, col = colourvector[i])
                 }
             }

             #loop to divide description from description and ease readability
             for (j in 0.5:(length(rownames(finalMFSet))+0.5)){
                 abline(h=j, lty=3, col="grey80")
             }

             #add vertical line to mark x = 0
             abline(v=0, lty=3, col = "black")

             #loop to create and colour stripcharts
             for (i in 1:length(finalBPSet)) {
                 if (i == 1) {
                     stripchart(finalBPSet[,i] ~ rownames(finalBPSet), method = "jitter", jitter = 0.2, pch = 19, las = 2, cex = 0.7, cex.axis = 0.8, group.names = rownames(finalBPSet), xlab = "Centered Log of Percent Reads", col = colourvector[i], main = "Biological Process")
                 }
                 else {
                     stripchart(finalBPSet[,i] ~ rownames(finalBPSet), method = "jitter", jitter = 0.2, pch = 19, las = 2, cex = 0.7, add = TRUE, col = colourvector[i])
                 }
             }

             #loop to divide description from description and ease readability
             for (j in 0.5:(length(rownames(finalBPSet))+0.5)){
                 abline(h=j, lty=3, col="grey80")
             }

             #add vertical line to mark x = 0
             abline(v=0, lty=3, col = "black")

             #loop to create and colour stripcharts
             for (i in 1:length(finalCCSet)) {
                 if (i == 1) {
                     stripchart(finalCCSet[,i] ~ rownames(finalCCSet), method = "jitter", jitter = 0.2, pch = 19, las = 2, cex = 0.7, cex.axis = 0.8, group.names = rownames(finalCCSet), xlab = "Centered log ratio", col = colourvector[i], main = "Cellular Component")
                 }
                 else {
                     stripchart(finalCCSet[,i] ~ rownames(finalCCSet), method = "jitter", jitter = 0.2, pch = 19, las = 2, cex = 0.7, add = TRUE, col = colourvector[i])
                 }
             }

             #loop to divide description from description and ease readability
             for (j in 0.5:(length(rownames(finalCCSet))+0.5)){
                 abline(h=j, lty=3, col="grey80")
             }

             #add vertical line to mark x = 0
             abline(v=0, lty=3, col = "black")
         }

         if (input$returnpdf) {
             pdf("goslim_stripchart.pdf", width = 25, height = nrow(data)/6.8)
            createStripcharts(data)
            dev.off()
         }
         createStripcharts(data)
     })

# this dynamically sets the height to adjust to the size of the window.
output$ebi_stripchart <- renderPlot(expr = {goslim_stripchart()},

# dynamically gets height of plotting window.
height = function() {
    session$clientData$output_ebi_stripchart_width
}

)
output$pdflink <- downloadHandler(
    filename <- "ebi_stripchart.pdf",
    content <- function(file) {
        file.copy("goslim_stripchart.pdf", file)
    }
)
output$pdf_note <- renderText({"Note: it may take several seconds to generate the PDF before download is available."})


  ##############################################################################
  #this will be implemented in a later version

  # #association plot
  #
  # output$associationplot <- renderPlot({
  #   d <- data()
  #   rhocut <- input$rhocutoff
  #
  #   if (is.null(d)) {
  #     return(NULL)
  #   }
  #   associationPlot <- function(x, cutoff = rhocut) {
  #     # Separate out the taxonomy column from the counts
  #     d <- x[, 0:(dim(x)[2] - 1)]
  #     taxon <- x[(dim(x)[2])]
  #
  #     # Generate ALDEx2 object so values are centered log ratio transformed,
  #     # then convert to propr object with a symmetric rho statistic matrix
  #     d.clr <- aldex.clr(d, mc.samples = 128, verbose = TRUE)
  #     d.sma.df <- aldex2propr(d.clr, how = "perb")
  #
  #     # Make a reference list with pairs related by a rho statistic
  #     # less than the cutoff
  #     d.sma.lo.rho <- d.sma.df["<", cutoff]
  #
  #     # **** propr and ALDEx stuff done ****
  #
  #     # igraph: Convert the connections into a graphical object using propr's
  #     # cytescape function to first generate a table of indexed pairs and
  #     # proportions
  #     g <- graph.data.frame(cytescape(d.sma.lo.rho), directed = FALSE)
  #
  #     # igraph: Find the clusters
  #     g.clust <- clusters(g)
  #
  #     # Make a table to examine the cluster membership by hand
  #     g.df <-
  #     data.frame(
  #       Systematic.name = V(g)$name,
  #       cluster = g.clust$membership,
  #       cluster.size = g.clust$csize[g.clust$membership]
  #     )
  #
  #     # Generate a set of clusters larger than some size
  #     # Minimum cluster size is 2 (obviously)
  #     big <- g.df[which(g.df$cluster.size >= 2), ]
  #     colnames(big) <- colnames(g.df)
  #
  #     # Get genera
  #     genera <- c()
  #     for (i in 1:dim(taxon)[1]) {
  #       genera <- c(genera, sapply(strsplit(as.character(taxon[i,]), "[[:punct:]]"),
  #       "[", 6))
  #     }
  #
  #     # igraph: Rename the cluster members by their genus name
  #     # gsub(pattern, replacement, strings, perl-syntax)
  #     V(g)$name <- gsub("(^[A-Za-z]{3}).+", "\\1",
  #     as.vector(genera[V(g)]), perl = TRUE)
  #
  #     plot.new()
  #
  #     # igraph:
  #     # vertex.size controls point and text color
  #     # vertex.color controls point color
  #     # vertex.frame controls point outline color
  #     plot(
  #       g,
  #       vertex.size = 5,
  #       vertex.color = rgb(0, 0, 0, 0.2),
  #       vertex.frame.color = "white"
  #     )
  #   }
  #   associationPlot(d)
  # })
  #
  # output$associationtext <- renderText({
  #   "This association plot uses Spearman's rho to measure the strength and direction of association between two variables.
  #   The cutoff values ranges from -1 to 1. If greater than 1000 features are detected, input will be required on the R console."
  # })

  ################################################################################

  # Relative abundance plots

  #dendrogram
  output$dendrogram <- renderPlot({

      validate(need(data.check() == FALSE, "These plots require a count table with the final column titled taxonomy. Ensure the input is correctly formatted."))

    x <- data.t()
    meta <- metadata()
    abund <- input$abundcutoffbarplot

    observeEvent(input$dendro_dblclick, {
      brush <- input$dendro_brush
      if (!is.null(brush)) {
        ranges$x <- c(brush$xmin, brush$xmax)
        ranges$y <- c(brush$ymin, brush$ymax)
      } else {
        ranges$x <- NULL
        ranges$y <- NULL
      }
    })

    #get filtered data if filtered
    if (is.null(vals$data)) {
      x <- x
    } else {

      x <- omicplotr.metadataFilter(x, meta, column = metaval$data, values = vals$data)
    }

    if (is.null(x$taxonomy)) {
      d <- x
    } else {
      d <- x[, seq_along(x) - 1]
      taxon <- x[(dim(x)[2])]
    }

    validate(need(x$taxonomy, "No taxonomy column detected"))

    # Get genera
    genera <- c()
    for(i in (seq_len(nrow(taxon)))) {
      genera <- c(genera, vapply(strsplit(as.character(taxon[i, ]), "[[:punct:]]"),
      "[", 6, FUN.VALUE=character(1)))
    }

    # sum counts by name
    d.agg <- aggregate(d, by=list(genera), FUN=sum)
    tax.agg <- d.agg$Group.1
    d.agg$Group.1 <- NULL

    # convert to abundances
    d.prop <- apply(d.agg, 2, function(x){x/sum(x)})

    #filters by abundance (slider bar)
    d.abund <- d.agg[apply(d.prop, 1, max) > abund,]
    tax.abund.u <- tax.agg[apply(d.prop, 1, max) > abund]

    if (any(d.abund == 0)) {
    d.abund <- t(zCompositions::cmultRepl(t(d.abund), label = 0, method = "CZM"))
    } else {
        d.abund <- d.abund
    }
    # get proportions of the filtered data for plotting below
    # in log-ratio speak, you are re-closing your dataset
    d.P.u <- apply(d.abund, 2, function(x){x/sum(x)})

    # order by OTU abundances
    new.order <- rownames(d.P.u)[order(apply(d.P.u, 1, sum), decreasing=TRUE)]
    tax.abund <- tax.abund.u[order(apply(d.P.u, 1, sum), decreasing=TRUE)]
    d.P <- d.P.u[new.order, ]
    d.clr <- apply(d.P, 2, function(x){log2(x) - mean(log2(x))})

    #distance matrix
    inp <- as.numeric(input$dismethod)
    method.m <- c("euclidean", "maximum", "manhattan")
    method.mchoice <- method.m[inp]

    dist.d.clr <- dist(t(d.clr), method=method.mchoice)

    #get clustering method for dendrogram
    num <- as.numeric(input$clustermethod)
    method <- c("complete", "single", "ward.D2")
    method.choice <- method[num]

    clust.d <- hclust(dist.d.clr, method=method.choice)

    plot.new()
    par(fig=c(0, 1, 0, 1), new=TRUE)

    plot(as.dendrogram(clust.d), main=NULL, cex=0.8, xlab="", xlim = ranges$x, ylim= ranges$y, xpd = TRUE)
  })

  output$dendrotext <- renderText({
    "Drag mouse and double click to zoom in. \nDouble click again to zoom out"
  })

  #taoxnomic distribution
  output$barplot <- renderPlot({

      #no need for message to be repeated, but still squashes error.
    validate(need(data.check() == FALSE, ""))

    x <- data.t()
    meta <- metadata()
    abund <- input$abundcutoffbarplot

    observeEvent(input$bp_dblclick, {
      bp_brush <- input$bp_brush
      if (!is.null(bp_brush)) {
        bp_ranges$x <- c(bp_brush$xmin, bp_brush$xmax)
        bp_ranges$y <- c(bp_brush$ymin, bp_brush$ymax)
      } else {
        bp_ranges$x <- NULL
        bp_ranges$y <- NULL
      }
    })

    #get filtered data if filtered
    if (is.null(vals$data)) {
      x <- x
    } else {

      x <- omicplotr.metadataFilter(x, meta, column = metaval$data, values = vals$data)
    }

    if (is.null(x$taxonomy)) {
      d <- x
    } else {
      d <- x[, 0:(dim(x)[2] - 1)]
      taxon <- x[(dim(x)[2])]
    }

    validate(need(x$taxonomy, ""))

    # Get genera
    genera <- c()
    for(i in (seq_len(nrow(taxon)))) {
      genera <- c(genera, vapply(strsplit(as.character(taxon[i, ]), "[[:punct:]]"),
      "[", 6, FUN.VALUE=character(1)))
    }

    # sum counts by name
    d.agg <- aggregate(d, by=list(genera), FUN=sum)
    tax.agg <- d.agg$Group.1
    d.agg$Group.1 <- NULL

    # convert to abundances
    d.prop <- apply(d.agg, 2, function(x){x/sum(x)})

    d.abund <- d.agg[apply(d.prop, 1, max) > abund,]
    tax.abund.u <- tax.agg[apply(d.prop, 1, max) > abund]

    if (any(d.abund == 0)) {
    d.abund <- t(zCompositions::cmultRepl(t(d.abund), label = 0, method = "CZM"))
    } else {
        d.abund <- d.abund
    }

    # get proportions of the filtered data for plotting below
    # in log-ratio speak, you are re-closing your dataset
    d.P.u <- apply(d.abund, 2, function(x){x/sum(x)})

    # order by OTU abundances
    new.order <- rownames(d.P.u)[order(apply(d.P.u, 1, sum), decreasing=TRUE)]
    tax.abund <- tax.abund.u[order(apply(d.P.u, 1, sum), decreasing=TRUE)]
    d.P <- d.P.u[new.order, ]
    d.clr <- apply(d.P, 2, function(x){log2(x) - mean(log2(x))})

    #distance matrix
    inp <- as.numeric(input$dismethod)
    method.m <- c("euclidean", "maximum", "manhattan")
    method.mchoice <- method.m[inp]

    dist.d.clr <- dist(t(d.clr), method=method.mchoice)

    #clustering
    num <- as.numeric(input$clustermethod)
    method <- c("complete", "single", "ward.D2")
    method.choice <- method[num]
    clust.d <- hclust(dist.d.clr, method=method.choice)

    # standard colour scheme (Jean Macklaim)
    colours <- c("steelblue3", "skyblue1", "indianred", "mediumpurple1",
    "olivedrab3", "pink", "#FFED6F", "mediumorchid3",
    "ivory2", "tan1", "aquamarine3", "#C0C0C0", "royalblue4",
    "mediumvioletred", "#999933", "#666699", "#CC9933", "#006666",
    "#3399FF", "#993300", "#CCCC99", "#666666", "#FFCC66",
    "#6699CC", "#663366", "#9999CC", "#CCCCCC", "#669999",
    "#CCCC66", "#CC6600", "#9999FF", "#0066CC", "#99CCCC",
    "#999999", "#FFCC00", "#009999", "#FF9900", "#999966",
    "#66CCCC", "#339966", "#CCCC33", "#EDEDED")

    plot.new()
    par(fig=c(0,1,0,1), new = TRUE)

    par(fig=c(0,0.80,0, 1), new = TRUE)

    barplot(d.P[,clust.d$order], names.arg = clust.d$labels, space=0, xlim = bp_ranges$x, ylim= bp_ranges$y, col=colours, las=2, axisnames=TRUE, border = NA, xpd = TRUE)

    par(fig=c(0.8,1, 0, 1), new=TRUE)
    par(xpd = TRUE)

    leg.col <- rev(colours[seq_len(nrow(d.P))])

    legend(x="center", legend=rev(tax.abund), col=leg.col, lwd=5, cex=0.8,
    border=NULL)
  })



################################################################################

# effect plots

#find points by string and colour them
point.colour <- eventReactive(input$update_points, {
  effect <- effect_input()
  point.colour <- input$point.colour

  row.num <- grep(point.colour, rownames(effect))

  points(effect$diff.win[row.num], effect$diff.btw[row.num], pch=19, col=rgb(1,0,0,0.8), cex=1)
})

#find points by string and colour them
BA.point.colour <- eventReactive(input$update_points, {
  effect <- effect_input()
  point.colour <- input$point.colour

  row.num <- grep(point.colour, rownames(effect))

  points(effect$rab.all[row.num], effect$diff.btw[row.num], pch=19, col=rgb(1,0,0,0.8), cex=1)
})

#effect plots after calculating
output$effectMW <- renderPlot({

    # ensure data is inputted
    validate(need(data() != "", ""))

  x.all <- aldex.obj()

 if (input$effectplot_ab) {

        if (is.null(x.all)){
          return(NULL)
        } else {
          aldex.plot(x.all, type="MW", test="welch", all.cex = 1.5, rare.cex = 1.5, called.cex = 1.5, xlab = "Dispersion", ylab = "Difference")
          title(main = "Effect Plot")
        }
}
  })

output$effectMA <- renderPlot({

    # ensure data is inputted
    validate(need(data() != "", ""))

  x.all <- aldex.obj()
 if (input$effectplot_ab) {
    if (is.null(x.all)){
      return(NULL)
    } else {
      aldex.plot(x.all, type="MA", test="welch", all.cex = 1.5, rare.cex = 1.5, called.cex = 1.5, xlab = "CLR abundance", ylab = "Difference between")
      title(main = "Bland-Altman Plot")
    }
}
})

#effect plots for inputted aldex table
output$table_effect <- renderPlot({

    # ensure data is inputted
    validate(need(data() != "", ""))

  effect <- effect_input()
  point.colour <- input$point.colour

  if (input$effectplot_ab2) {
    if (is.null(effect)){
      return(NULL)
    } else {
      plot(effect$diff.win, effect$diff.btw, pch = 19, col=rgb(0,0,0,0.1), cex=0.4, xlab = "Difference within", ylab = "Difference between")
      title(main = "Effect Plot")
    }
    if (input$update_points == 0){
      return()
    } else {
      point.colour()
    }
  }
})

#effect plots for inputted aldex table
output$table_bland <- renderPlot({

    # ensure data is inputted
    validate(need(data() != "", ""))

  effect <- effect_input()
  point.colour <- input$point.colour

  if (input$effectplot_ab2) {
    if (is.null(effect)){
      return(NULL)
    } else {
      plot(effect$rab.all, effect$diff.btw, pch = 19, col=rgb(0,0,0,0.1), cex=0.4, xlab = "CLR abundance", ylab = "Difference within")
      title(main = "Bland-Altman Plot")
    }
    if (input$update_points == 0){
      return()
    } else {
      BA.point.colour()
    }
  }
})

#display information of hovered point
output$ma_hovertext <- renderUI({
  x.all <- aldex.obj()

  row <- nearPoints(x.all, input$ma_hover, xvar = "rab.all", yvar = "diff.btw")

  x <- paste("Feature Id: ", rownames(row))
  y <- paste("Median Log2 relative abundance: ", round(row$rab.all, digits =3))
  z <- paste("Between condition difference size: ", round(row$diff.btw, digits = 3))
  e <- paste("Effect size: ", round(row$effect, digits = 3))
  a <- ""
  HTML(paste(x, y, z, e, a, sep = "<br/>"))
})

#display information of hovered point
output$mw_hovertext <- renderUI({

    # ensure data is inputted
    validate(need(data() != "", ""))

  x.all <- aldex.obj()
  a <- d.clr()

  row <- nearPoints(x.all, input$mw_hover, xvar = "diff.win", yvar = "diff.btw", maxpoints = 1)

  if (is.null(input$mw_hover)) {
    HTML("Hover over MW plot with mouse")

  } else {

    p.v <- round(row$we.eBH, digits = 4)
    ps <- format(p.v, scientific=TRUE)
    x <- paste("Feature Id: ", rownames(row))
    y <- paste("Within condition difference size: ", round(row$diff.win, digits =3))
    z <- paste("Between condition difference size: ", round(row$diff.btw, digits = 3))
    e <- paste("Effect size: ", round(row$effect, digits = 3))
    p <- paste("Benjami Hochberg corrected p-value:", ps)
    a <- ""
    HTML(paste(x, y, z, e, p, a, sep = "<br/>"))
  }
})

#strip chart of expected CLR values for each sample per condition of hovered
# point
output$stripchart <- renderPlot({

    # ensure data is inputted
    validate(need(data() != "", "You need data. Input your data."))
    x.all <- aldex.obj()
    d.clr <- d.clr()

if (!is.null(input$mw_hover)){

feature <- nearPoints(x.all, input$mw_hover, xvar = "diff.win", yvar = "diff.btw", maxpoints = 1)

# only plot if one feature select
if (dim(feature)[1] == 1) {
aldex.plotFeature(d.clr, rownames(feature), pooledOnly=TRUE, densityOnly=TRUE)
}
} else {

    #HTML("Hover over MW plot with mouse")

    validate(need(input$mw_hover != "", "Hover over a point on the graph."))

}

})

#displays row name for plots when inputting ALDEx2 table
output$featurename<- renderUI({
  x.all <- effect_input()

  row <- nearPoints(x.all, input$mw_hover2, xvar = "diff.win", yvar = "diff.btw")

  ba.row <- nearPoints(x.all, input$ma_hover, xvar = "rab.all", yvar = "diff.btw")

  feature <- rownames(row)

  feature2 <- rownames(ba.row)

  number <- c(rownames(x.all[feature,]), rownames(x.all[feature2,]))
  return(number)
})

output$effectwarning <- renderText({
  "See the GitHub Wiki for more information on how to select conditions"
})
}