app.R

# ==============================================
# Shiny App: ExomDash
# Final Version 4.0 - Portfolio Ready
# ==============================================

# --- 1. LOAD PACKAGES ---
library(shiny)
library(shinydashboard)
library(data.table)
library(DT)
library(ggplot2)
library(maftools)
library(ggrepel)
library(dplyr)
library(future)
library(promises)
library(plotly)
library(survival)
library(survminer)
library(pheatmap)
library(Rtsne) 
library(umap)
library(ggpubr)
library(tidyr)
library(grid)

# --- 2. SETTINGS AND HELPER FUNCTIONS ---
plan(multisession)
options(shiny.maxRequestSize = 200*1024^2)

exportPlot <- function(id, plotFun) {
  downloadHandler(
    filename = function() paste0(id, ".png"),
    content = function(file) {
      p_obj <- plotFun()
      if (inherits(p_obj, "pheatmap")) {
        png(file, width = 12, height = 10, units = "in", res = 300)
        grid::grid.draw(p_obj$gtable)
        dev.off()
      } else {
        ggsave(file, plot = p_obj, width = 10, height = 8, dpi = 300)
      }
    }
  )
}

exportTable <- function(id, dataFun) {
  downloadHandler(
    filename = function() paste0(id, ".csv"),
    content = function(file) {
      write.csv(dataFun(), file, row.names = FALSE)
    }
  )
}

# --- ACCESS PASSWORD ---
# Change "YOUR_PASSWORD" to your desired password
passwd <- "YOUR_PASSWORD"

# ==============================================
# 3. USER INTERFACE (UI)
# ==============================================

# --- UI Part 1: LOGIN SCREEN ---
login_ui <- fluidPage(
  titlePanel("Restricted Access to ExomDash"),
  wellPanel(
    fluidRow(
      column(4,
             passwordInput("password", "Please enter the password:"),
             actionButton("login_button", "Login", class = "btn-primary"),
             br(), br(),
             p("This is a protected dashboard. Access is restricted to authorized collaborators.")
      )
    )
  )
)

