## Installation of Packages

my_packages <- c("shiny", "shinyalert", "data.table", "ggplot2", "shinyWidgets", "shinythemes", "dplyr", "wesanderson", "shinycssloaders", "forcats", "magick", "corrplot", "cowplot", "ggcorrplot", "shinyalert") 
## Find the non-installed ones 
not_installed <- my_packages[!(my_packages %in% installed.packages()[ , "Package"])]
## install them
if(length(not_installed)) install.packages(not_installed)                             

## Load Packages 

library(shiny)
library(shinyalert)
library(data.table)
library(ggplot2)
library(shinyWidgets)
library(shinythemes)
library(dplyr)
library(forcats)
library(wesanderson)
library(shinycssloaders)
library(magick)
library(corrplot)
library(ggcorrplot)
library(cowplot)


## Download Reference Distributions File 

epigen <- fread("./Training_Proteins.csv", check.names = F)
epigen = as.data.frame(epigen)
predicted_variables <- read.csv("./Predictors.csv")
tmp <- read.csv("./Shiny_Mean_Methylation_Dataset_Proteins.csv")
out <- fread("./Methylation_Scores_for_GSw3_Proteins.csv", check.names = F)
out = as.data.frame(out)

timeoutSeconds <- 180

inactivity <- sprintf("function idleTimer() {
                      var t = setTimeout(logout, %s);
                      window.onmousemove = resetTimer; // catches mouse movements
                      window.onmousedown = resetTimer; // catches mouse movements
                      window.onclick = resetTimer;     // catches mouse clicks
                      window.onscroll = resetTimer;    // catches scrolling
                      window.onkeypress = resetTimer;  //catches keyboard actions
                      
                      function logout() {
                      Shiny.setInputValue('timeOut', '%ss')
                      }
                      
                      function resetTimer() {
                      clearTimeout(t);
                      t = setTimeout(logout, %s);  // time is in milliseconds (1000 is 1 second)
                      }
                      }
                      idleTimer();", timeoutSeconds*1000, timeoutSeconds, timeoutSeconds*1000)

## User Interface Portion 
ui <- fluidPage(                tags$head(
  tags$style( ".title {margin-left: 550px; width: 500px}", 
              HTML(".shiny-notification {
                   position:fixed;
                   top: calc(50%);
                   left: calc(50%);
                   }
                   "
                   
              ),
              ".modal-body {padding: 20px}
              .modal-header {background-color: #8ae9ff}
              .modal { text-align: center}"
  )
              ),  tags$div(class="title", titlePanel("MethylDetectR")), 
  tags$script(inactivity), theme = shinytheme("united"), 
  sidebarLayout(
    sidebarPanel(
      useShinyalert(),
      fileInput("file1", "Choose CSV File",
                multiple = TRUE,
                accept = c("text/csv",
                           "text/comma-separated-values,text/plain",
                           ".csv")),
      
      fileInput("file2", "Upload Case/Control File",
                multiple = TRUE,
                accept = c("text/csv",
                           "text/comma-separated-values,text/plain",
                           ".csv")), 
    
      
      conditionalPanel(condition = "input.tabs == 1", 
    
      selectInput("var", 
                  label = "Choose an Epigenetic Predictor to Display",
                  choices = list("Biological Age" = c("Biological Age"), "Lifestyle" = c("Alcohol", 
                                                                   "Body Fat %",
                                                                   "Body Mass Index",
                                                                   "HDL Cholesterol",
                                                                   "Smoking",
                                                                   "Waist:Hip Ratio"), "Proteins" =c("CCL11",
                                                                                                     "CD6",
                                                                                                     "CRTAM",
                                                                                                     "CXCL10",
                                                                                                     "CXCL11",
                                                                                                     "CXCL9",
                                                                                                     "EN.RAGE",
                                                                                                     "FcRL2",
                                                                                                     "FGF.21",
                                                                                                     "G.CSF",
                                                                                                     "GDF.8",
                                                                                                     "GZMA",
                                                                                                     "HGF",
                                                                                                     "IL.18R1",
                                                                                                     "IL6",
                                                                                                     "MATN3",
                                                                                                     "MDGA1",
                                                                                                     "N.CDase",
                                                                                                     "NEP",
                                                                                                     "NMNAT1",
                                                                                                     "NTRK3",
                                                                                                     "OSM",
                                                                                                     "SIGLEC1",
                                                                                                     "SKR3",
                                                                                                     "SMPD1",
                                                                                                     "TGF.alpha",
                                                                                                     "TNFB",
                                                                                                     "TNFSF14",
                                                                                                     "VEGFA")),
                  multiple = F, selected = NULL),
      uiOutput("var_ui"),
      uiOutput("df_cases_choice")),
      
      
      conditionalPanel(condition = "input.tabs == 2", 
                       
                       selectizeInput("var1", 
                                   label = "Choose up to 10 Epigenetic Predictors to Display",
                                   choices = list("Lifestyle" = c("Alcohol", 
                                                                                    "Body Fat %",
                                                                                    "Body Mass Index",
                                                                                    "HDL Cholesterol",
                                                                                    "Smoking",
                                                                                    "Waist:Hip Ratio"), "Proteins" =c("CCL11",
                                                                                                                      "CD6",
                                                                                                                      "CRTAM",
                                                                                                                      "CXCL10",
                                                                                                                      "CXCL11",
                                                                                                                      "CXCL9",
                                                                                                                      "EN.RAGE",
                                                                                                                      "FcRL2",
                                                                                                                      "FGF.21",
                                                                                                                      "G.CSF",
                                                                                                                      "GDF.8",
                                                                                                                      "GZMA",
                                                                                                                      "HGF",
                                                                                                                      "IL.18R1",
                                                                                                                      "IL6",
                                                                                                                      "MATN3",
                                                                                                                      "MDGA1",
                                                                                                                      "N.CDase",
                                                                                                                      "NEP",
                                                                                                                      "NMNAT1",
                                                                                                                      "NTRK3",
                                                                                                                      "OSM",
                                                                                                                      "SIGLEC1",
                                                                                                                      "SKR3",
                                                                                                                      "SMPD1",
                                                                                                                      "TGF.alpha",
                                                                                                                      "TNFB",
                                                                                                                      "TNFSF14",
                                                                                                                      "VEGFA")),
                                   multiple = T,options = list(maxItems = 10), selected = NULL),
                       
                       
                       radioButtons("choice_percentile", "Choose to Display ID-level Percentiles or Cases vs. Controls", choices = c("Individual Level", "Cases vs. Controls")),
                       uiOutput("var_ui2"),
                       uiOutput("df_cases_choice2")),
      
      
      conditionalPanel(condition = "input.tabs == 3", 
                       
                       selectizeInput("var2", 
                                      label = "Choose up to 10 Epigenetic Predictors to Display",
                                      choices = list("Lifestyle" = c("Alcohol", 
                                                                     "Body Fat %",
                                                                     "Body Mass Index",
                                                                     "HDL Cholesterol",
                                                                     "Smoking",
                                                                     "Waist:Hip Ratio"), "Proteins" =c("CCL11",
                                                                                                       "CD6",
                                                                                                       "CRTAM",
                                                                                                       "CXCL10",
                                                                                                       "CXCL11",
                                                                                                       "CXCL9",
                                                                                                       "EN.RAGE",
                                                                                                       "FcRL2",
                                                                                                       "FGF.21",
                                                                                                       "G.CSF",
                                                                                                       "GDF.8",
                                                                                                       "GZMA",
                                                                                                       "HGF",
                                                                                                       "IL.18R1",
                                                                                                       "IL6",
                                                                                                       "MATN3",
                                                                                                       "MDGA1",
                                                                                                       "N.CDase",
                                                                                                       "NEP",
                                                                                                       "NMNAT1",
                                                                                                       "NTRK3",
                                                                                                       "OSM",
                                                                                                       "SIGLEC1",
                                                                                                       "SKR3",
                                                                                                       "SMPD1",
                                                                                                       "TGF.alpha",
                                                                                                       "TNFB",
                                                                                                       "TNFSF14",
                                                                                                       "VEGFA")),
                                      multiple = T, options = list(maxItems = 10), selected = NULL),
                       
                uiOutput("radioButtons_case_control"),
                uiOutput("df_cases_choice3")),
      
      conditionalPanel(condition = "input.tabs == 4", 
                       NULL
      ), 
      
      
    selectInput("sexInput", "Sex", choices = c("NULL", "Male", "Female"), selected = NULL),
    uiOutput("ageslider")),
    
    mainPanel(tabsetPanel(id="tabs", tabPanel("Density Plot", value = 1, withSpinner(plotOutput("plot"), type = 8), br(), uiOutput("downloadinput"), br(), br(), br(), uiOutput("message1")),
                          tabPanel("Percentile Plot", value = 2, withSpinner(plotOutput("results"),type=8), textOutput("text"),br(), br(), br(), br(), uiOutput("message2")), 
                          tabPanel("Correlation between Predictors", value = 3, withSpinner(plotOutput("corrplot"),type = 8), br(), br(), br(), uiOutput("message3")), 
                          tabPanel("Correlation between Predictors and Traits", value = 4, withSpinner(plotOutput("corrplot2"), type = 8),  br(), br(), br(), uiOutput("message4"))))
  ))
