## 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("Epigenetic Age", "Lifestyle or Biochemical" = c("Alcohol", 
                                                                                                           "Body Fat %",
                                                                                                           "Body Mass Index",
                                                                                                           "HDL Cholesterol",
                                                                                                           "Smoking",
                                                                                                           "Waist:Hip Ratio"), "Proteins" = c(
                                                                                                             "ADAMTS",
                                                                                                             "Adiponectin",
                                                                                                             "Afamin",
                                                                                                             "Alpha-L-iduronidase",
                                                                                                             "Aminoacylase-1",
                                                                                                             "B2-microglobulin",
                                                                                                             "BCAM",
                                                                                                             "BMP-1",
                                                                                                             "CCL11",
                                                                                                             "CCL17",
                                                                                                             "CCL18",
                                                                                                             "CCL21",
                                                                                                             "CCL22",
                                                                                                             "CCL25 C-C",
                                                                                                             "CD163",
                                                                                                             "CD209 antigen",
                                                                                                             "CD48 antigen",
                                                                                                             "CD6",
                                                                                                             "CDL5",
                                                                                                             "CHIT-1",
                                                                                                             "CLEC11A e1",
                                                                                                             "CLEC11A e2",
                                                                                                             "Coagulation factor VII",
                                                                                                             "Complement C4",
                                                                                                             "Complement C5a",
                                                                                                             "Complement c9",
                                                                                                             "Contactin-4",
                                                                                                             "CRP",
                                                                                                             "CRTAM",
                                                                                                             "CXCL10",
                                                                                                             "CXCL10 soma",
                                                                                                             "CXCL11",
                                                                                                             "CXCL11 soma",
                                                                                                             "CXCL9",
                                                                                                             "E-selectin",
                                                                                                             "Ectodysplasin-A",
                                                                                                             "EN.RAGE",
                                                                                                             "ENPP7",
                                                                                                             "ESM-1",
                                                                                                             "EZR",
                                                                                                             "FAP",
                                                                                                             "FCER2",
                                                                                                             "FCGR3A",
                                                                                                             "FcRL2",
                                                                                                             "FGF.21",
                                                                                                             "G.CSF",
                                                                                                             "Galectin-4",
                                                                                                             "GDF.8",
                                                                                                             "GHR",
                                                                                                             "GPIba",
                                                                                                             "Granulysin",
                                                                                                             "Granzyme A",
                                                                                                             "GZMA",
                                                                                                             "HCII",
                                                                                                             "HGF",
                                                                                                             "HGFA",
                                                                                                             "HGFI",
                                                                                                             "ICAM5",
                                                                                                             "IGFBP-1",
                                                                                                             "IGFBP-4",
                                                                                                             "Insulin receptor",
                                                                                                             "Interleukin-19",
                                                                                                             "L-selectin",
                                                                                                             "LFT",
                                                                                                             "LGALS3BP",
                                                                                                             "LY9",
                                                                                                             "Lymphotoxin abeta",
                                                                                                             "Lysozyme C",
                                                                                                             "MIA",
                                                                                                             "MMP-1",
                                                                                                             "MMP-12",
                                                                                                             "MMP-2",
                                                                                                             "MMP-9",
                                                                                                             "MMP.1",
                                                                                                             "MRC2",
                                                                                                             "Myeloperoxidase",
                                                                                                             "N.CDase",
                                                                                                             "NCAM-120",
                                                                                                             "NEP",
                                                                                                             "NMNAT1",
                                                                                                             "NOTCH1",
                                                                                                             "NRTK3",
                                                                                                             "NTRK3",
                                                                                                             "OSM",
                                                                                                             "Osteomodulin",
                                                                                                             "PAPP-A",
                                                                                                             "PIGR",
                                                                                                             "RARRES2",
                                                                                                             "Resistin",
                                                                                                             "S100-A9",
                                                                                                             "Semaphorin-3E",
                                                                                                             "SERPIN A3",
                                                                                                             "SHBG",
                                                                                                             "SIGLEC1",
                                                                                                             "SKR3",
                                                                                                             "SLITRK5",
                                                                                                             "SMPD1",
                                                                                                             "Stanniocalcin-1",
                                                                                                             "Testican-2",
                                                                                                             "TGF.alpha",
                                                                                                             "THBS2",
                                                                                                             "TNFRSF17",
                                                                                                             "TNFRSF1B",
                                                                                                             "TPO",
                                                                                                             "Trypsin-2",
                                                                                                             "Tryptase beta-2",
                                                                                                             "VCAM1",
                                                                                                             "VEGFA",
                                                                                                             "WFIKKN2")),
                                   multiple = F, selected = NULL),
                       uiOutput("var_ui"),
                       uiOutput("df_cases_choice"), 
                        actionButton("press1", "Press Here For Information on Panel and Useful Links"),
      br(),
      br(),
      actionButton("version", "Please Cite The Following Papers - Press Here"),
      br(),
      br(),
      actionButton("format", "Press Here For File Formats and Useful Links"),
      br(),
      br()),

      
      conditionalPanel(condition = "input.tabs == 2", 
                       actionButton("press2", "Press Here For Information on Panel and Useful Links"),
                       br(), 
                       br(),
                       selectizeInput("var1", 
                                   label = "Choose Epigenetic Predictors to Display",
                                   choices = list("Lifestyle or Biochemical" = c("Alcohol","Body Fat %",
                                                                                    "Body Mass Index",
                                                                                    "HDL Cholesterol",
                                                                                    "Smoking",
                                                                                    "Waist:Hip Ratio"), "Proteins"= 
                                                    
                                                    c("ADAMTS",
                                                      "Adiponectin",
                                                      "Afamin",
                                                      "Alpha-L-iduronidase",
                                                      "Aminoacylase-1",
                                                      "B2-microglobulin",
                                                      "BCAM",
                                                      "BMP-1",
                                                      "CCL11",
                                                      "CCL17",
                                                      "CCL18",
                                                      "CCL21",
                                                      "CCL22",
                                                      "CCL25 C-C",
                                                      "CD163",
                                                      "CD209 antigen",
                                                      "CD48 antigen",
                                                      "CD6",
                                                      "CDL5",
                                                      "CHIT-1",
                                                      "CLEC11A e1",
                                                      "CLEC11A e2",
                                                      "Coagulation factor VII",
                                                      "Complement C4",
                                                      "Complement C5a",
                                                      "Complement c9",
                                                      "Contactin-4",
                                                      "CRP",
                                                      "CRTAM",
                                                      "CXCL10",
                                                      "CXCL10 soma",
                                                      "CXCL11",
                                                      "CXCL11 soma",
                                                      "CXCL9",
                                                      "E-selectin",
                                                      "Ectodysplasin-A",
                                                      "EN.RAGE",
                                                      "ENPP7",
                                                      "ESM-1",
                                                      "EZR",
                                                      "FAP",
                                                      "FCER2",
                                                      "FCGR3A",
                                                      "FcRL2",
                                                      "FGF.21",
                                                      "G.CSF",
                                                      "Galectin-4",
                                                      "GDF.8",
                                                      "GHR",
                                                      "GPIba",
                                                      "Granulysin",
                                                      "Granzyme A",
                                                      "GZMA",
                                                      "HCII",
                                                      "HGF",
                                                      "HGFA",
                                                      "HGFI",
                                                      "ICAM5",
                                                      "IGFBP-1",
                                                      "IGFBP-4",
                                                      "Insulin receptor",
                                                      "Interleukin-19",
                                                      "L-selectin",
                                                      "LFT",
                                                      "LGALS3BP",
                                                      "LY9",
                                                      "Lymphotoxin abeta",
                                                      "Lysozyme C",
                                                      "MIA",
                                                      "MMP-1",
                                                      "MMP-12",
                                                      "MMP-2",
                                                      "MMP-9",
                                                      "MMP.1",
                                                      "MRC2",
                                                      "Myeloperoxidase",
                                                      "N.CDase",
                                                      "NCAM-120",
                                                      "NEP",
                                                      "NMNAT1",
                                                      "NOTCH1",
                                                      "NRTK3",
                                                      "NTRK3",
                                                      "OSM",
                                                      "Osteomodulin",
                                                      "PAPP-A",
                                                      "PIGR",
                                                      "RARRES2",
                                                      "Resistin",
                                                      "S100-A9",
                                                      "Semaphorin-3E",
                                                      "SERPIN A3",
                                                      "SHBG",
                                                      "SIGLEC1",
                                                      "SKR3",
                                                      "SLITRK5",
                                                      "SMPD1",
                                                      "Stanniocalcin-1",
                                                      "Testican-2",
                                                      "TGF.alpha",
                                                      "THBS2",
                                                      "TNFRSF17",
                                                      "TNFRSF1B",
                                                      "TPO",
                                                      "Trypsin-2",
                                                      "Tryptase beta-2",
                                                      "VCAM1",
                                                      "VEGFA",
                                                      "WFIKKN2")),
                                   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", 
                      
                       actionButton("press3", "Press Here For Information on Panel and Useful Links"),
                    br(),
                    br(),
                       selectizeInput("var2", 
                                      label = "Choose Epigenetic Predictors to Display",
                                      choices = list("Lifestyle or Biochemical" = c("Alcohol","Body Fat %",
                                                                                    "Body Mass Index",
                                                                                    "HDL Cholesterol",
                                                                                    "Smoking",
                                                                                    "Waist:Hip Ratio"), "Proteins"= 
                                                       
                                                       c("ADAMTS",
                                                         "Adiponectin",
                                                         "Afamin",
                                                         "Alpha-L-iduronidase",
                                                         "Aminoacylase-1",
                                                         "B2-microglobulin",
                                                         "BCAM",
                                                         "BMP-1",
                                                         "CCL11",
                                                         "CCL17",
                                                         "CCL18",
                                                         "CCL21",
                                                         "CCL22",
                                                         "CCL25 C-C",
                                                         "CD163",
                                                         "CD209 antigen",
                                                         "CD48 antigen",
                                                         "CD6",
                                                         "CDL5",
                                                         "CHIT-1",
                                                         "CLEC11A e1",
                                                         "CLEC11A e2",
                                                         "Coagulation factor VII",
                                                         "Complement C4",
                                                         "Complement C5a",
                                                         "Complement c9",
                                                         "Contactin-4",
                                                         "CRP",
                                                         "CRTAM",
                                                         "CXCL10",
                                                         "CXCL10 soma",
                                                         "CXCL11",
                                                         "CXCL11 soma",
                                                         "CXCL9",
                                                         "E-selectin",
                                                         "Ectodysplasin-A",
                                                         "EN.RAGE",
                                                         "ENPP7",
                                                         "ESM-1",
                                                         "EZR",
                                                         "FAP",
                                                         "FCER2",
                                                         "FCGR3A",
                                                         "FcRL2",
                                                         "FGF.21",
                                                         "G.CSF",
                                                         "Galectin-4",
                                                         "GDF.8",
                                                         "GHR",
                                                         "GPIba",
                                                         "Granulysin",
                                                         "Granzyme A",
                                                         "GZMA",
                                                         "HCII",
                                                         "HGF",
                                                         "HGFA",
                                                         "HGFI",
                                                         "ICAM5",
                                                         "IGFBP-1",
                                                         "IGFBP-4",
                                                         "Insulin receptor",
                                                         "Interleukin-19",
                                                         "L-selectin",
                                                         "LFT",
                                                         "LGALS3BP",
                                                         "LY9",
                                                         "Lymphotoxin abeta",
                                                         "Lysozyme C",
                                                         "MIA",
                                                         "MMP-1",
                                                         "MMP-12",
                                                         "MMP-2",
                                                         "MMP-9",
                                                         "MMP.1",
                                                         "MRC2",
                                                         "Myeloperoxidase",
                                                         "N.CDase",
                                                         "NCAM-120",
                                                         "NEP",
                                                         "NMNAT1",
                                                         "NOTCH1",
                                                         "NRTK3",
                                                         "NTRK3",
                                                         "OSM",
                                                         "Osteomodulin",
                                                         "PAPP-A",
                                                         "PIGR",
                                                         "RARRES2",
                                                         "Resistin",
                                                         "S100-A9",
                                                         "Semaphorin-3E",
                                                         "SERPIN A3",
                                                         "SHBG",
                                                         "SIGLEC1",
                                                         "SKR3",
                                                         "SLITRK5",
                                                         "SMPD1",
                                                         "Stanniocalcin-1",
                                                         "Testican-2",
                                                         "TGF.alpha",
                                                         "THBS2",
                                                         "TNFRSF17",
                                                         "TNFRSF1B",
                                                         "TPO",
                                                         "Trypsin-2",
                                                         "Tryptase beta-2",
                                                         "VCAM1",
                                                         "VEGFA",
                                                         "WFIKKN2")),
                                      multiple = T, options = list(maxItems = 10), selected = NULL),
                       
                uiOutput("radioButtons_case_control"),
                uiOutput("df_cases_choice3")),
      
      conditionalPanel(condition = "input.tabs == 4", 
                       actionButton("press4", "Press Here For Information on Panel and Useful Links"),
                       br(),
                       br()
                       
      ), 
      
      
    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) {
  options(shiny.maxRequestSize = 1000*1024^2)
  
  
    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:ncol(df1())], 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)
    }) 
  
  
  url <- a("'MethylDetectR - Calculate Your Scores'", href="https://shiny.igmm.ed.ac.uk/Calculate_Your_Scores/", target = "_blank")
  url_demo <- a("'MethylDetectR' Demo", href="https://shiny.igmm.ed.ac.uk/MethylDetectR_Demo/", target = "_blank")
  url2 <- a("'MethylDetectR' Website", href="https://www.ed.ac.uk/centre-genomic-medicine/research-groups/marioni-group/methyldetectr", target = "_blank")
  url3 <- a("Example Input and Output Files in Zenodo", href="https://zenodo.org/record/5869576#.YeaBaf7P3IU", target = "_blank")
  url4 <- a("'MethylDetectR' Paper", href="https://wellcomeopenresearch.org/articles/5-283", target = "_blank")
  
  p1 <- a("McCartney & Hillary - Lifestyle/Biochemical", href="https://genomebiology.biomedcentral.com/articles/10.1186/s13059-018-1514-1", target = "_blank")
  p2 <- a("Gadd - Proteins", href="https://elifesciences.org/articles/71802", target = "_blank")
  p3 <- a("Zhang - Epigenetic Age", href="https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-019-0667-1", target = "_blank")
  p4 <- a("Bernabeu - Epigenetic Age", href="https://www.biorxiv.org/content/10.1101/2022.09.08.507115v1", target = "_blank")
  p5 <- a("Hillary & Marioni - MethylDetectR", href="https://wellcomeopenresearch.org/articles/5-283", target = "_blank")
  
  
  observeEvent(input$format,{ 
    # show pop-up ...
    showModal(modalDialog(
      title = "How to Format Your Files",
      div("Example input files are displayed in the Zenodo repository (link at bottom).", style="font-size:90%"),
      br(),
      div("The file containting DNAm-based scores can be generated using the 'MethylDetectR - Calculate Your Scores' application, or an R script provided in Zenodo.", style="font-size:90%"),
      br(),
      div("This is a .csv file with the first column containing IDs for individuals as in the methylation file. The next two columns: 'Sex' and 'True Age' may contain missing values if the user does not wish to upload these data. The remaining columns detail DNAm-based scores.", style="font-size:90%"),
      br(),
      div("If 'True Age' is provided, we will use this to subset the data by age range in the age slider. If these data are missing, predicted age (Bernabeu predictor) will be used.", style="font-size:90%"),
      br(),
      div("An optional case vs. control file as a .csv file should contain an 'ID' column with IDs for individuals as in the methylation file and any binary traits of interest for the remainder of the columns with controls coded as '0' and cases as '1'.", style="font-size:90%"),
      br(),
      div("Applications will time out after three minutes of inactivity. Important information on data privacy and protection can be accessed in our website using the link below.", style="font-size:90%"),
      br(),
      div("Please use the following links to access the other parts of our platform.", style="font-size:90%"),
      br(),
      tagList(url),
      br(),
      tagList(url_demo),
      br(),
      tagList(url2),
      br(),
      tagList(url3),
      br(),
      tagList(url4),
      easyClose = TRUE,
      footer = NULL))
    
  })
  
  observeEvent(input$press1,{ 
    # show pop-up ...
    showModal(modalDialog(
      title = "Information on Panel 1",
      div("In this panel, we can look at scores for individuals in your dataset. For example, we estimate Person X is 60 years old. We also estimate that they score in the 10th percentile for alcohol consumption in your dataset. You can also see how scores vary in cases and controls for binary variables of interest in your dataset.", style="font-size:90%"),
      br(),
      br(),
      div("This suggests that Person X is a light drinker relative to others in the dataset. However, these predictions are not always accurate at an individual level, they work well at a population level. Please refer to Panel 4 to see how well the predictors perform in Generation Scotland according to different age ranges and sex.", style="font-size:90%"),
      br(),
      br(),
      div("These scores might help us gain information about someone's characteristics and health profile but they need to be more accurate at an individual level to be used clinically.", style="font-size:90%"),
      br(),
      div("Applications will time out after three minutes of inactivity. Important information on data privacy and protection can be accessed in our website using the link below.", style="font-size:90%"),
      br(),
      div("Please use the following links to access the other parts of our platform.", style="font-size:90%"),
      br(),
      tagList(url),
      br(),
      tagList(url_demo),
      br(),
      tagList(url2),
      br(),
      tagList(url3),
      br(),
      tagList(url4),
      
      easyClose = TRUE,
      footer = NULL))
    
  })
  
  observeEvent(input$version,{ 
    # show pop-up ...
    showModal(modalDialog(
      title = "Please Cite These Papers If You Use This App",
      div(HTML(paste0("1. Lifestyle and Biochemical Trait Predictors: McCartney and Hillary ",  em("et al."), " (2018). PMID: 30257690")),style="font-size:90%"),
      tagList(p1),
      br(),
      br(),
      div(HTML(paste0("2. Protein Trait Predictors: Gadd ",  em("et al."), " (2022). PMID: 35023833")),style="font-size:90%"),
      tagList(p2),
      br(),
      br(),
      div(HTML(paste0("3. Epigenetic Age: Zhang ",  em("et al."), " (2019). PMID: 30257690")),style="font-size:90%"),
      tagList(p3),
      br(),
      br(),
      div(HTML(paste0("4. Epigenetic Age: Bernabeu ",  em("et al."), " (2022). PMID: pending")),style="font-size:90%"),
      tagList(p4),
      br(),
      br(),
      div(HTML(paste0("5. MethylDetectR: Hillary and Marioni (2021). PMID: 33969230")),style="font-size:90%"),
      tagList(p5),
      br(),
      br(),
      easyClose = TRUE,
      footer = NULL))
    
  })
  
  
  observeEvent(input$press1,{ 
    # show pop-up ...
    showModal(modalDialog(
      title = "Information on Panel 1",
      div("To develop the predictors in this app, we can for instance look at patterns of DNA methylation in smokers and non-smokers. We can identify general differences between these groups and make algorithms using these differences to estimate smoking status of other individuals.", style="font-size:90%"),
      br(),
      br(),
      div("In this panel, we can look at scores for individuals in your dataset. For example, we estimate Person X is 60 years old. We also estimate that they score in the 10th percentile for alcohol consumption in your dataset. You can also see how scores vary in cases and controls for binary variables of interest in your dataset.", style="font-size:90%"),
      br(),
      br(),
      div("This suggests that Person X is a light drinker relative to others in the dataset. However, these predictions are not always accurate at an individual level, they work well at a population level.", style="font-size:90%"),
      br(),
      br(),
      div("These scores might help us gain information about someone's characteristics and health profile but they need to be more accurate at an individual level to be used clinically.", style="font-size:90%"),
      br(),
      div("Applications will time out after three minutes of inactivity. Important information on data privacy and protection can be accessed in our website using the link below.", style="font-size:90%"),
      br(),
      div("Please use the following links to access the other parts of our platform.", style="font-size:90%"),
      br(),
      tagList(url),
      br(),
      tagList(url_demo),
      br(),
      tagList(url2),
      br(),
      tagList(url3),
      br(),
      tagList(url4),
      
      easyClose = TRUE,
      footer = NULL))
    
  })
  
  
  observeEvent(input$press3,{ 
    # show pop-up ...
    showModal(modalDialog(
      title = "Information on Panel 3",
      div("Here, we can see correlations between methylation-based scores, side-by-side, for different groups:", style="font-size:90%"),
      br(),
      div("1). Scores in your input dataset versus scores in Generation Scotland.", style="font-size:87%"),
       div("2). Cases in your dataset, for a variable of interest, and Generation Scotland participants.", style="font-size:87%"),
      div("3). Controls in your dataset, for a variable of interest, and Generation Scotland participants.", style="font-size:87%"),
      div("4). The cases versus the controls.", style="font-size:87%"),
      br(),
      br(),
      div("As an example, we might expect to see high correlations between some predictors such as those for body fat % and body mass index. All plots can be updated by age range and/or sex.", style="font-size:89%"),
      br(),
      br(),
      div("If there are unexpected differences in correlations between groups, it may mean that the performance of  predictors also differ between these groups or that the groups have different characteristics.", style="font-size:89%"),
      br(),
      br(),
      div("Applications will time out after three minutes of inactivity. Important information on data privacy and protection can be accessed in our website using the link below.", style="font-size:90%"),
      br(),
      div("Please use the following links to access the other parts of our platform.", style="font-size:90%"),
      br(),
      tagList(url),
      br(),
      tagList(url_demo),
      br(),
      tagList(url2),
      br(),
      tagList(url3),
      br(),
      tagList(url4),
      
      easyClose = TRUE,
      footer = NULL))
    
  })
  
  observeEvent(input$press4,{ 
    # show pop-up ...
    showModal(modalDialog(
      title = "Information on Panel 4",
      div("This panel shows how well methylation-based scores for various traits correlate with true values in Generation Scotland.", style="font-size:90%"),
      br(),
      br(),
      div("For example, methylation-based scores for age i.e. predicted age correlate 99% with true age in Generation Scotland when no subset is made by age range and/or sex. Here, the Zhang predictor for age is used and not the Bernabeu predictor as the latter used the test dataset shown in Panel 4 in their training sample whereas the former did not. This does not affect any other features or panels in the MethylDetectR platform.", style="font-size:90%"),
      br(),
      div("You can see how these correlations or 'performance indicators' change when you subset by age range and/or sex.", style="font-size:90%"),
      br(),
      div("This is important because if you see that the correlations become poor in a given subset, then this may also apply to your own dataset using the same subset criteria but this depends on how similar your dataset is to Generation Scotland.", style="font-size:90%"),
      br(),  
      div("Applications will time out after three minutes of inactivity. Important information on data privacy and protection can be accessed in our website using the link below.", style="font-size:90%"),
      br(),
      div("Please use the following links to access the other parts of our platform.", style="font-size:90%"),
      br(),
      tagList(url),
      br(),
      tagList(url_demo),
      br(),
      tagList(url2),
      br(),
      tagList(url3),
      br(),
      tagList(url4),
      
      easyClose = TRUE,
      footer = NULL))
    
  })
  
  observeEvent(input$press2,{ 
    # show pop-up ...
    showModal(modalDialog(
      title = "Information on Panel 2",
      div("This panel lets you see where a particular individual lies along your input dataset with respect to their methylation-based scores.", style="font-size:90%"),
      br(),
      br(),
      div("As an example, Person X may be a non-smoker and light-drinker. Their methylation-based scores put them in the 10th percentile for smoking (maybe passive smoking or imperfect performance of predictor) and 5th percentile for alcohol. The scores show that they have a favourable lifestyle profile relative to others in the dataset.", style="font-size:86%"),
      br(),
      br(),
      div("You can also see the average (median) percentile ranks for cases in your dataset. As an example, your variable of interest is people with lung disease vs. people without lung disease. For smoking scores, the median percentile rank of cases is high, it is the 75th percentile. An interquartile range of the 70th-80th percentile tells us that 50% of all cases lie between the 70th and 80th percentile, covering a lot of higher smoking scores compared to controls.", style="font-size:87%"),
      br(),
      br(),
      div("Applications will time out after three minutes of inactivity. Important information on data privacy and protection can be accessed in our website using the link below.", style="font-size:90%"),
      br(),
      div("Please use the following links to access the other parts of our platform.", style="font-size:90%"),
      br(),
      tagList(url),
      br(),
      tagList(url_demo),
      br(),
      tagList(url2),
      br(),
      tagList(url3),
      br(),
      tagList(url4),
      
      easyClose = TRUE,
      footer = NULL))
    
  })
  
  
  output$message1 <- renderUI(HTML(paste(
    "Important Note: If you have not already generated epigenetic scores, please do so using the 'Calculate Your Scores app' and refer to this app for further details"
  )))
  
  output$message2 <- renderUI(HTML(paste(
    "Important Note: If you have not already generated epigenetic scores, please do so using the 'Calculate Your Scores app' and refer to this app for further details"
  )))
  
  output$message3 <- renderUI(HTML(paste(
    "Important Note: If you have not already generated epigenetic scores, please do so using the 'Calculate Your Scores app' and refer to this app for further details"
  )))
  
  output$message4 <- renderUI(HTML(paste(
    "Important Note: If you have not already generated epigenetic scores, please do so using the 'Calculate Your Scores app' and refer to this app for further details"
  )))
  
  
  
  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),2)
    max.epigen <- round(max(epigen$`True Age`, na.rm = T),2)
    
    min.uploaded = round(min(c(df1()$'True Age', df1()$'Epigenetic Age (Bernabeu)'), na.rm = T),2)
    max.uploaded = round(max(c(df1()$'True Age', df1()$'Epigenetic Age (Bernabeu)'), na.rm = T),2)
    
    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()$'Epigenetic Age (Bernabeu)' >= input$ageInput[1] & df1()$'Epigenetic Age (Bernabeu)' <= input$ageInput[2]),]  
      return(test1) }  
      
      else  {
        test2 <- df1()[which(df1()$'Epigenetic Age (Bernabeu)' >= input$ageInput[1] & df1()$'Epigenetic Age (Bernabeu)' <= 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$'Epigenetic Age (Bernabeu)' >= input$ageInput[1] & has_biological_age$'Epigenetic Age (Bernabeu)' <= 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$'Epigenetic Age (Bernabeu)' >= input$ageInput[1] & has_biological_age$'Epigenetic Age (Bernabeu)' <= 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",
                   "ADAMTS"="ADAMTS",
                   "Adiponectin"="Adiponectin",
                   "Afamin"="Afamin",
                   "Alpha-L-iduronidase"="Alpha-L-iduronidase",
                   "Aminoacylase-1"="Aminoacylase-1",
                   "B2-microglobulin"="B2-microglobulin",
                   "BCAM"="BCAM",
                   "BMP-1"="BMP-1",
                   "CCL11"="CCL11",
                   "CCL17"="CCL17",
                   "CCL18"="CCL18",
                   "CCL21"="CCL21",
                   "CCL22"="CCL22",
                   "CCL25 C-C"="CCL25 C-C",
                   "CD163"="CD163",
                   "CD209 antigen"="CD209 antigen",
                   "CD48 antigen"="CD48 antigen",
                   "CD6"="CD6",
                   "CDL5"="CDL5",
                   "CHIT-1"="CHIT-1",
                   "CLEC11A e1"="CLEC11A e1",
                   "CLEC11A e2"="CLEC11A e2",
                   "Coagulation factor VII"="Coagulation factor VII",
                   "Complement C4"="Complement C4",
                   "Complement C5a"="Complement C5a",
                   "Complement c9"="Complement c9",
                   "Contactin-4"="Contactin-4",
                   "CRP"="CRP",
                   "CRTAM"="CRTAM",
                   "CXCL10"="CXCL10",
                   "CXCL10 soma"="CXCL10 soma",
                   "CXCL11"="CXCL11",
                   "CXCL11 soma"="CXCL11 soma",
                   "CXCL9"="CXCL9",
                   "E-selectin"="E-selectin",
                   "Ectodysplasin-A"="Ectodysplasin-A",
                   "EN.RAGE"="EN.RAGE",
                   "ENPP7"="ENPP7",
                   "ESM-1"="ESM-1",
                   "EZR"="EZR",
                   "FAP"="FAP",
                   "FCER2"="FCER2",
                   "FCGR3A"="FCGR3A",
                   "FcRL2"="FcRL2",
                   "FGF.21"="FGF.21",
                   "G.CSF"="G.CSF",
                   "Galectin-4"="Galectin-4",
                   "GDF.8"="GDF.8",
                   "GHR"="GHR",
                   "GPIba"="GPIba",
                   "Granulysin"="Granulysin",
                   "Granzyme A"="Granzyme A",
                   "GZMA"="GZMA",
                   "HCII"="HCII",
                   "HGF"="HGF",
                   "HGFA"="HGFA",
                   "HGFI"="HGFI",
                   "ICAM5"="ICAM5",
                   "IGFBP-1"="IGFBP-1",
                   "IGFBP-4"="IGFBP-4",
                   "Insulin receptor"="Insulin receptor",
                   "Interleukin-19"="Interleukin-19",
                   "L-selectin"="L-selectin",
                   "LFT"="LFT",
                   "LGALS3BP"="LGALS3BP",
                   "LY9"="LY9",
                   "Lymphotoxin abeta"="Lymphotoxin abeta",
                   "Lysozyme C"="Lysozyme C",
                   "MIA"="MIA",
                   "MMP-1"="MMP-1",
                   "MMP-12"="MMP-12",
                   "MMP-2"="MMP-2",
                   "MMP-9"="MMP-9",
                   "MMP.1"="MMP.1",
                   "MRC2"="MRC2",
                   "Myeloperoxidase"="Myeloperoxidase",
                   "N.CDase"="N.CDase",
                   "NCAM-120"="NCAM-120",
                   "NEP"="NEP",
                   "NMNAT1"="NMNAT1",
                   "NOTCH1"="NOTCH1",
                   "NRTK3"="NRTK3",
                   "NTRK3"="NTRK3",
                   "OSM"="OSM",
                   "Osteomodulin"="Osteomodulin",
                   "PAPP-A"="PAPP-A",
                   "PIGR"="PIGR",
                   "RARRES2"="RARRES2",
                   "Resistin"="Resistin",
                   "S100-A9"="S100-A9",
                   "Semaphorin-3E"="Semaphorin-3E",
                   "SERPIN A3"="SERPIN A3",
                   "SHBG"="SHBG",
                   "SIGLEC1"="SIGLEC1",
                   "SKR3"="SKR3",
                   "SLITRK5"="SLITRK5",
                   "SMPD1"="SMPD1",
                   "Stanniocalcin-1"="Stanniocalcin-1",
                   "Testican-2"="Testican-2",
                   "TGF.alpha"="TGF.alpha",
                   "THBS2"="THBS2",
                   "TNFRSF17"="TNFRSF17",
                   "TNFRSF1B"="TNFRSF1B",
                   "TPO"="TPO",
                   "Trypsin-2"="Trypsin-2",
                   "Tryptase beta-2"="Tryptase beta-2",
                   "VCAM1"="VCAM1",
                   "VEGFA"="VEGFA",
                   "WFIKKN2"="WFIKKN2")
    

    if(input$var %in% "Epigenetic Age"){ 
      
      plot(0,0, axes = F, xlab = "", ylab= "", pch = NA)
      text(0,0,paste("Epigenetic Age (Bernabeu) estimate:", round(df1()$'Epigenetic Age (Bernabeu)'[id1()],1), "years*", sep = " "), cex = 1.7)
      text(0,0,paste("\n", "\n", "\n", "Epigenetic Age (Zhang) estimate:", round(df1()$'Epigenetic Age (Zhang)'[id1()],1), "years*", sep = " "), cex = 1.7)
      text(0,0, paste("\n","\n", "\n", "\n", "\n", "\n", "\n","\n","\n", "\n", "\n", "\n", "\n","\n", "\n","\n", "*All methylation-based predictors in 'MethylDetectR' can make inaccurate", "\n", "predictions at an individual level, limiting clinical utility. However, they can work well at the", "\n", "population level and will improve with refined prediction methods and larger-scale studies."), cex = 0.96)  
      
             
        
    } 
    
    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 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 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 at least two 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$`Epigenetic Age (Zhang)`, 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","Epigenetic Age (Zhang)", "Epigenetic Age (Bernabeu)")], 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)