# --- UI Part 2: MAIN DASHBOARD ---
main_ui <- navbarPage(
  "ExomDash",
  
  tabPanel("Load Data",
           sidebarLayout(
             sidebarPanel(
               h4("1. Upload Required Files"),
               fileInput("maf_file", "MAF File (.txt or .maf)", accept = c(".txt", ".maf")),
               fileInput("clinical_file", "Clinical Data (.csv or .tsv)", accept = c(".csv", ".tsv")),
               tags$hr(),
               h4("2. Upload Optional Files"),
               fileInput("cnv_file", "CNV Data (for Oncoplot, .csv)", accept = c(".csv")),
               fileInput("cov_file", "Coverage Data (mosdepth, .csv)", accept = c(".csv"))
             ),
             mainPanel(
               div(style="text-align: center;",
                   tags$img(src = "logo.jpg", width = "400px"),
                   hr()
               ),
               h3("ExomDash: From Variant to Clinical Insight"),
               tags$blockquote(
                 p("The exome represents the essential chapters of our genome, containing the recipes for all the proteins that govern our body. Small variations in these recipes, sometimes the change of a single letter (SNP), make us unique, but can also be the key to understanding various diseases. The scientist's task is to act as a genomic detective, and ExomDash is the tool developed for this mission: an interactive platform to navigate the complexity of the exome, identify the variants that matter, and connect these genomic discoveries to the clinical story of each patient. Born from the need to transform a sea of data into clear insights, ExomDash aims to accelerate the discovery and understanding of the genomic landscape.")
               )
             )
           )
  ),
  
  tabPanel("Overview",
           fluidRow(
             valueBoxOutput("nsamples_box"),
             valueBoxOutput("ngenes_box"),
             valueBoxOutput("nmutations_box")
           ),
           fluidRow(
             box(title = "Mutation Summary", width = 12, solidHeader = TRUE, status = "primary",
                 plotOutput("maf_summary_plot", height = "500px"),
                 downloadButton("download_maf_summary_plot", "Download PNG")
             )
           )
  ),
  
  tabPanel("Clinical Analysis",
           sidebarLayout(
             sidebarPanel(
               h4("Exploratory Analysis"),
               selectInput("clinical_var_selector", "Select a clinical variable:", choices = NULL)
             ),
             mainPanel(fluidRow(
               box(title = "Frequency", width = 6, plotOutput("clinical_bar_plot"), downloadButton("download_clinical_bar", "Download PNG")),
               box(title = "Proportion", width = 6, plotOutput("clinical_donut_plot"), downloadButton("download_clinical_donut", "Download PNG"))
             ))
           )
  ),
  
  navbarMenu("Mutation Landscape",
             tabPanel("Oncoplot",
                      sidebarLayout(
                        sidebarPanel(
                          h4("Plot Options"),
                          numericInput("oncoplot_top_n", "Top N Genes:", 20, min = 5, max = 50),
                          checkboxInput("show_barcode", "Show sample names", value = TRUE),
                          tags$hr(),
                          h4("Clinical Annotations"),
                          selectInput("clinical_features", "Select features:", choices = NULL, multiple = TRUE)
                        ),
                        mainPanel(
                          plotOutput("oncoplot", height = "800px"),
                          downloadButton("download_oncoplot", "Download PNG")
                        )
                      )
             ),
             tabPanel("TMB (Tumor Mutational Burden)",
                      sidebarLayout(
                        sidebarPanel(
                          h4("TMB Analysis Options"),
                          p("Tumor Mutational Burden (TMB) is calculated as the number of mutations per Megabase (Mb)."),
                          tags$hr(),
                          numericInput("exome_size", "Exome/Panel Size (Mb):", value = 38, min = 1, max = 3000),
                          checkboxInput("tmb_log_scale", "Use logarithmic scale (Y-axis)", value = TRUE),
                          tags$hr(),
                          h4("Clinical Comparisons"),
                          selectInput("tmb_clinical_var", "Compare TMB by variable:", choices = NULL),
                          selectInput("tmb_stat_test", "Statistical test:", 
                                      choices = c("None" = "Nenhum", 
                                                  "Overall (Kruskal-Wallis)" = "kruskal.test", 
                                                  "Pairwise (Wilcoxon)" = "wilcox.test"))
                        ),
                        mainPanel(
                          fluidRow(
                            column(8, 
                                   h4("TMB Distribution by Clinical Group"),
                                   plotOutput("tmb_boxplot")
                            ),
                            column(4, 
                                   h4("TMB vs. Mutation Count"),
                                   plotOutput("tmb_correlation_plot")
                            )
                          ),
                          fluidRow(
                            column(6, downloadButton("download_tmb_boxplot", "Download Boxplot")),
                            column(6, downloadButton("download_tmb_correlation", "Download Correlation Plot"))
                          ),
                          tags$hr(),
                          h4("TMB Data Table"),
                          DT::DTOutput("tmb_table"),
                          downloadButton("download_tmb_table", "Download Table (CSV)")
                        )
                      )
            ),
             tabPanel("Ti/Tv Summary",
                      h3("Transitions and Transversions (Ti/Tv) Summary"),
                      plotOutput("titv_plot", height = "600px"),
                      downloadButton("download_titv_plot", "Download PNG")
             ),
             tabPanel("Rainfall Plot",
                      sidebarLayout(
                        sidebarPanel(
                          h4("Sample Selection"),
                          selectInput("rainfall_sample_selector", "Select a sample:", choices = NULL)
                        ),
                        mainPanel(
                          h3("Rainfall Plot by Sample"),
                          plotOutput("rainfall_plot", height = "600px"),
                          downloadButton("download_rainfall_plot", "Download PNG")
                        )
                      )
             )
  ),
  
  navbarMenu("Interactions & Pathways",
             tabPanel("Somatic Interactions",
                      sidebarLayout(
                        sidebarPanel(
                          numericInput("somatic_top_n", "Analyze Top N Genes:", 25, min=10),
                          p("Shows co-occurring (green) or mutually exclusive (gold) gene pairs.")
                        ),
                        mainPanel(
                          h3("Somatic Interactions Among Top Mutated Genes"),
                          plotOutput("somatic_interactions_plot", height = "700px"),
                          downloadButton("download_somatic_plot", "Download PNG")
                        )
                      )
             ),
             tabPanel("Metabolic Pathways (Oncoplot)",
                      sidebarLayout(
                        sidebarPanel(
                          numericInput("pathway_top_genes", "Top N Genes:", 10, min=5, max=50)
                        ),
                        mainPanel(
                          h3("Oncoplot Grouped by Curated Signaling Pathways"),
                          plotOutput("pathway_plot", height="700px"),
                          downloadButton("download_pathway_plot", "Download PNG")
                        )
                      )
             )
  ),

  navbarMenu("Advanced Analyses",
    tabPanel("Survival Analysis",
             sidebarLayout(
               sidebarPanel(
                 selectInput("surv_gene", "Select Gene:", choices = NULL),
                 p("This analysis requires 'time' and 'Status' columns in the clinical data.")
               ),
               mainPanel(
                 plotOutput("survPlot"),
                 downloadButton("downloadSurvPlot", "Download PNG")
               )
             )
    ),
    tabPanel("Population Frequency (gnomAD)",
             sidebarLayout(
               sidebarPanel(
                 h4("Population Allele Frequency Analysis"),
                 p("This analysis explores the frequency of variants (from gnomAD) found in your cohort, grouped by a clinical variable."),
                 tags$hr(),
                 selectizeInput("gnomad_gene_selector", "Select Gene(s):", choices = NULL, multiple = TRUE),
                 selectInput("gnomad_clinical_selector", "Group by clinical variable:", choices = NULL),
                 checkboxInput("gnomad_log_scale", "Use logarithmic scale (Y-axis)", value = FALSE),
                 tags$hr(),
                 selectInput("gnomad_stat_test", "Statistical test:", 
                             choices = c("None" = "Nenhum", 
                                         "Overall (Kruskal-Wallis)" = "kruskal.test", 
                                         "Pairwise (Wilcoxon)" = "wilcox.test")),
                 tags$hr(),
                 p("The table below shows the filtered data used to generate the plot.")
               ),
               mainPanel(
                 plotOutput("gnomad_plot", height = "600px"),
                 downloadButton("download_gnomad_plot", "Download Plot (PNG)"),
                 tags$hr(),
                 DT::DTOutput("gnomad_table"),
                 downloadButton("download_gnomad_table", "Download Table (CSV)")
               )
             )
    ),
    tabPanel("Pathway Analysis",
             sidebarLayout(
               sidebarPanel(
                 h4("Oncoplot by Signaling Pathways"),
                 p("Visualize mutations grouped by oncogenic pathways and correlate them with clinical data."),
                 tags$hr(),
                 numericInput("pathway_oncoplot_top_n", "Top N Genes to display:", 20, min = 10, max = 100),
                 checkboxInput("pathway_show_samples", "Show sample names", value = FALSE),
                 tags$hr(),
                 selectizeInput("pathway_clinical_features", "Select clinical features for annotation:", 
                                choices = NULL, multiple = TRUE)
               ),
               mainPanel(
                 plotOutput("pathway_oncoplot_clinical", height = "800px"),
                 downloadButton("download_pathway_oncoplot_clinical", "Download Plot (PNG)")
               )
             )
    ),
    tabPanel("Exploratory Analysis",
             sidebarLayout(
                sidebarPanel(
                    h4("1. Analysis Options"),
                    selectInput("multiMethod", "Reduction Method:", choices = c("PCA", "t-SNE", "UMAP")),
                    tags$hr(),
                    h4("2. Grouping Options"),
                    selectInput("clinVar", "Group by clinical variable:", choices = NULL),
                    uiOutput("dynamic_filter_ui")
                ),
                mainPanel(
                    plotOutput("multiPlot", height = "600px"),
                    downloadButton("downloadMultiPlot", "Download PNG")
                )
             )
    )
  ),
  
  tabPanel("Coverage Quality",
     mainPanel(width = 12,
        h4("Mean Coverage per Sample"),
        plotOutput("covBoxPlot"),
        downloadButton("downloadCovBox", "Download PNG")
     )
  ),
  
  tabPanel("About", icon = icon("info-circle"),
           fluidPage(
             titlePanel("About ExomDash"),
             hr(),
             fluidRow(
               column(8,
                 h3("Abstract"),
                 p("This dashboard was created for the interactive analysis and visualization of somatic mutation data (MAF format), typically generated from whole-exome sequencing (WES) or gene panels. It integrates mutation data, copy number variations (CNV), and clinical information to generate visualizations and analyses that allow for the exploration of a sample cohort's genomic landscape."),
                 hr(),
                 h3("Input File Formats"),
                 tags$b("1. MAF File:"),
                 p("A text file (TSV) containing variant annotations. It must follow the MAF format and contain at least the following columns: ", tags$code("Hugo_Symbol"), ",", tags$code("Tumor_Sample_Barcode"), ",", tags$code("Variant_Classification"), ",", tags$code("Chromosome"), ",", tags$code("Start_Position"), "."),
                 tags$b("2. Clinical Data:"),
                 p("A text file (CSV or TSV) where the first column must be ", tags$code("Tumor_Sample_Barcode"), ", corresponding to the samples in the MAF file. Subsequent columns can contain any clinical variable (e.g., 'age', 'molecular_subtype', 'vital_status')."),
                 tags$b("3. CNV Data (Optional):"),
                 p("A CSV file with 3 required columns:", tags$code("Hugo_Symbol"), "(gene),", tags$code("Tumor_Sample_Barcode"), "(sample), and", tags$code("cn"), "(alteration, with 'Amp' for amplification or 'Del' for deletion)."),
                 hr(),
                 h3("How to Cite"),
                 p("If you use this dashboard in your work, please cite it as follows:"),
                 tags$blockquote("Falco, M. L. (2025). ExomDash: An Interactive Dashboard for Exomic Data Analysis. (Version 4.0). Available at: https://github.com/mlfalco-bioinfo/exomdash"),
                 hr(),
                 h3("License"),
                 p("This project is distributed under the [License Name, e.g., MIT] License. See the LICENSE file in the repository for more details."),
                 hr(),
                 h3("Acknowledgments"),
                 p("Acknowledgments to CNPq and Fundação Araucária.")
               ),
               column(4, 
                  tags$img(src = "logo.jpg", width = "200px", style="display: block; margin-left: auto; margin-right: auto;"),
                  hr(),
                  h4("Developer"),
                  p("Mateus L. Falco"),
                  p(tags$b("Email:"), tags$a(href="mailto:mlfalco.bioinfo@gmail.com", "mlfalco.bioinfo@gmail.com")),
                  p(tags$b("GitHub:"), tags$a(href="https://github.com/mlfalco-bioinfo", target="_blank", "mlfalco-bioinfo")),
                  hr(),
                  h4("Reproducibility"),
                  p("For details about the R package versions used in this session, click the button below."),
                  actionButton("show_session_info", "Show Session Info", icon = icon("laptop-code"), style="width: 100%;")
               )
             )
           )
  )
)