## Server Portion 

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

  
    predictors <- predicted_variables$Predictors
  

  df1 <- 
    reactive({ 
      req(input$file1)
      df1 <- fread(input$file1$datapath, check.names = F)
      df1 = as.data.frame(df1)
      return(df1)
    }) 
  
 

  
  
  observeEvent(input$file1, { 
    means1 = sapply(df1()[,4:36], is.numeric)
    means1 = as.data.frame(means1)
    
    if(nrow(df1()) == 0){ 
      shinyalert("Sorry, An Error Has Occurred","There Are No Individuals In The Dataset!",type = "warning",  showCancelButton = TRUE, inputId = "foo1")
    } else if(length(which(means1[,1] %in% "FALSE")) > 0) 
    {
      # show pop-up ...
      shinyalert("Sorry, An Error Has Occurred","Non-Numeric Values Are Detected For the Predictors!",type = "warning",  showCancelButton = TRUE, inputId = "foo1")
      observeEvent(input$foo1, {
        session$reload()
      })
    } else if(length(which(!names(df1()) %in% names(epigen))) > 0) { 
      # show pop-up ...
      shinyalert("Sorry, An Error Has Occurred", paste("The Following Column Names Are Not Present in the Correct Format:", paste(names(epigen)[-which(names(epigen) %in% names(df1()))], collapse = ", "), sep = " "),type = "warning",  showCancelButton = TRUE, inputId = "foo")
      observeEvent(input$foo, {
        session$reload()
      })
    } else { NULL } } )
  
  
  df2 <- 
    reactive({ req(input$file1)
      df2 <- read.csv(input$file1$datapath, check.names=F)
      return(df2)
    }) 
  
  df_cases <- 
    reactive({ req(input$file2)
      df_cases1 <- read.csv(input$file2$datapath, check.names=F)
      return(df_cases1)
    }) 
  
  
  
  output$message1 <- renderUI(HTML(paste(
    "Please remember to cite McCartney & Hillary", em("et al."), "(2018),", "Zhang", em("et al."), "(2019)", "and Gadd",  em("et al."), "(2020)"
  )))
  
  output$message2 <- renderUI(HTML(paste(
    "Please remember to cite McCartney & Hillary", em("et al."), "(2018),", "Zhang", em("et al."), "(2019)", "and Gadd",  em("et al."), "(2020)"
  )))
  
  output$message3 <- renderUI(HTML(paste(
    "Please remember to cite McCartney & Hillary", em("et al."), "(2018),", "Zhang", em("et al."), "(2019)", "and Gadd",  em("et al."), "(2020)"
  )))
  
  output$message4 <- renderUI(HTML(paste(
    "Please remember to cite McCartney & Hillary", em("et al."), "(2018),", "Zhang", em("et al."), "(2019)", "and Gadd",  em("et al."), "(2020)"
  )))
  
  
  
  output$df_cases_choice <- renderUI({
    selectInput("Case_Input",
                "Choose a Case/Control Variable", 
                choices = c("NULL", names(df_cases()[2:ncol(df_cases())])),
                multiple = F, selected = NULL) })
  
  
  output$df_cases_choice2 <- renderUI({
    selectInput("Case_Input2",
                "Choose a Case/Control Variable", 
                choices = c("NULL", names(df_cases()[2:ncol(df_cases())])),
                multiple = F, selected = NULL) })
  
  output$df_cases_choice3 <- renderUI({
    selectInput("Case_Input3",
                "Choose a Case/Control Variable", 
                choices = c("NULL", names(df_cases()[2:ncol(df_cases())])),
                multiple = F, selected = NULL) })
  
  
  output$radioButtons_case_control <- renderUI({ 
    if(is.null(input$file2)) { 
      NULL
    } else if(is.null(input$Case_Input3) || input$Case_Input3 %in% "NULL") {
      NULL
    } else { 
      radioButtons("case_control", "Subset by Case/Controls or Plot Case/Controls side by side", choices = c("Subset by Cases", "Subset by Controls", "Show Cases vs. Controls"))
    }
  })
  
  
  output$ageslider <- renderUI({
    
    min.epigen <- round(min(epigen$`True Age`, na.rm = T),1)
    max.epigen <- round(max(epigen$`True Age`, na.rm = T),1)
    
    min.uploaded = round(min(c(df1()$'True Age', df1()$'Biological Age'), na.rm = T),1)
    max.uploaded = round(max(c(df1()$'True Age', df1()$'Biological Age'), na.rm = T),1)
    
    min.value <- ifelse(min.epigen < min.uploaded, min.epigen, min.uploaded)
    max.value <- ifelse(max.epigen > max.uploaded, max.epigen, max.uploaded)
    
    
    sliderInput("ageInput", "Age", min = min.value, max = max.value, value = c(min.value,max.value))
  })
  
  
  output$var_ui <- renderUI({
    selectizeInput("IDInput",
                "ID", 
                choices = df1()$ID, multiple = F, options = list(maxOptions = length(df1()$ID))) }) 
  
  
  output$var_ui2 <- renderUI({
    selectizeInput("IDInput2",
                   "ID", 
                   choices = df1()$ID, multiple = F, options = list(maxOptions = length(df1()$ID))) }) 
  
  
  dat <- reactive({
    if(all(is.na(df1()[,"True Age"]))){
      
    if(input$sexInput %in% "NULL") { 
      test1 <- df1()[which(df1()$'Biological Age' >= input$ageInput[1] & df1()$'Biological Age' <= input$ageInput[2]),]  
      return(test1) }  
      
      else  {
        test2 <- df1()[which(df1()$'Biological Age' >= input$ageInput[1] & df1()$'Biological Age' <= input$ageInput[2] & df1()$Sex %in% input$sexInput),]
       return(test2) } 
    }
      else { 
         if(input$sexInput %in% "NULL") { 
           has_true_age <- df1()[which(!is.na(df1()$'True Age')), ]
           trueage <- has_true_age[which(has_true_age$'True Age' >= input$ageInput[1] & has_true_age$'True Age' <= input$ageInput[2]),]  
           
          has_biological_age <- df1()[which(is.na(df1()$'True Age')), ]
          biologicalage <- has_biological_age[which(has_biological_age$'Biological Age' >= input$ageInput[1] & has_biological_age$'Biological Age' <= input$ageInput[2]),]  
          
          combinedage <- rbind(trueage, biologicalage)
         
           return(combinedage) }  
        
        else  { 
             has_true_age <- df1()[which(!is.na(df1()$'True Age')), ]
           trueage <- has_true_age[which(has_true_age$'True Age' >= input$ageInput[1] & has_true_age$'True Age' <= input$ageInput[2] & has_true_age$Sex %in% input$sexInput),]  
           
           has_biological_age <- df1()[which(is.na(df1()$'True Age')), ]
           biologicalage <- has_biological_age[which(has_biological_age$'Biological Age' >= input$ageInput[1] & has_biological_age$'Biological Age' <= input$ageInput[2] & has_biological_age$Sex %in% input$sexInput),]  
           
           combinedage_sex <- rbind(trueage, biologicalage)
           return(combinedage_sex)
         
         }
       }
    
  }) 
  

  cases_split <- reactive({ 
    cases.split <- merge(dat(), df_cases(), by = "ID")
    return(cases.split)
    })
  
  cases_split_cases <- reactive({ 
    cases.split.cases <- merge(dat(), df_cases()[which(df_cases()[,as.character(input$Case_Input3)] == 1),], by = "ID")
    return(cases.split.cases)
  })
  
  cases_split_controls <- reactive({ 
    cases.split.controls <- merge(dat(), df_cases()[which(df_cases()[,as.character(input$Case_Input3)] == 0),], by = "ID")
    return(cases.split.controls)
  })

  
  
  dat_epigen <- reactive({
    if (input$sexInput %in% "NULL") { 
      epigen1 <- epigen[which(epigen$'True Age' >= input$ageInput[1] & epigen$'True Age' <= input$ageInput[2]),]  
      return(epigen1)}  else  
        epigen2 <- epigen[which(epigen$'True Age' >= input$ageInput[1] & epigen$'True Age' <= input$ageInput[2] & epigen$Sex %in% input$sexInput),]
      return(epigen2)
  })
  
  
  
  
  dat2 <- reactive({
    dat_input <- dat()[which(dat()$ID %in% input$IDInput),]
    full_set <- rbind(dat_input,epigen)
    return(full_set)
  })
  
  
  
  
  id <- reactive({ individual <- which(dat()$ID %in% input$IDInput2)
  return(individual)
  })
  
  id1 <- reactive({ individual1 <- which(df1()$ID %in% input$IDInput)
  return(individual1)
  })
  
  
  id2 <- reactive({ individual2 <- dat()[which(dat()$ID %in% input$IDInput),"ID"] 
  return(individual2)
  })
  
  
  dat_for_plot <- reactive({ 
    dat_plot <- rbind(df1()[which(df1()$ID %in% input$IDInput),], dat())
    return(dat_plot)
  })
  
  
  dat_for_percentile <- reactive({ 
    dat_percentile <- rbind(df1()[which(df1()$ID %in% input$IDInput2),], dat())
    return(dat_percentile)
  })
  
  cases_split_plot <- reactive({ 
cases <- merge(dat_for_plot(), df_cases(), by = "ID")    
return(cases)
  })
  

  
  table_input_rbind <- reactive({ 
    reference = cbind(dat_for_plot()$ID, as.data.frame(apply(dat_for_plot()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))) 
    row.names(reference) <- reference[,1]
    reference[,1] <- NULL 
    reference <- as.data.frame(t(reference))
    return(reference)
  })
  
  id_plot <- reactive({ 
    id_ggplot <- dat_for_plot()[which(dat_for_plot()$ID %in% input$IDInput), "ID"] 
    return(id_ggplot) 
  })
  
  id_plot_panel2 <- reactive({ 
    id_ggplot <- dat_for_percentile()[which(dat_for_percentile()$ID %in% input$IDInput2), "ID"] 
    return(id_ggplot) 
  })
  
  table_input <- reactive({ if(nrow(dat()) > 0) { 
    reference = cbind(dat()$ID, as.data.frame(apply(dat()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))) 
    row.names(reference) <- reference[,1]
    reference[,1] <- NULL 
    reference <- as.data.frame(t(reference))
    return(reference)
  } else print("0")
  })
  
  dat_for_table <- reactive({ 
    dat_table <- rbind(df1()[which(df1()$ID %in% input$IDInput2),], dat())
    return(dat_table)
  })
  
  
  id_table <- reactive({ 
    id_table1 <- which(dat_for_table()$ID %in% input$IDInput2)
    return(id_table1) 
  })
  
  
  id_table2 <- reactive({ 
    id_table2 <- which(dat_for_plot()$ID %in% input$IDInput2)
    return(id_table2) 
  })
  
  
  table_rbind <- reactive({
    reference = cbind(dat_for_table()$ID, as.data.frame(apply(dat_for_table()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))) 
    row.names(reference) <- reference[,1]
    reference[,1] <- NULL 
    reference <- as.data.frame(t(reference))
    return(reference) }) 
  
  table <- reactive({ 
    reference = cbind(dat()$ID, as.data.frame(apply(dat()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))) 
    row.names(reference) <- reference[,1]
    reference[,1] <- NULL 
    reference <- as.data.frame(t(reference))
    return(reference)
  })
  
  
  case_control_percentile <- reactive({ 
    
    case_control <- table()[,which(colnames(table()) %in% df_cases()[which(df_cases()[,as.character(input$Case_Input2)] == 1), "ID"])]
    return(case_control)
    })
  
  median_iqr_case_control <- reactive({ 
    
    median.cases = as.data.frame(apply(case_control_percentile(),1,function(x) median(x, na.rm = T)))
    q25.cases = as.data.frame(apply(case_control_percentile(),1,function(x) quantile(x, prob = 0.25, na.rm = T)))
    q75.cases = as.data.frame(apply(case_control_percentile(),1,function(x) quantile(x, prob = 0.75, na.rm = T)))
    total = cbind(median.cases, cbind(q25.cases, q75.cases))
    names(total) <- c("Median", "Q25", "Q75")
    total1 = total[which(row.names(total) %in% as.character(input$var1)),]
    return(total1)
    })
  
  
  for_table <- reactive({ 
    tmp <- cases_split()[-which(cases_split()$ID %in% input$IDInput),]
        tmp <- tmp[which(!is.na(tmp[,as.character(input$Case_Input)])),]
    return(tmp) 
    })
  
  for_table_plot <- reactive({ 
    tmp <- cases_split_plot()[-which(cases_split_plot()$ID %in% input$IDInput),]
    tmp <- tmp[which(!is.na(tmp[,as.character(input$Case_Input)])),]
  })
  
  for_table_cases_rbind <- reactive({ 
    dat_plot_cases <- rbind(cases_split()[which(cases_split()$ID %in% input$IDInput),], for_table()[which(for_table()[,as.character(input$Case_Input)] == 1),])
    return(dat_plot_cases)
  })
  
  for_table_cases_rbind_plot <- reactive({ 
    dat_plot_cases.2 <- rbind(cases_split_plot()[which(cases_split_plot()$ID %in% input$IDInput),], for_table_plot()[which(for_table_plot()[,as.character(input$Case_Input)] == 1),])
    return(dat_plot_cases.2)
  })
  
  
  for_table_controls_rbind <- reactive({ 
    dat_plot_controls <- rbind(cases_split()[which(cases_split()$ID %in% input$IDInput),], for_table()[which(for_table()[,as.character(input$Case_Input)] == 0),])
    return(dat_plot_controls)
  })
  
  for_table_controls_rbind_plot <- reactive({ 
    dat_plot_controls.2 <- rbind(cases_split_plot()[which(cases_split_plot()$ID %in% input$IDInput),], for_table_plot()[which(for_table_plot()[,as.character(input$Case_Input)] == 0),])
    return(dat_plot_controls.2)
  })

  
  
  cases_table <- reactive({ 
    if(nrow(for_table_cases_rbind()) > 1) { 
    reference = cbind(for_table_cases_rbind()$ID, as.data.frame(apply(for_table_cases_rbind()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))) 
    row.names(reference) <- reference[,1]
    reference[,1] <- NULL 
    reference <- as.data.frame(t(reference))
    return(reference)
    } else { 
    reference = cbind(for_table_cases_rbind()$ID, apply(for_table_cases_rbind()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))
    reference = as.data.frame(reference)
    reference[,1] <- NULL
    names(reference) <- for_table_cases_rbind()$ID
    reference[,1] <- 0
    return(reference)
      }
  })
  
  controls_table <- reactive({ 
    if(nrow(for_table_controls_rbind()) > 1) { 
      reference = cbind(for_table_controls_rbind()$ID, as.data.frame(apply(for_table_controls_rbind()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))) 
      row.names(reference) <- reference[,1]
      reference[,1] <- NULL 
      reference <- as.data.frame(t(reference))
      return(reference)
    } else { 
      reference = cbind(for_table_controls_rbind()$ID, apply(for_table_controls_rbind()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))
      reference = as.data.frame(reference)
      reference[,1] <- NULL
      names(reference) <- for_table_controls_rbind()$ID
      reference[,1] <- 0
      return(reference)
    }
  })
  
  cases_table_plot <- reactive({ 
    if(nrow(for_table_cases_rbind_plot()) > 1) { 
      reference = cbind(for_table_cases_rbind_plot()$ID, as.data.frame(apply(for_table_cases_rbind_plot()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))) 
      row.names(reference) <- reference[,1]
      reference[,1] <- NULL 
      reference <- as.data.frame(t(reference))
      return(reference)
    } else { 
      reference = cbind(for_table_cases_rbind_plot()$ID, apply(for_table_cases_rbind_plot()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))
      reference = as.data.frame(reference)
      reference[,1] <- NULL
      names(reference) <- for_table_cases_rbind_plot()$ID
      reference[,1] <- 0
      return(reference)
    }
  })
  
  
  controls_table_plot <- reactive({ 
    if(nrow(for_table_controls_rbind_plot()) > 1) { 
      reference = cbind(for_table_controls_rbind_plot()$ID, as.data.frame(apply(for_table_controls_rbind_plot()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))) 
      row.names(reference) <- reference[,1]
      reference[,1] <- NULL 
      reference <- as.data.frame(t(reference))
      return(reference)
    } else { 
      reference = cbind(for_table_controls_rbind_plot()$ID, apply(for_table_controls_rbind_plot()[,as.character(predictors)], 2, function(x) percent_rank(x)*100))
      reference = as.data.frame(reference)
      reference[,1] <- NULL
      names(reference) <- for_table_controls_rbind_plot()$ID
      reference[,1] <- 0
      return(reference)
    }
  })
  
  
  
  id_cases <- reactive({ individual <- for_table_cases_rbind()[which(for_table_cases_rbind()$ID %in% input$IDInput), "ID"]
  return(individual)
  })
  
  id_controls <- reactive({ individual.2 <- for_table_controls_rbind()[which(for_table_controls_rbind()$ID %in% input$IDInput), "ID"]
  return(individual.2)
  })
  
  rounded_rbind <- reactive({ rounded_value_rbind <- round(table_rbind()[,id_table()]) 
  return(rounded_value_rbind)
  })
  
  rounded <- reactive({ rounded_value <- round(table()[,id()]) 
  return(rounded_value)
  })
  
  
  
  output$plot <- renderPlot({
    
    
    pick <- switch(input$var, 
                   "Alcohol" = "Alcohol",
                   "Smoking" = "Smoking", 
                   "Body Mass Index" = "Body Mass Index",
                   "Body Fat %" = "Body Fat %",
                   "HDL Cholesterol" = "HDL Cholesterol",
                   "Waist:Hip Ratio" = "Waist:Hip Ratio",
                   "CCL11" = "CCL11",
                   "CD6" = "CD6",
                   "CRTAM" = "CRTAM",
                   "CXCL10" = "CXCL10",
                   "CXCL11" = "CXCL11",
                   "CXCL9"="CXCL9",
                   "EN.RAGE"="EN.RAGE",
                   "FcRL2"="FcRL2",
                   "FGF.21"="FGF.21",
                   "G.CSF"="G.CSF",
                   "GDF.8"="GDF.8",
                   "GZMA"="GZMA",
                   "HGF"="HGF",
                   "IL.18R1"="IL.18R1",
                   "IL6"="IL6",
                   "MATN3"="MATN3",
                   "MDGA1"="MDGA1",
                   "N.CDase"="N.CDase",
                   "NEP"="NEP",
                   "NMNAT1"="NMNAT1",
                   "NTRK3"="NTRK3",
                   "OSM"="OSM",
                   "SIGLEC1"="SIGLEC1",
                   "SKR3"="SKR3",
                   "SMPD1"="SMPD1",
                   "TGF.alpha"="TGF.alpha",
                   "TNFB"="TNFB",
                   "TNFSF14"="TNFSF14",
                   "VEGFA"="VEGFA")
    

    if(input$var %in% "Biological Age"){ 
      
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0,paste("We predict that you are","\n", round(df1()$'Biological Age'[id1()],1), "years old!", sep = " "), cex = 1.7)
    } 
    
    else if(nrow(dat()) <= 1){ 
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0,"There are too few individuals selected in dataset for density plot!", cex = 1.7)
    } 
    else if(is.null(input$Case_Input) || input$Case_Input %in% "NULL") { 
    
      
      if(length(which(dat()$ID %in% input$IDInput))==0) { 
        
        
        p = ggplot(dat_for_plot(), aes(dat_for_plot()[,as.character(pick)])) + 
          geom_density(aes(fill = "blue"),alpha=0.2) + ylab("Density") 
        
        p= p + theme(legend.title = element_text(size = 13), legend.text = element_text(size=12), axis.text = element_text(size = 14), axis.title.y = element_text(size = 13), axis.title.x=element_text(size = 14, margin=margin(t=10,r=0,b=10,l=0)))
        to_fill <- tibble(
          x = ggplot_build(p)$data[[1]]$x,
          y = ggplot_build(p)$data[[1]]$y)
        
        minimum <- min(dat_for_plot()[,as.character(pick)], na.rm = T)
        maximum <- max(dat_for_plot()[,as.character(pick)], na.rm = T)
        
        
        p + geom_vline(xintercept = dat_for_plot()[,as.character(pick)][which(dat_for_plot()[,"ID"]== id_plot())], linetype = "dashed") + theme(legend.title = element_text(hjust = 0.5)) + geom_segment(aes(x=minimum, xend = maximum, y=0, yend = 0), size=1, arrow = arrow(length = unit(0.6,"cm"))) +
          xlab(paste("Methylation-based Scores for", as.character(pick))) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.line = element_line("black")) + scale_fill_discrete(name = "Dataset",labels = c(paste(paste(paste("Input", "(n = ", sep = " "), nrow(dat()), ")", sep = ""),"\n Percentile vs Input: ", ifelse(nrow(dat()) > 0, as.character(round(table_input_rbind()[as.character(pick),which(colnames(table_input_rbind()) %in% id_plot())]),2), "NA")))) } 
      
      
      else {
        p = ggplot(dat(), aes(dat()[,as.character(pick)])) + 
          geom_density(aes(fill = "blue"),alpha=0.2) + ylab("Density") 
        
        p= p + theme(legend.title = element_text(size = 13), legend.text = element_text(size=12), axis.text = element_text(size = 14), axis.title.y = element_text(size = 13), axis.title.x=element_text(size = 14, margin=margin(t=10,r=0,b=10,l=0)))
        to_fill <- tibble(
          x = ggplot_build(p)$data[[1]]$x,
          y = ggplot_build(p)$data[[1]]$y)
        
        minimum <- min(dat()[,as.character(pick)], na.rm = T)
        maximum <- max(dat()[,as.character(pick)], na.rm = T)
        
        
        p  +  geom_vline(xintercept = dat()[,as.character(pick)][which(dat()[,"ID"]== id2())], linetype = "dashed") + theme(legend.title = element_text(hjust = 0.5)) + geom_segment(aes(x=minimum, xend = maximum, y=0, yend = 0), size=1, arrow = arrow(length = unit(0.6,"cm"))) +
          xlab(paste("Methylation-based Scores for", as.character(pick))) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.line = element_line("black")) + scale_x_continuous(limits = c(minimum,maximum)) + scale_fill_discrete(name = "Dataset",labels = c(paste(paste(paste("Input", "(n = ", sep = " "), nrow(dat()), ")", sep = ""),"\n Percentile vs Input: ", ifelse(nrow(dat()) > 0, as.character(round(table_input()[as.character(pick),which(colnames(table_input()) %in% id2())])), "0")))) 
        
      } 
      
    }
    
    
    else {
    
      if(length(which(!is.na(df_cases()[,as.character(input$Case_Input)]))) < 2 || length(which(df_cases()[,as.character(input$Case_Input)] %in% "NA")) > (nrow(df_cases())-2) ) 
      { 
        plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
        text(0,0,"There are too few individuals with case/control information for density plot!", cex = 1.7)
      } else { 
    if(length(which(cases_split()$ID %in% input$IDInput))==0) { 
      
      
      p = ggplot(cases_split_plot()[which(!is.na(cases_split_plot()[,as.character(input$Case_Input)])),], aes(cases_split_plot()[which(!is.na(cases_split_plot()[,as.character(input$Case_Input)])),as.character(pick)], fill = factor(cases_split_plot()[which(!is.na(cases_split_plot()[,as.character(input$Case_Input)])),as.character(input$Case_Input)]))) + 
        geom_density(alpha=0.2) + ylab("Density") 
      
      p= p + theme(legend.title = element_text(size = 13), legend.text = element_text(size=12), axis.text = element_text(size = 14), axis.title.y = element_text(size = 13), axis.title.x=element_text(size = 14, margin=margin(t=10,r=0,b=10,l=0)))
      to_fill <- tibble(
        x = ggplot_build(p)$data[[1]]$x,
        y = ggplot_build(p)$data[[1]]$y)
      
      minimum <- min(cases_split_plot()[which(!is.na(cases_split_plot()[,as.character(input$Case_Input)])),as.character(pick)], na.rm = T)
      maximum <- max(cases_split_plot()[which(!is.na(cases_split_plot()[,as.character(input$Case_Input)])),as.character(pick)], na.rm = T)
      
      
      p + geom_vline(xintercept = dat_for_plot()[,as.character(pick)][which(dat_for_plot()[,"ID"]== id_plot())], linetype = "dashed") + theme(legend.title = element_text(hjust = 0.5)) + geom_segment(aes(x=minimum, xend = maximum, y=0, yend = 0), size=1, arrow = arrow(length = unit(0.6,"cm"))) +
        xlab(paste("Methylation-based Scores for", as.character(pick))) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.line = element_line("black")) + scale_x_continuous(limits = c(minimum-0.001,maximum+0.001)) +
        scale_fill_discrete(name = "Group",labels = rev(c(paste(paste(paste("Cases", "(n = ", sep = " "), length(which(for_table_plot()[,as.character(input$Case_Input)] == 1)), ")", sep = ""),"\n Percentile vs Cases: ", as.character(round(cases_table_plot()[as.character(pick),which(colnames(cases_table_plot()) %in% id_plot())]))), paste(paste(paste("Controls", "(n = ", sep = " "),length(which(for_table_plot()[,as.character(input$Case_Input)] == 0)), ")", sep = ""), "\n Percentile vs Controls: ", as.character(round(controls_table_plot()[as.character(pick),which(colnames(controls_table_plot()) %in% id_plot())])), sep= ""))))



    } 
    
    else {
      p = ggplot(cases_split()[which(!is.na(cases_split()[,as.character(input$Case_Input)])),], aes(cases_split()[which(!is.na(cases_split()[,as.character(input$Case_Input)])),as.character(pick)], fill = factor(cases_split()[which(!is.na(cases_split()[,as.character(input$Case_Input)])),as.character(input$Case_Input)]))) + 
        geom_density(alpha = 0.2) + ylab("Density") 
      
      p= p + theme(legend.title = element_text(size = 13), legend.text = element_text(size=12), axis.text = element_text(size = 14), axis.title.y = element_text(size = 13), axis.title.x=element_text(size = 14, margin=margin(t=10,r=0,b=10,l=0)))
      to_fill <- tibble(
        x = ggplot_build(p)$data[[1]]$x,
        y = ggplot_build(p)$data[[1]]$y)
      
      minimum <- min(cases_split()[which(!is.na(cases_split()[,as.character(input$Case_Input)])),as.character(pick)], na.rm = T)
      maximum <- max(cases_split()[which(!is.na(cases_split()[,as.character(input$Case_Input)])),as.character(pick)], na.rm = T)
      
      
      p +  geom_vline(xintercept = dat()[,as.character(pick)][which(dat()[,"ID"]== id2())], linetype = "dashed") + theme(legend.title = element_text(hjust = 0.5)) + geom_segment(aes(x=minimum, xend = maximum, y=0, yend = 0), size=1, arrow = arrow(length = unit(0.6,"cm"))) +
        xlab(paste("Methylation-based Scores for", as.character(pick))) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.line = element_line("black")) + scale_x_continuous(limits = c(minimum,maximum)) +
        scale_fill_discrete(name = "Group",labels = rev(c(paste(paste(paste("Cases", "(n = ", sep = " "), length(which(for_table_cases_rbind()[,as.character(input$Case_Input)] == 1)), ")", sep = ""),"\n Percentile vs Cases: ", as.character(round(cases_table()[as.character(pick),which(colnames(cases_table()) %in% id_plot())]))), paste(paste(paste("Controls", "(n = ", sep = " "),length(which(for_table_controls_rbind()[,as.character(input$Case_Input)] == 0)), ")", sep = ""), "\n Percentile vs Controls: ", as.character(round(controls_table()[as.character(pick),which(colnames(controls_table()) %in% id_plot())])), sep= ""))))
      
    } 
      
      
      }
    }
  })
  
  table_rbind2 <- reactive({ t2 <- table_rbind()[which(row.names(table_rbind()) %in% input$var1),] 
  return(t2)
  }) 
  
  table2 <- reactive({ tab2 <- table()[which(row.names(table()) %in% input$var1),] 
  return(tab2)
  }) 
  
  rounded_rbind2 <- reactive({ rounded_value_rbind <- round(table_rbind2()[,id_table()]) 
  return(rounded_value_rbind)
  })
  

  
  rounded2 <- reactive({ rounded_value <- round(table2()[,id()]) 
  return(rounded_value)
  })
  
  
  predictions <- reactive({ p <- as.character(predictors[which(predictors %in% input$var1)])
  return(p)
  }) 
  predictions2 <- reactive({ p <- as.character(predictors[which(predictors %in% input$var2)])
  return(p)
  }) 
  
  output$results <- renderPlot({ 
    
    if(input$choice_percentile %in% "Individual Level") { 
    
    if(length(predictions()) < 1) { 
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0,"Please select up to 10 predictors to view percentile ranks!", cex = 1.7)
    }
    else if(nrow(dat()) <= 2) { 
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0, paste("Number of individuals selected of", nrow(dat()), "is too few for analysis!"), cex = 1.7)
      }
    else if(length(which(dat()$ID %in% input$IDInput2))==0)  { 
      ggplot(table_rbind2(), aes(x = table_rbind2()[,id_table()], y = fct_rev(row.names(table_rbind2())), colour = row.names(table_rbind2()))) + geom_point(size = 14) + 
        geom_text(label = rounded_rbind2(), colour = "white", size = 6.3) + theme(axis.text = element_text(size = 14), axis.title = element_text(size = 14), axis.title.y=element_text(margin = margin(t=0,r=15,b=0,l=0))) +
        xlab(paste("Percentile for",id_plot())) + ylab("Epigenetic Predictor") + 
        theme(legend.position = "none") + xlim(c(0,100)) + scale_y_discrete(limits = levels(row.names(table_rbind2())))} 
    
    
    else { 
          ggplot(table2(), aes(x = table2()[,id()], y = fct_rev(row.names(table2())), colour = row.names(table2()))) + geom_point(size = 14) + 
      geom_text(label = rounded2(), colour = "white", size = 6.3) + theme(axis.text = element_text(size = 14), axis.title = element_text(size = 14), axis.title.y=element_text(margin = margin(t=0,r=15,b=0,l=0))) +
      xlab(paste("Percentile for",names(table2())[id()])) + ylab("Epigenetic Predictor") + 
      theme(legend.position = "none") + xlim(c(0,100)) + scale_y_discrete(limits = levels(row.names(table2())))
    } 
      
    }  else { 
      
      if(length(predictions()) < 1) { 
        plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
        text(0,0,"Please select up to 10 predictors to view percentile ranks!", cex = 1.7)
      }
  
      else if(is.null(input$Case_Input2) || input$Case_Input2 %in% "NULL") { 
        plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
        text(0,0,"Please upload and select a Case/Control Variable to subset data!", cex = 1.7)
        
        }
      else if(length(which(cases_split()[,as.character(input$Case_Input2)] == 1)) <= 1){ 
        plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
        text(0,0, "Number of individuals selected is too few for analysis!", cex = 1.7)
      }
      else { 
      ggplot(median_iqr_case_control(), aes(x = median_iqr_case_control()$Median, y = fct_rev(row.names(median_iqr_case_control())), colour = row.names(median_iqr_case_control()))) + geom_point(size = 14) + 
        geom_errorbarh(aes(xmax = median_iqr_case_control()$Q75, xmin = median_iqr_case_control()$Q25, height = 0.2)) + geom_text(label = round(median_iqr_case_control()$Median), colour = "white", size = 6.3) + 
        theme(axis.text = element_text(size = 14), axis.title = element_text(size = 14), axis.title.y=element_text(margin = margin(t=0,r=15,b=0,l=0))) +
        xlab("Median Percentile for Cases [Interquartile Range]") + ylab("Epigenetic Predictor") + 
        theme(legend.position = "none") + xlim(c(0,100)) + scale_y_discrete(limits = levels(row.names(median_iqr_case_control()))) 
    
        } 
      }
      })

  
  output$corrplot <- renderPlot({ 
    if(length(predictions2()) < 2) { 
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0,"Please select between 2 and 10 predictors to view their relationships!", cex = 1.7)
      }
    else if (nrow(dat_epigen()) <= 10 ) {
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0,"Sorry, the number of reference individuals must be above 10!", cex = 1.7) 
    } else {
      
      
    if(is.null(input$Case_Input3) || input$Case_Input3 %in% "NULL")  {
    
    if(nrow(dat_epigen()) > 10 & nrow(dat()) >= 3) { 
      cormatrix_gs <- cor(dat_epigen()[,predictions2()])
      cormatrix_input <- cor(dat()[,predictions2()])
      corrplot_gs <- ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T)  + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
      corrplot_input <- ggcorrplot(cormatrix_input,type="lower", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Input Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
      legend_corr <- (get_legend(ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
      corr_together <- plot_grid(plot_grid(corrplot_input, corrplot_gs, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
      print(corr_together)} else { 
        cormatrix_gs <- cor(dat_epigen()[,predictions2()])
        corrplot_gs <- ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation",colors = rev(c("#6b81ff", "white", "#E46726")), lab = T)  + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
        legend_corr <- (get_legend(ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation",colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
        x <- matrix(nrow = 1, ncol = 2)
        x[1,1] <- 0
        x[1,2] <- 0
        x <- as.data.frame(x)
        p <- ggplot(x, aes(x[,1],x[,2]))
        p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                       panel.background = element_blank(), axis.line = element_line(colour = "white"))
        p <- p + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
        p <- p + annotate("text", label = "Sorry, the number of individuals in the input dataset \n must be at least 3 \n \n Correlation matrix cannot be calculated!", 0,0, cex = 5.4)
        plot_together <- plot_grid(plot_grid(p, corrplot_gs, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
        print(plot_together)
      }
      
    } else if(input$case_control %in% "Subset by Cases") {
      if(nrow(cases_split_cases()) >= 3){
      cormatrix_gs <- cor(dat_epigen()[,predictions2()])
      cormatrix_input <- cor(cases_split_cases()[,predictions2()])
      corrplot_gs <- ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T)  + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
      corrplot_input <- ggcorrplot(cormatrix_input,type="lower", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle(paste0("Input Dataset ", paste(as.character(input$Case_Input3), "- Cases"))) + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
      legend_corr <- (get_legend(ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
      corr_together <- plot_grid(plot_grid(corrplot_input, corrplot_gs, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
      print(corr_together)
      } else { 
        cormatrix_gs <- cor(dat_epigen()[,predictions2()])
        corrplot_gs <- ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation",colors = rev(c("#6b81ff", "white", "#E46726")), lab = T)  + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
        legend_corr <- (get_legend(ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation",colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
        x <- matrix(nrow = 1, ncol = 2)
        x[1,1] <- 0
        x[1,2] <- 0
        x <- as.data.frame(x)
        p <- ggplot(x, aes(x[,1],x[,2]))
        p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                       panel.background = element_blank(), axis.line = element_line(colour = "white"))
        p <- p + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
        p <- p + annotate("text", label = "Sorry, the number of cases in the input dataset \n must be at least 3 \n \n Correlation matrix cannot be calculated!", 0,0, cex = 5.4)
        plot_together <- plot_grid(plot_grid(p, corrplot_gs, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
        print(plot_together)
        }
    } else if(input$case_control %in% "Subset by Controls") { 
      if(nrow(cases_split_controls()) >= 3) { 
      cormatrix_gs <- cor(dat_epigen()[,predictions2()])
      cormatrix_input <- cor(cases_split_controls()[,predictions2()])
      corrplot_gs <- ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T)  + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
      corrplot_input <- ggcorrplot(cormatrix_input,type="lower", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle(paste0("Input Dataset ", paste(as.character(input$Case_Input3), "- Controls"))) + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
      legend_corr <- (get_legend(ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
      corr_together <- plot_grid(plot_grid(corrplot_input, corrplot_gs, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
      print(corr_together)
      } else { 
        cormatrix_gs <- cor(dat_epigen()[,predictions2()])
        corrplot_gs <- ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation",colors = rev(c("#6b81ff", "white", "#E46726")), lab = T)  + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
        legend_corr <- (get_legend(ggcorrplot(cormatrix_gs,type="lower", legend.title = "Correlation",colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Generation Scotland Dataset") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
        x <- matrix(nrow = 1, ncol = 2)
        x[1,1] <- 0
        x[1,2] <- 0
        x <- as.data.frame(x)
        p <- ggplot(x, aes(x[,1],x[,2]))
        p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                       panel.background = element_blank(), axis.line = element_line(colour = "white"))
        p <- p + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
        p <- p + annotate("text", label = "Sorry, the number of controls in the input dataset \n must be at least 3 \n \n Correlation matrix cannot be calculated!", 0,0, cex = 5.4)
        plot_together <- plot_grid(plot_grid(p, corrplot_gs, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
        print(plot_together)
        
        }
      } else  {
        if(nrow(cases_split_cases()) >= 3 && nrow(cases_split_controls()) >= 3) {
      cormatrix_cases <- cor(cases_split_plot()[which(cases_split_plot()[,as.character(input$Case_Input3)] == 1),predictions2()])
      cormatrix_controls <- cor(cases_split_plot()[which(cases_split_plot()[,as.character(input$Case_Input3)] == 0),predictions2()])
      corrplot_cases <- ggcorrplot(cormatrix_cases,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T)  + ggplot2::ggtitle(paste0("Input Dataset ", paste(as.character(input$Case_Input3), "- Cases"))) + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
      corrplot_controls <- ggcorrplot(cormatrix_controls,type="lower", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle(paste0("Input Dataset ", paste(as.character(input$Case_Input3), "- Controls"))) + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
      legend_corr <- (get_legend(ggcorrplot(cormatrix_cases,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Input Dataset - Cases") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
      corr_together <- plot_grid(plot_grid(corrplot_cases, corrplot_controls, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
      print(corr_together)  
        } else if(!nrow(cases_split_cases()) >= 3 && nrow(cases_split_controls()) >= 3) {
          cormatrix_controls <- cor(cases_split_plot()[which(cases_split_plot()[,as.character(input$Case_Input3)] == 0),predictions2()])
          corrplot_controls <- ggcorrplot(cormatrix_controls,type="lower", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle(paste0("Input Dataset ", paste(as.character(input$Case_Input3), "- Controls"))) + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
          legend_corr <- (get_legend(ggcorrplot(cormatrix_controls,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Input Dataset - Controls") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
          x <- matrix(nrow = 1, ncol = 2)
          x[1,1] <- 0
          x[1,2] <- 0
          x <- as.data.frame(x)
          p <- ggplot(x, aes(x[,1],x[,2]))
          p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                         panel.background = element_blank(), axis.line = element_line(colour = "white"))
          p <- p + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
          p <- p + annotate("text", label = "Sorry, the number of cases in the input dataset \n must be at least 3 \n \n Correlation matrix cannot be calculated!", 0,0, cex = 5.4)
          plot_together <- plot_grid(plot_grid(p, corrplot_controls, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
          print(plot_together)
        } else if(nrow(cases_split_cases()) >= 3 && !nrow(cases_split_controls()) >= 3) { 
          cormatrix_cases <- cor(cases_split_plot()[which(cases_split_plot()[,as.character(input$Case_Input3)] == 1),predictions2()])
          corrplot_cases <- ggcorrplot(cormatrix_cases,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T)  + ggplot2::ggtitle(paste0("Input Dataset ", paste(as.character(input$Case_Input3), "- Cases"))) + ggplot2::theme(plot.title = element_text(hjust =0.5, size = 16.5)) + ggplot2::theme(legend.position = "none")
          legend_corr <- (get_legend(ggcorrplot(cormatrix_cases,type="lower", legend.title = "Correlation", colors = rev(c("#6b81ff", "white", "#E46726")), lab = T) + ggplot2::ggtitle("Input Dataset - Cases") + ggplot2::theme(plot.title = element_text(hjust =0.5)) + ggplot2::theme(legend.position = "right") + ggplot2::theme(legend.title.align = 0.5)))
          x <- matrix(nrow = 1, ncol = 2)
          x[1,1] <- 0
          x[1,2] <- 0
          x <- as.data.frame(x)
          p <- ggplot(x, aes(x[,1],x[,2]))
          p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                         panel.background = element_blank(), axis.line = element_line(colour = "white"))
          p <- p + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
          p <- p + annotate("text", label = "Sorry, the number of controls in the input dataset \n must be at least 3 \n \n Correlation matrix cannot be calculated!", 0,0, cex = 5.4)
          plot_together <- plot_grid(plot_grid(p, corrplot_cases, ncol = 2),legend_corr, rel_widths = c(3, 0.4))
          print(plot_together)
          
        } else { 
          x <- matrix(nrow = 1, ncol = 2)
          x[1,1] <- 0
          x[1,2] <- 0
          x <- as.data.frame(x)
          p <- ggplot(x, aes(x[,1],x[,2]))
          p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                         panel.background = element_blank(), axis.line = element_line(colour = "white"))
          p <- p + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
          p <- p + annotate("text", label = "Sorry, the number of cases in the input dataset \n must be at least 3 \n \n Correlation matrix cannot be calculated!", 0,0, cex = 5.4) 
          
          x1 <- matrix(nrow = 1, ncol = 2)
          x1[1,1] <- 0
          x1[1,2] <- 0
          x1 <- as.data.frame(x1)
          q <- ggplot(x1, aes(x1[,1],x1[,2]))
          q <- q + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                         panel.background = element_blank(), axis.line = element_line(colour = "white"))
          q <- q + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
          q <- q + annotate("text", label = "Sorry, the number of controls in the input dataset \n must be at least 3 \n \n Correlation matrix cannot be calculated!", 0,0, cex = 5.4)
          plot_together <- plot_grid(plot_grid(p, q, ncol = 2))
          print(plot_together)
          
          }
      }
    }
      
    
  })
  
  output$corrplot2 <- renderPlot({
    if(is.null(input$file1)){
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0,"Please upload dataset to launch the app!", cex = 1.7) 
      
    } else if (nrow(dat_epigen()) <= 10 ) {
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0,"Sorry, the number of reference individuals must be above 10!", cex = 1.7) 
    } else { 
      ## Correlation between Predictors and Respective Traits in GS 
      
      corr<-merge(tmp, out, by.x = "Basename", by.y="ID")
      corr = corr[which(corr$Basename %in% dat_epigen()$ID),]
      corrmatrix <- matrix(nrow = 7, ncol = 2)
      corrmatrix[1,2] <- cor(corr$age, corr$`Biological Age`, use = "complete.obs")
      corrmatrix[2,2] <- cor(corr$Alcohol[corr$usual %in% 2], corr$units[corr$usual %in% 2], use = "complete.obs")
      corrmatrix[3,2] <- cor(corr$`Body Fat %`, corr$body_fat, use = "complete.obs")
      corrmatrix[4,2] <- cor(corr$`Body Mass Index`, corr$bmi, use = "complete.obs")
      corrmatrix[5,2] <- cor(corr$`HDL Cholesterol`, corr$HDL_cholesterol, use = "complete.obs")
      corrmatrix[6,2] <-cor(corr$Smoking, corr$pack_years, use = "complete.obs")
      corrmatrix[7,2] <-cor(corr$`Waist:Hip Ratio`, corr$whr, use = "complete.obs")
      
      corrmatrix[1,1] <- "Age"
      corrmatrix[2,1] <- "Alcohol"
      corrmatrix[3,1] <- "Body Fat %"
      corrmatrix[4,1] <- "Body Mass Index"
      corrmatrix[5,1] <- "HDL Cholesterol"
      corrmatrix[6,1] <- "Smoking"
      corrmatrix[7,1] <- "Waist:Hip Ratio"
      
      corrmatrix <- as.data.frame(corrmatrix)
      names(corrmatrix) <- c("Trait", "r")
      corrmatrix$r <- as.numeric(as.character(corrmatrix$r))
      row.names(corrmatrix) <- corrmatrix$Trait
      corrmatrix$Trait <- NULL
      
      
      
      corr_gs <- ggcorrplot(as.matrix(corrmatrix), lab = TRUE, hc.order = F) + ggplot2::theme(legend.position = "bottom" ,legend.key.width = unit(3, "cm")) + 
        scale_fill_gradient2(breaks = c(0,0.2,0.4,0.6,0.8,1.0), limit=c(0, 1), low = "white", mid = "cyan",high = "blue") + 
        scale_y_discrete(labels =expression(paste(italic("r"), " coefficient", sep = " "))) + ggtitle("Correlation between Traits and their Epigenetic Predictors in Generation Scotland") + theme(plot.title = element_text(hjust = 0.5))
      print(corr_gs)
    }   
  })
  
  thedata2 <- reactive({
    n =nrow(df2())
    withProgress(message = "Generating Input Percentiles for Download", value = 0, {
      for(i in 1:n){ incProgress(1/n, detail = paste(round(((i/n)*100)), "% complete"))
      } 
    reference_data = cbind(df2()$ID, as.data.frame(apply(df2()[,as.character(predictors)], 2, function(x) percent_rank(x)*100)))
    })
    names(reference_data)[1] <- "ID"
    reference_data <- merge(reference_data, df2()[,c("ID","Biological Age")], by = "ID")
    reference_data <- cbind(reference_data$ID,round(reference_data[,2:ncol(reference_data)], 2))
    names(reference_data)[1] <- "ID"
    return(reference_data)
  })
  
  
  
  output$downloadinput <- renderUI({
    if(!is.null(input$file1)) {
      downloadButton("download1", "Download Percentile Ranks for Individuals within Input Dataset")
    }
  })
  
 

  output$dto1 <- renderDataTable({thedata2()})
  output$download1 <- downloadHandler(
    filename = function(){paste("Percentiles_Comparing_Input_including_age", ".csv", sep= "")}, 
    content = function(file1){
      fwrite(thedata2(), file1, row.names = F,dateTimeAs = "write.csv")
    }
  )
  
  observeEvent(input$timeOut, { 
    print(paste0("Session (", session$token, ") timed out at: ", Sys.time()))
    showModal(modalDialog(
      title = "Timeout",
      paste("Session Time Due To", input$timeOut, "Inactivity -", Sys.time()),
      footer = NULL
    ))
    session$close()
    stopApp()
  })
  
  
  
}


shinyApp(ui=ui,server=server)
