Shiny.Rmd

---
title: "second"
author: "Owen Wen"
date: "2023-05-21"
output: html_document
---

```{r}
# List of packages to install
pkg_list <- c("dplyr", "stringr", "caret", "plotly", 
              "shiny", "bslib", "shinydashboard", "kernlab")

# Bioconductor packages
bioconductor_packages <- c("GEOquery", "Biobase", "limma")

# Function to install packages if not already installed
install_packages <- function(pkg_list){
  for(pkg in pkg_list){
    if(!pkg %in% installed.packages()){
      install.packages(pkg, dependencies = TRUE)
      print(paste(pkg, "has been installed successfully!"))
    } else {
      print(paste(pkg, "is already installed."))
    }
  }
}

# Function to install Bioconductor packages if not already installed
install_bioconductor_packages <- function(pkg_list){
  if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

  for(pkg in pkg_list){
    if(!pkg %in% installed.packages()){
      BiocManager::install(pkg)
      print(paste(pkg, "has been installed successfully!"))
    } else {
      print(paste(pkg, "is already installed."))
    }
  }
}

# Call the functions with your list of packages
install_packages(pkg_list)
install_bioconductor_packages(bioconductor_packages)

library(GEOquery)
library(dplyr)
library(stringr)
library(Biobase)
library(caret)
library(limma)
library(plotly)
library(shiny)
library(bslib)
library(shinydashboard)
library(caret)
library(kernlab)

set.seed(3888)

ui = dashboardPage(
  dashboardHeader(title = "NephroPulse"),

  
  dashboardSidebar(
    
    sidebarMenu(
      menuItem("Dashboard", tabName = "dashboard", icon = icon("home")),
      menuItem("Short-term rejection", tabName = "binary", icon = icon("notes-medical")),
      menuItem("CADI score", tabName = "cadi", icon = icon("notes-medical")),
      menuItem("i-IFTA score", tabName = "ifta", icon = icon("notes-medical")),
      menuItem("Our Prediction Model", tabName = "info", icon = icon("info-circle"))
    )
  ),
  
  dashboardBody(
    tags$head(
      tags$link(rel = "stylesheet", type = "text/css", href = "customtitle.css")
    ),
    
    tags$head(tags$style(HTML("
                              .bold-text {
                                font-weight: bold;
                              }"
                              ))),
    
    tabItems(
      tabItem(tabName = "dashboard",
              h1("Dashboard"),
              br(),
              p("Welcome to the Kidney Transplantation Risk Calculator!"),
              p("This app helps you predict the risk of kidney transplantation for individual patients. 
                To use the app, you will need to upload a CSV file containing patient data."),
              p("Once the data has been uploaded, the app will automaticaly predict the risk for that patient. Please note that the app can predict risk for only one patient at a time. Details of the structure of the .csv file required are provided in each tab."),
              p("For more information about the model, please visit the 'Our Prediction Model' tab."),
              br(),
              h3("Quick Summary of Features"),
              fluidRow(column(4,
                              h4("Short-term rejection"),
                              p("Our model classifies the kidney transplant as experiencing acute rejection or short-term stability based on 
                                patient gene expression profiles.")),
                       column(4,
                              h4("CADI score"),
                              p("The Chronic Allograft Damage Index (CADI) provides a measure of chronic damage that the kidney is experiencing and 
                                can help assess long-term graft survival.")),
                       column(4,
                              h4("i-IFTA score"),
                              p("The Interstitial Fibrosis and Tubular Atrophy (i-IFTA) score is a measure of inflammation or scarring in the kidney."))),
              br(),
              p("Each of our models requires a different set of genes to perform the prediction, the details of which are provided within .csv files wihtin each tab.
                Please ensure that your .csv patient file contains one row representing the patient, with the columns representing the gene expression values.
                You can find a sample patient file below."),
              downloadLink("sample_patient", "Sample Patient File")),
      tabItem(tabName = "binary",
              h1("Short-term rejection"),
              fluidRow(box(status = "primary", fileInput("binary_file", "Choose .csv file", accept = ".csv"))),
              p("Your gene expression profile is most consistent with:"),
              fluidRow(box(tags$div(style = "text-align: center;
                                           font-size: 18px;
                                           font-weight:bold", textOutput("binary_outcome")))),
              br(),
              p(class = "bold-text", "What does your prediction mean?"),
              fluidRow(column(width = 5, 
                           HTML("<p><b><i>Acute rejection</b></i> is a condition where the transplanted organ is
                           recognized as a 'foreign' entity and gets rejected by the immune system. This can
                                occur anytime between one week after the transplantation to 3 months after the
                                transplantation, and is a normal immune response.")),
                       column(width = 5,
                              HTML("<p>If the condition of the transplant is <b><i>stable</i></b>, it is not being
                                   actively rejected by the immune system."))),
              br(),
              p("You can find information on the genes required for this prediction below."),
              downloadLink("binary_download", "Download CSV"),
              br(),
              p("You can also find a template file below."),
              downloadLink("column_names_binary", "Download CSV")),
      tabItem(tabName = "cadi",
              h1("CADI score"),
              fluidRow(box(status = "primary", fileInput("cadi_file", "Choose .csv file", accept = ".csv"))),
              p("You predicted Chronic Allograft Damage Index is:"),
              fluidRow(box(tags$div(style = "text-align: center;
                                    font-size: 18px;
                                    font-weight: bold", textOutput("cadi_outcome")))),
              br(),
              HTML("<p>The <b>Chronic Allograft Damage Index</b>, or <b>CADI score</b>, provides a measure of 
                   chronic damage to the kidney and is closely correlated with later graft function and outcome.
                   The score ranges from a minimum of 0 to a maximum of 18, with a higher score corresponding to
                   a less favourable outcome for allograft status. It is calculated as a sum of
                   six histopathological parameters:"),
              HTML("<ul>
                      <li>interstitial inflammation</li>
                      <li>tubular atrophy</li>
                      <li>interstitial fibrosis</li>
                      <li>arterial fibrointimal thickening</li>
                      <li>glomerular mesangial matrix increase</li>
                      <li>percentage of globally sclerosed glomeruli</li>
                   </ul>"),
              br(),
              p("You can find information on the genes required for this prediction below."),
              downloadLink("cadi_download", "Download CSV"),
              br(),
              p("You can also find a template file below."),
              downloadLink("column_names_cadi", "Download CSV"),
              HTML("<p><b><i>Please note, this prediction is supplementary to our main short-term rejection prediction, and must be used in
                   conjunction with it.</i></b></p>")),
      tabItem(tabName = "ifta",
              h1("i-IFTA score"),
              fluidRow(box(status = "primary", fileInput("ifta_file", "Choose .csv file", accept = ".csv"))),
              fluidRow(box(tags$div(style = "text-align: center;
                                    font-size: 18px;
                                    font-weight: bold", textOutput("ifta_outcome")))),
              br(),
              HTML("<p>The <b>I-IFTA</b> score is a measure of inflammation in the kidney and has been associated
                   with decreased allograft survival. The score ranges from a minimum of 0 to a maximum of 3."),
              HTML("<ul>
                      <li><b>0: </b>no inflammation</li>
                      <li><b>1: </b>mild inflammation</li>
                      <li><b>2: </b>moderate inflammation</li>
                      <li><b>3: </b>severe inflammation</li
                   </ul>"),
              br(),
              p("You can find information on the genes required for this prediction below."),
              downloadLink("ifta_download", "Download CSV"),
              br(),
              p("You can also find a template file below."),
              downloadLink("column_names_ifta", "Download CSV"),
              HTML("<p><b><i>Please note, this prediction is supplementary to our main short-term rejection prediction, and must be used in
                   conjunction with it.</i></b></p>")),
      tabItem(tabName = "info",
  h1("Our Prediction Model"),
  tabsetPanel(
    tabPanel("Short-term Rejection Model", 
      h2("Short-term rejection model"),
      tabsetPanel(
        tabPanel("Selection of k", 
                 plotlyOutput("k_selection_cv"),
                 plotlyOutput("k_selection_ba"),
          HTML("<p>The model employs k-Nearest Neighbors (kNN) classification, al algorithm that works by classifying a data point based on the majority class of its 'k' nearest neighbors. The optimal 'k' of 1 for our case was determined based on specific characteristics of the data and the model's performance.")),
        tabPanel("Balanced Accuracy & cross validation", 
                 plotlyOutput("cv_binary_boxplot"),
          HTML("<p>To validate the model's performance and generalize it to unseen data, it uses cross-validation (CV) with 5 folds and 50 repeats. This resampling procedure involves dividing the data into 5 subsets (folds). For each 
repetition, the model is trained on 4 subsets and tested on the remaining one. This process is repeated 50 
times, with each fold used as the testing set 10 times over the 50 repetitions.")),
        tabPanel("Sensitivity",  
                 plotlyOutput("sensitivity_binary"),
          HTML("<p>The ability of a model to correctly predict that a patient's kidney is experiencing acute rejection, given that the patient is actually experiencing acute rejection. For our model, this measure was 0.96, which means that 96% of the time, the model can correctly identify patients experiencing acute rejection.")),
        tabPanel("Specificity",  
                 plotlyOutput("specificity_binary"),
          HTML("<p>The ability of a model to correctly predict that a patient's kidney is stable, given that the patient's kidney is actually stable. For our model, this measure was 0.97, which means that 97% of the time, the model can correctly identify patients with stable kidneys."))
      )
    ),
    tabPanel("CADI Score Model", 
      h2("CADI score cross validation"),
      tabsetPanel(
        tabPanel("General", 
          HTML("<p>To ensure the model's robustness and assess its performance, cross-validation (CV) with 5 folds and 50 repeats is used. In this process, the dataset is partitioned into 5 subsets. The model is then trained on 4 of 
these subsets, and the remaining subset is used for testing. This is done 50 times, each with a different testing subset, hence 50 repeats.")),
        tabPanel("Mean Absolute Error (MAE)",   
                 plotlyOutput("cadi_mae_plot"),
          HTML("<p>The average difference between the model's predicted CADI score andthe actual score. For our model, this measure was 2.21, which means that the predicted CADI score differs from the actual score by 2.21 points on average. This is the most straightforward measure of our model's accuracy")),
        tabPanel("Root Mean Square Error (RMSE)",    
                 plotlyOutput("cadi_rmse_plot"),
          HTML("<p>A measure of error similar to MAE, but it penalises bigger differences between the model's predicted value and the actual value. For our model, this measure was 2.56.")),
        tabPanel("R Squared",    
                 plotlyOutput("cadi_rsquared_plot"),
          HTML("<p>A measure of how well our model emulates the relationship between gene expression levels and the CADI score. The value ranges from 0 to 1, and the closer the value is to 1, the better the model
can perform. For our model, this value was 0.36."))
      )
    ),
    tabPanel("i-IFTA Score Model",
      h2("i-IFTA score cross validation"),
      tabsetPanel(
        tabPanel("General", 
          HTML("<p>To ensure the model's robustness and assess its performance, cross-validation (CV) with 5 folds and 50 repeats is used. In this process, the dataset is partitioned into 5 subsets. The model is then trained on 4 of 
these subsets, and the remaining subset is used for testing. This is done 50 times, each with a different testing subset, hence 50 repeats.")),
        tabPanel("Mean Absolute Error (MAE)", 
                 plotlyOutput("ifta_mae_plot"),
          HTML("<p>The average difference between the model's predicted i-IFTA score and the actual score. For our model, this measure was 0.68, which means that the predicted CADI score differs from the actual score by 0.68 points on average. This is the most straightforward measure of our model's accuracy.")),
        tabPanel("Root Mean Square Error (RMSE)", 
                 plotlyOutput("ifta_rmse_plot"),
          HTML("<p>A measure of error similar to MAE, but it penalises bigger differences between the model's predicted value and the actual value. For our model, this measure was 0.9.")),
        tabPanel("R Squared", 
                 plotlyOutput("ifta_rsquared_plot"),
          HTML("<p>A measure of how well our model emulates the relationship between gene expression levels and
the i-IFTA score. The value ranges from 0 to 1, and the closer the value is to 1, the better the model can perform. For our model, this value was 0.21."))
      )
    ),
  tabPanel("Resource", 
          h2("Additional Resources"),
              HTML("<p>The datasets used to train this model were all obtained from the <a href='https://www.ncbi.nlm.nih.gov/gds/'>Gene Expression Omnibus database</a>.</p>"),
              HTML("<p>Two datasets were used to train the short-term rejection model: <a href='https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE26578'>GSE26578</a> and <a href='https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14328'>GSE14328</a>.</p>"),
              HTML("<p>For both the CADI and IFTA models, <a href='https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE25902'>GSE25902</a> was used.</p>")
        )
  )
)


    
  )
)
)


server = function(input, output) {
  binary_model = readRDS("model/binary_model.rds")
  ifta_model = readRDS("model/ifta_model.rds")
  cadi_model = readRDS("model/cadi_model.rds")
  
  cv_binary_metrics = read.csv("model/cv_binary_metrics.csv")
  cv_cadi = read.csv("model/cv_cadi.csv")
  cv_ifta = read.csv("model/cv_ifta.csv")
  # When reading from the CSV file
  acc_long <- read.csv("model/acc_long.csv")
  ba_long <- read.csv("model/ba_long.csv")
  
  # Converting 'variable' column into an ordered factor
  acc_long$variable <- factor(acc_long$variable, levels = paste("k =", 1:50))
  ba_long$variable <- factor(ba_long$variable, levels = paste("k =", 1:50))

  
  output$binary_outcome = renderText({
    if (is.null(input$binary_file)) {
      return("Please upload a .csv file to see results.")
    }
    
    file_path = input$binary_file$datapath
    testing_data = read.csv(file_path)
    testing_data = testing_data[, 2:411]
    
    out = class::knn(train = binary_model$train_data, test = testing_data, cl = binary_model$cl, k = binary_model$k)
    
    paste(out)
  })
  
  output$ifta_outcome = renderText({
    if (is.null(input$ifta_file)) {
      return("Please upload a .csv file to see results.")
    }
    
    file_path = input$ifta_file$datapath
    testing_data = read.csv(file_path)
    testing_data = testing_data[, 2:501]
    
    out = round(predict(ifta_model, testing_data)[1])
    
    paste(out)
  })
  
  output$cadi_outcome = renderText({
    if (is.null(input$cadi_file)) {
      return("Please upload a .csv file to see results.")
    }
    
    file_path = input$cadi_file$datapath
    testing_data = read.csv(file_path)
    testing_data = testing_data[, 2:501]
    
    out = round(predict(cadi_model, testing_data)[1])
    
    paste(out)
  })
  
  output$binary_download = downloadHandler(
    filename = function() {
      paste("overview/binary_genes.csv", sep = "")
    },
    content = function(binary_genes) {
      file.copy("overview/binary_genes.csv", binary_genes)
    })
  
  output$sample_patient = downloadHandler(
    filename = function() {
      paste("tests/test_binary.csv", sep = "")
    },
    content = function(sample_patient_file) {
      file.copy("tests/test_binary.csv", sample_patient_file)
    })
  
  output$cadi_download = downloadHandler(
    filename = function() {
      paste("overview/cadi_genes.csv", sep = "")
    },
    content = function(cadi_genes) {
      file.copy("overview/cadi_genes.csv", cadi_genes)
    })
  
  output$ifta_download = downloadHandler(
    filename = function() {
      paste("overview/ifta_genes.csv", sep = "")
    },
    content = function(ifta_genes) {
      file.copy("overview/ifta_genes.csv", ifta_genes)
    })
  
  output$cv_binary_boxplot <- renderPlotly({
  tryCatch({
    cv_binary_metrics_df <- as.data.frame(cv_binary_metrics)  # convert matrix to data frame
    cv_binary_metrics_df$id <- rownames(cv_binary_metrics_df)  # add an id column

    long = cv_binary_metrics_df %>%
    select(c(id, cv_binary, balanced_accs)) %>%
    reshape2::melt(id.vars = "id")

    p = ggplot(long, aes(y = value)) + geom_boxplot() +
    facet_grid(~ variable, labeller = labeller(variable = c(cv_binary = "Cross Validation Accuracy", balanced_accs = "Balanced Accuracy"))) + xlab("") + 
    ylab("Accuracy") + ggtitle("Short-term rejection prediction") + theme(plot.title = element_text(hjust = 0.5)) +
      scale_x_discrete(labels = c()) + theme_bw() + theme(text = element_text(size = 17))
    plotly::ggplotly(p)
  }, error = function(e) {
    print(as.character(e$message))
  })
})
  
  output$k_selection_cv <- renderPlotly({
    p_acc <- ggplot(acc_long, aes(x = variable, y = value)) + 
  geom_boxplot() +
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x = "k", y = "Accuracy", title = "Accuracy for different k")
    plotly::ggplotly(p_acc)
  })
  
  output$k_selection_ba <- renderPlotly({
  p_ba <- ggplot(ba_long, aes(x = variable, y = value)) + 
  geom_boxplot() +
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x = "k", y = "Balanced Accuracy", title = "Balanced Accuracy for different k")
    plotly::ggplotly(p_ba)
  })
  
  output$sensitivity_binary <- renderPlotly({
  tryCatch({
    cv_binary_metrics_df <- as.data.frame(cv_binary_metrics)  # convert matrix to data frame
    cv_binary_metrics_df$id <- rownames(cv_binary_metrics_df)  # add an id column

    long = cv_binary_metrics_df %>%
    select(c(id, se)) %>%
    reshape2::melt(id.vars = "id")

    p = ggplot(long, aes(y = value)) + geom_boxplot() +
    facet_grid(~ variable, labeller = labeller(variable = c(se = "Sensitivity"))) + xlab("") + 
    ylab("Accuracy") + ggtitle("Short-term rejection Sensitivity") + theme(plot.title = element_text(hjust = 0.5)) +
      scale_x_discrete(labels = c()) + theme_bw() + theme(text = element_text(size = 17))
    plotly::ggplotly(p)
  }, error = function(e) {
    print(as.character(e$message))
  })
})

output$specificity_binary <- renderPlotly({
  tryCatch({
    cv_binary_metrics_df <- as.data.frame(cv_binary_metrics)  # convert matrix to data frame
    cv_binary_metrics_df$id <- rownames(cv_binary_metrics_df)  # add an id column

    long = cv_binary_metrics_df %>%
    select(c(id, sp)) %>%
    reshape2::melt(id.vars = "id")

    p = ggplot(long, aes(y = value)) + geom_boxplot() +
    facet_grid(~ variable, labeller = labeller(variable = c(sp = "Specificity"))) + xlab("") + 
    ylab("Accuracy") + ggtitle("Short-term rejection Specificity") + theme(plot.title = element_text(hjust = 0.5)) +
      scale_x_discrete(labels = c()) + theme_bw() + theme(text = element_text(size = 17))
    plotly::ggplotly(p)
  }, error = function(e) {
    print(as.character(e$message))
  })
})

  
 output$cadi_mae_plot = renderPlotly({
  long = cv_cadi %>% select(X, mae) %>% as.data.frame() %>% reshape2::melt(id.vars = "X")
  p = ggplot(long, aes(y = value)) + geom_boxplot() +
    ggtitle("CADI Mean Absolute Error") + 
    theme(plot.title = element_text(hjust = 0.5)) + theme_bw() + theme(text = element_text(size = 17))
  plotly::ggplotly(p)
})

output$cadi_rmse_plot = renderPlotly({
  long = cv_cadi %>% select(X, rmse) %>% as.data.frame() %>% reshape2::melt(id.vars = "X")
  p = ggplot(long, aes(y = value)) + geom_boxplot() +
    ggtitle("CADI Root Mean Squared Error") + 
    theme(plot.title = element_text(hjust = 0.5)) + theme_bw() + theme(text = element_text(size = 17))
  plotly::ggplotly(p)
})

output$cadi_rsquared_plot = renderPlotly({
  long = cv_cadi %>% select(X, rsquared) %>% as.data.frame() %>% reshape2::melt(id.vars = "X")
  p = ggplot(long, aes(y = value)) + geom_boxplot() +
    ggtitle("CADI R Squared over") + 
    theme(plot.title = element_text(hjust = 0.5)) + theme_bw() + theme(text = element_text(size = 17))
  plotly::ggplotly(p)
})

output$ifta_mae_plot = renderPlotly({
  long = cv_ifta %>% select(X, mae) %>% as.data.frame() %>% reshape2::melt(id.vars = "X")
  p = ggplot(long, aes(y = value)) + geom_boxplot() +
    ggtitle("i-IFTA Mean Absolute Error") + 
    theme(plot.title = element_text(hjust = 0.5)) + theme_bw() + theme(text = element_text(size = 17))
  plotly::ggplotly(p)
})

output$ifta_rmse_plot = renderPlotly({
  long = cv_ifta %>% select(X, rmse) %>% as.data.frame() %>% reshape2::melt(id.vars = "X")
  p = ggplot(long, aes(y = value)) + geom_boxplot() +
    ggtitle("i-IFTA Root Mean Squared Error") + 
    theme(plot.title = element_text(hjust = 0.5)) + theme_bw() + theme(text = element_text(size = 17))
  plotly::ggplotly(p)
})

output$ifta_rsquared_plot = renderPlotly({
  long = cv_ifta %>% select(X, rsquared) %>% as.data.frame() %>% reshape2::melt(id.vars = "X")
  p = ggplot(long, aes(y = value)) + geom_boxplot() +
    ggtitle("i-IFTA R Squared") + 
    theme(plot.title = element_text(hjust = 0.5)) + theme_bw() + theme(text = element_text(size = 17))
  plotly::ggplotly(p)
})

  
  output$column_names_binary = downloadHandler(
    filename = function() {
      paste("overview/column_names.csv", sep = "")
    },
    content = function(column_names_binary) {
      file.copy("overview/column_names_binary", column_names_binary)
    })
  
  output$column_names_cadi = downloadHandler(
    filename = function() {
      paste("overview/column_names_cadi.csv", sep = "")
    },
    content = function(column_names_cadi) {
      file.copy("overview/column_names_cadi.csv", column_names_cadi)
    })
  
  output$column_names_ifta = downloadHandler(
    filename = function() {
      paste("overview/column_names_ifta.csv", sep = "")
    },
    content = function(column_names_ifta) {
      file.copy("overview/column_names_ifta.csv", column_names_ifta)
    })
}
shinyApp(ui = ui, server = server)
```