# --- Final UI: Controller ---
ui <- fluidPage(
  uiOutput("app_ui")
)

# ==============================================
# 4. SERVER LOGIC
# ==============================================
server <- function(input, output, session) {

  # --- AUTHENTICATION LOGIC ---
  user_authenticated <- reactiveVal(FALSE)

  output$app_ui <- renderUI({
    if (user_authenticated()) {
      return(main_ui)
    } else {
      return(login_ui)
    }
  })

  observeEvent(input$login_button, {
    if (!is.null(input$password) && input$password == passwd) {
      user_authenticated(TRUE)
    } else {
      showModal(modalDialog(
        title = "Access Error",
        "Incorrect password. Please try again.",
        easyClose = TRUE,
        footer = modalButton("Close")
      ))
    }
  })
  
  # --- MAIN DASHBOARD LOGIC (RUNS ONLY AFTER SUCCESSFUL LOGIN) ---
  
  maf_obj_base <- reactiveVal(NULL)
  cnv_data_raw <- reactiveVal(NULL)

  observeEvent(list(input$maf_file, input$clinical_file), {
    req(user_authenticated(), input$maf_file)
    showNotification("Reading and processing MAF and Clinical files...", type = "message", duration = NULL, id = "loading_maf")
    maf_path <- input$maf_file$datapath
    clin_path <- if (!is.null(input$clinical_file)) input$clinical_file$datapath else NULL
    future({
      clin_df <- if (!is.null(clin_path)) fread(clin_path) else NULL
      read.maf(maf = maf_path, clinicalData = clin_df, verbose = FALSE)
    }) %...>%
      maf_obj_base()
  })

  observeEvent(input$cnv_file, {
    req(user_authenticated())
    showNotification("Reading CNV file...", type = "message", duration = 4)
    df <- fread(input$cnv_file$datapath)
    cnv_data_raw(df)
  })
  
  maf_object <- reactive({
    req(user_authenticated())
    maf_obj <- maf_obj_base()
    req(maf_obj)
    cnv_df <- cnv_data_raw()
    if (!is.null(cnv_df)) {
      cnv_df_formatted <- cnv_df %>%
        select(any_of(c("Hugo_Symbol", "Tumor_Sample_Barcode", "cn"))) %>%
        rename(Variant_Classification = cn)
      maf_obj@cnv.data <- cnv_df_formatted
    }
    removeNotification("loading_maf")
    showNotification("Data ready for analysis!", type = "message", duration = 5)
    return(maf_obj)
  })

  coverage_data <- eventReactive(input$cov_file, { 
    req(user_authenticated(), input$cov_file)
    fread(input$cov_file$datapath) 
  })
  
  observe({
    req(user_authenticated(), maf_object())
    
    genes <- getGeneSummary(maf_object())$Hugo_Symbol
    samples <- as.character(getSampleSummary(maf_object())$Tumor_Sample_Barcode)
    
    updateSelectInput(session, "rainfall_sample_selector", choices = samples, selected = samples[1])
    updateSelectInput(session, "surv_gene", choices = genes, selected = genes[1])
    updateSelectizeInput(session, "gnomad_gene_selector", choices = genes, server = TRUE)
    
    if (!is.null(maf_object()@clinical.data)){
        clin_choices <- colnames(maf_object()@clinical.data)
        clin_choices <- clin_choices[clin_choices != "Tumor_Sample_Barcode"]
        
        if (length(clin_choices) > 0) {
          updateSelectInput(session, "clinical_var_selector", choices = clin_choices, selected = clin_choices[1])
          updateSelectInput(session, "clinical_features", choices = clin_choices)
          updateSelectizeInput(session, "pathway_clinical_features", choices = clin_choices)
          updateSelectInput(session, "clinVar", choices = clin_choices, selected = clin_choices[1])
          updateSelectInput(session, "gnomad_clinical_selector", choices = clin_choices, selected = clin_choices[1])
          updateSelectInput(session, "tmb_clinical_var", choices = clin_choices, selected = clin_choices[1])
        }
    }
  })

  observeEvent(input$show_session_info, {
    req(user_authenticated())
    showModal(modalDialog(
      title = "R Session Info",
      verbatimTextOutput("session_info_output"),
      easyClose = TRUE,
      footer = modalButton("Close")
    ))
  })
  
  output$session_info_output <- renderPrint({
    sessionInfo()
  })

  # --- TAB LOGIC ---
  
  # -- Overview --
  output$maf_summary_text <- renderPrint({ req(maf_object()); maf_object() })
  output$nsamples_box <- renderValueBox({ req(maf_object()); valueBox(getSampleSummary(maf_object())[, .N], "Samples", icon = icon("vials"), color = "purple") })
  output$ngenes_box <- renderValueBox({ req(maf_object()); valueBox(getGeneSummary(maf_object())[, .N], "Mutated Genes", icon = icon("dna"), color = "blue") })
  output$nmutations_box <- renderValueBox({ req(maf_object()); valueBox(maf_object()@data[,.N], "Total Mutations", icon=icon("bolt"), color="yellow") })
  
  maf_summary_plot_reactive <- reactive({ 
    req(maf_object())
    plotmafSummary(maf = maf_object(), rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE, fs = 1.2) 
  })
  output$maf_summary_plot <- renderPlot({ maf_summary_plot_reactive() })
  output$download_maf_summary_plot <- exportPlot("MAF_Summary", maf_summary_plot_reactive)
  
  # -- Clinical Analysis --
  clinical_plot_data <- reactive({ 
    req(maf_object(), input$clinical_var_selector)
    as.data.frame(maf_object()@clinical.data) %>% 
      count(!!sym(input$clinical_var_selector), name = "count") %>% 
      rename(category = 1) %>% 
      filter(!is.na(category), category != "") 
  })
  clinical_bar_plot_obj <- reactive({ 
    plot_data <- clinical_plot_data()
    req(nrow(plot_data) > 0)
    ggplot(plot_data, aes(y = count, x = reorder(as.character(category), count))) + 
      geom_bar(stat = "identity", fill = "steelblue", alpha=0.8) + 
      coord_flip() +
      geom_text(aes(label = count), hjust = -0.3, fontface = "bold", size = 6, color="black") +
      scale_y_continuous(expand = expansion(mult = c(0, 0.18))) +
      labs(y = "Count", x = input$clinical_var_selector) + 
      theme_minimal(base_size = 18)
  })
  clinical_donut_plot_obj <- reactive({
    df <- clinical_plot_data()
    req(nrow(df) > 0)
    df <- df %>% mutate(fraction = count / sum(count), ymax = cumsum(fraction), ymin = c(0, head(ymax, n = -1)), 
                        labelPosition = (ymax + ymin) / 2, 
                        label = paste0(category, "\n (", scales::percent(fraction, accuracy = 0.1), ")"))
    ggplot(df, aes(ymax = ymax, ymin = ymin, xmax = 4, xmin = 3, fill = as.character(category))) + 
      geom_rect() + 
      geom_text(aes(x = 3.5, y = labelPosition, label = label), size = 6, fontface="bold", color="white") +
      coord_polar(theta = "y") + 
      xlim(c(2, 4.5)) + 
      theme_void() + 
      theme(legend.position = "none")
  })
  output$clinical_bar_plot <- renderPlot({ clinical_bar_plot_obj() })
  output$clinical_donut_plot <- renderPlot({ clinical_donut_plot_obj() })
  output$download_clinical_bar <- exportPlot("Clinical_Bar_Plot", clinical_bar_plot_obj)
  output$download_clinical_donut <- exportPlot("Clinical_Donut_Plot", clinical_donut_plot_obj)
  
  # -- Mutation Landscape --
  oncoplot_reactive <- reactive({
      req(maf_object(), input$oncoplot_top_n)
      clin_features <- if (is.null(input$clinical_features) || input$clinical_features == "") NULL else input$clinical_features
      oncoplot(maf = maf_object(), top = input$oncoplot_top_n, 
               showTumorSampleBarcodes = input$show_barcode, clinicalFeatures = clin_features,
               fontSize = 1.2, gene_mar = 8, legendFontSize = 1.4, annotationFontSize = 1.4)
  })
  output$oncoplot <- renderPlot({ oncoplot_reactive() })
  output$download_oncoplot <- exportPlot("Oncoplot", oncoplot_reactive)
  
  tmb_results <- reactive({
    req(maf_object(), input$exome_size)
    tmb(maf = maf_object(), captureSize = input$exome_size, logScale = FALSE)
  })
  tmb_clinical_data <- reactive({
    req(tmb_results(), maf_object()@clinical.data, input$tmb_clinical_var)
    tmb_df <- tmb_results()
    clinical_df <- as.data.frame(maf_object()@clinical.data)
    left_join(tmb_df, clinical_df, by = "Tumor_Sample_Barcode") %>%
      filter(!is.na(.data[[input$tmb_clinical_var]]) & .data[[input$tmb_clinical_var]] != "")
  })
  tmb_boxplot_obj <- reactive({
    plot_data <- tmb_clinical_data()
    req(plot_data, nrow(plot_data) > 0)
    p <- ggplot(plot_data, aes(x = .data[[input$tmb_clinical_var]], y = total_per_mb, fill = .data[[input$tmb_clinical_var]])) +
      geom_boxplot(outlier.shape = NA, alpha = 0.7) +
      geom_jitter(width = 0.1, alpha = 0.7) +
      labs(x = input$tmb_clinical_var, y = "TMB (mutations/Mb)") +
      theme_minimal(base_size = 18) +
      theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1, size=14),
            plot.title = element_text(size=20, face="bold"))
    if (input$tmb_log_scale) {
      p <- p + scale_y_log10() + ylab("TMB (mutations/Mb) - Log10 Scale")
    }
    if (input$tmb_stat_test != "Nenhum") {
      n_groups <- length(unique(plot_data[[input$tmb_clinical_var]]))
      if (input$tmb_stat_test == "wilcox.test" && n_groups > 1) {
        comparisons_list <- combn(unique(as.character(plot_data[[input$tmb_clinical_var]])), 2, simplify = FALSE)
        p <- p + stat_compare_means(method = "wilcox.test", comparisons = comparisons_list, label = "p.signif")
      } else if (input$tmb_stat_test == "kruskal.test" && n_groups > 2) {
        p <- p + stat_compare_means(method = "kruskal.test", label.y.npc = 0.9, label.x.npc = 0.5, size=6)
      }
    }
    p
  })
  tmb_correlation_plot_obj <- reactive({
    ggplot(tmb_results(), aes(x = total, y = total_per_mb)) +
      geom_point(alpha = 0.8, size=3) +
      geom_smooth(method = "lm", se = FALSE) +
      labs(x = "Total Mutation Count", y = "TMB (mutations/Mb)") +
      theme_minimal(base_size = 18)
  })
  output$tmb_boxplot <- renderPlot({ tmb_boxplot_obj() })
  output$tmb_correlation_plot <- renderPlot({ tmb_correlation_plot_obj() })
  output$tmb_table <- DT::renderDT({
    req(tmb_clinical_data())
    tmb_clinical_data() %>%
      select(Tumor_Sample_Barcode, total, total_per_mb, !!sym(input$tmb_clinical_var)) %>%
      rename(Sample = Tumor_Sample_Barcode, N_Mutations = total, TMB = total_per_mb) %>%
      arrange(desc(TMB)) %>%
      DT::datatable(options = list(pageLength = 5, scrollX = TRUE), filter="top", rownames=FALSE)
  })
  output$download_tmb_boxplot <- exportPlot("TMB_Boxplot", tmb_boxplot_obj)
  output$download_tmb_correlation <- exportPlot("TMB_Correlation", tmb_correlation_plot_obj)
  output$download_tmb_table <- exportTable("TMB_Data", tmb_clinical_data)

  titv_plot_reactive <- reactive({ 
    req(maf_object())
    titv_results <- titv(maf = maf_object(), plot = FALSE, useSyn = TRUE)
    plotTiTv(res = titv_results) 
  })
  output$titv_plot <- renderPlot({ titv_plot_reactive() })
  output$download_titv_plot <- exportPlot("TiTv_Plot", titv_plot_reactive)
  
  rainfall_plot_reactive <- reactive({ 
    req(maf_object(), input$rainfall_sample_selector)
    rainfallPlot(maf = maf_object(), tsb = input$rainfall_sample_selector, detectChangePoints = TRUE, pointSize = 0.6) 
  })
  output$rainfall_plot <- renderPlot({ rainfall_plot_reactive() })
  output$download_rainfall_plot <- exportPlot("Rainfall_Plot", rainfall_plot_reactive)
  
  # -- Interactions & Pathways --
  somatic_plot_reactive <- reactive({ 
    req(maf_object(), input$somatic_top_n)
    somaticInteractions(maf = maf_object(), top = input$somatic_top_n, pvalue = c(0.05, 0.1)) 
  })
  output$somatic_interactions_plot <- renderPlot({ somatic_plot_reactive() })
  output$download_somatic_plot <- exportPlot("Somatic_Interactions_Plot", somatic_plot_reactive)

  pathway_plot_reactive <- reactive({ 
    req(maf_object(), input$pathway_top_genes)
    oncoplot(maf = maf_object(), top = input$pathway_top_genes, pathways = "sigpw", fontSize = 1.2, gene_mar = 8) 
  })
  output$pathway_plot <- renderPlot({ pathway_plot_reactive() })
  output$download_pathway_plot <- exportPlot("Metabolic_Pathways_Plot", pathway_plot_reactive)
  
  # -- Advanced Analyses --
  survPlot_reactive <- reactive({
    req(maf_object(), input$surv_gene)
    tryCatch({
      mafSurvival(maf = maf_object(), genes = input$surv_gene, time = 'time', Status = 'Status', isTCGA = TRUE)
    }, error = function(e){
      plot.new(); text(0.5, 0.5, "Error: Please check if clinical data\ncontains 'time' and 'Status' columns.")
    })
  })
  output$survPlot <- renderPlot({ survPlot_reactive() })
  output$downloadSurvPlot <- exportPlot("Survival_Analysis", survPlot_reactive)

  gnomad_data_reactive <- reactive({
    req(maf_object(), input$gnomad_gene_selector, input$gnomad_clinical_selector)
    maf_df <- maf_object()@data
    validate(need("gnomADe_AF" %in% names(maf_df), "Column 'gnomADe_AF' not found in MAF file."))
    maf_df %>%
      mutate(gnomADe_AF = suppressWarnings(as.numeric(as.character(gnomADe_AF)))) %>%
      filter(!is.na(gnomADe_AF)) %>%
      filter(Hugo_Symbol %in% input$gnomad_gene_selector) %>%
      left_join(as.data.frame(maf_object()@clinical.data), by = "Tumor_Sample_Barcode") %>%
      filter(!is.na(.data[[input$gnomad_clinical_selector]]) & .data[[input$gnomad_clinical_selector]] != "")
  })
  gnomad_plot_obj <- reactive({
    plot_data <- gnomad_data_reactive()
    validate(need(is.data.frame(plot_data) && nrow(plot_data) > 0, "No mutations with gnomAD_AF data found for the selected filters."))
    
    p <- ggplot(plot_data, aes(x = .data[[input$gnomad_clinical_selector]], y = gnomADe_AF, fill = .data[[input$gnomad_clinical_selector]])) +
      geom_violin(trim = FALSE, alpha = 0.4) +
      geom_jitter(width = 0.1, alpha = 0.6, height=0) +
      stat_summary(fun = median, geom = "point", shape = 18, size = 5, color = "black") +
      labs(title = paste("Population Frequency (gnomAD) of Mutations in", paste(input$gnomad_gene_selector, collapse=", ")),
           x = input$gnomad_clinical_selector, y = "gnomAD Allele Frequency") +
      theme_minimal(base_size = 18) +
      theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1, size=14),
            plot.title = element_text(size=20, face="bold"))
      
    if (input$gnomad_log_scale) {
      p <- p + scale_y_log10(labels = scales::scientific) + ylab("gnomAD Allele Frequency (Log10 Scale)")
    }
    
    if (input$gnomad_stat_test != "Nenhum") {
      n_groups <- length(unique(plot_data[[input$gnomad_clinical_selector]]))
      if (input$gnomad_stat_test == "wilcox.test" && n_groups > 1) {
        comparisons_list <- combn(unique(as.character(plot_data[[input$gnomad_clinical_selector]])), 2, simplify = FALSE)
        p <- p + stat_compare_means(method = "wilcox.test", comparisons = comparisons_list, label = "p.signif", size=6)
      } else if (input$gnomad_stat_test == "kruskal.test" && n_groups > 2) {
        p <- p + stat_compare_means(method = "kruskal.test", label.y.npc = 0.95, label.x.npc = 0.5, size=7)
      }
    }
    p
  })
  
  gnomad_data_for_table <- reactive({
    df <- gnomad_data_reactive()
    req(is.data.frame(df) && nrow(df) > 0)
    df %>%
      select(Hugo_Symbol, Tumor_Sample_Barcode, Variant_Classification, HGVSp_Short, !!sym(input$gnomad_clinical_selector), gnomADe_AF) %>%
      arrange(gnomADe_AF)
  })
  output$gnomad_table <- DT::renderDT(server = TRUE, {
    DT::datatable(gnomad_data_for_table(),
                  options = list(pageLength = 10, scrollX = TRUE), 
                  filter = "top", 
                  rownames = FALSE)
  })
  
  output$gnomad_plot <- renderPlot({ gnomad_plot_obj() })
  output$download_gnomad_plot <- exportPlot("gnomAD_Frequency_Plot", gnomad_plot_obj)
  output$download_gnomad_table <- exportTable("gnomAD_Frequency_Data", gnomad_data_for_table)
  
  pathway_oncoplot_clinical_obj <- reactive({
    req(maf_object(), input$pathway_oncoplot_top_n)
    clin_features <- if (is.null(input$pathway_clinical_features) || input$pathway_clinical_features == "") NULL else input$pathway_clinical_features
    oncoplot(maf = maf_object(), top = input$pathway_oncoplot_top_n, pathways = "sigpw", 
             clinicalFeatures = clin_features, showTumorSampleBarcodes = input$pathway_show_samples,
             fontSize = 1.2, gene_mar = 8, legendFontSize = 1.4, annotationFontSize = 1.4)
  })
  output$pathway_oncoplot_clinical <- renderPlot({ pathway_oncoplot_clinical_obj() })
  output$download_pathway_oncoplot_clinical <- exportPlot("Pathway_Oncoplot_Clinical", pathway_oncoplot_clinical_obj)
  
  output$dynamic_filter_ui <- renderUI({
    req(maf_object(), input$clinVar)
    choices <- maf_object()@clinical.data %>%
      pull(!!sym(input$clinVar)) %>%
      unique() %>%
      na.omit()
    choices <- choices[choices != ""]
    selectizeInput("dynamic_filter_values", "Filter by observations:", 
                   choices = choices, selected = choices, multiple = TRUE)
  })
  exploratory_data_reactive <- reactive({
    req(maf_object(), input$clinVar, input$dynamic_filter_values)
    clin_data_filtered <- as.data.frame(maf_object()@clinical.data) %>%
      filter(!!sym(input$clinVar) %in% input$dynamic_filter_values)
    samples_to_keep <- clin_data_filtered$Tumor_Sample_Barcode
    maf_filtered <- subsetMaf(maf = maf_object(), tsb = samples_to_keep, mafObj = TRUE)
    if (getSampleSummary(maf_filtered)[,.N] < 3) return(NULL)
    mut_matrix <- t(maf2matrix(maf_filtered))
    if (ncol(mut_matrix) < 2) return(NULL)
    set.seed(123)
    if(input$multiMethod == "PCA"){
      res <- prcomp(mut_matrix, scale. = TRUE)
      res_df <- as.data.frame(res$x)
      res_df$Tumor_Sample_Barcode <- rownames(mut_matrix)
      res_df <- left_join(res_df, clin_data_filtered, by="Tumor_Sample_Barcode")
    } else if(input$multiMethod == "t-SNE"){
      res <- Rtsne(mut_matrix, perplexity=min(3, floor((nrow(mut_matrix)-1)/3)))
      res_df <- as.data.frame(res$Y)
      colnames(res_df) <- c("tSNE1", "tSNE2")
      res_df$Tumor_Sample_Barcode <- rownames(mut_matrix)
      res_df <- left_join(res_df, clin_data_filtered, by="Tumor_Sample_Barcode")
    } else if(input$multiMethod == "UMAP"){
      res <- umap(mut_matrix)
      res_df <- as.data.frame(res$layout)
      colnames(res_df) <- c("UMAP1", "UMAP2")
      res_df$Tumor_Sample_Barcode <- rownames(mut_matrix)
      res_df <- left_join(res_df, clin_data_filtered, by="Tumor_Sample_Barcode")
    }
    return(res_df)
  })
  multiPlot_obj <- reactive({
    df <- exploratory_data_reactive()
    validate(need(is.data.frame(df), "Insufficient data for analysis.\n(Requires >2 samples and >1 mutated gene in the filter)"))
    
    if(input$multiMethod == "PCA"){
      p <- ggplot(df, aes(x = PC1, y = PC2, color = .data[[input$clinVar]]))
      x_label <- "Principal Component 1"
      y_label <- "Principal Component 2"
    } else if(input$multiMethod == "t-SNE"){
      p <- ggplot(df, aes(x = tSNE1, y = tSNE2, color = .data[[input$clinVar]]))
      x_label <- "t-SNE 1"
      y_label <- "t-SNE 2"
    } else if(input$multiMethod == "UMAP"){
      p <- ggplot(df, aes(x = UMAP1, y = UMAP2, color = .data[[input$clinVar]]))
      x_label <- "UMAP 1"
      y_label <- "UMAP 2"
    }
    p + geom_point(size = 4, alpha = 0.8) +
      labs(title = paste(input$multiMethod, "of samples based on mutation profile"),
           x = x_label, y = y_label, color = input$clinVar) +
      theme_minimal(base_size = 18) +
      guides(color = guide_legend(override.aes = list(size=5)))
  })
  output$multiPlot <- renderPlot({ multiPlot_obj() })
  output$downloadMultiPlot <- exportPlot("Exploratory_Analysis_Plot", multiPlot_obj)
  
  covBoxPlot_reactive <- reactive({
    req(coverage_data())
    ggplot(coverage_data(), aes(x=reorder(Sample, MeanCoverage), y=MeanCoverage)) +
      geom_boxplot(fill="steelblue") +
      theme_minimal(base_size = 18) +
      ylab("Mean Coverage") + xlab("Sample") +
      coord_flip()
  })
  output$covBoxPlot <- renderPlot({ covBoxPlot_reactive() })
  output$downloadCovBox <- exportPlot("Coverage_Boxplot", covBoxPlot_reactive)
}

# ==============================================
# 5. EXECUTE APPLICATION
# ==============================================
shinyApp(ui = ui, server = server)