if (!require("rtpcr")) {install.packages("rtpcr")
library(rtpcr)}
library(shiny)
library(multcompView)
library(dplyr)
library(ggplot2)
library(emmeans)
library(multcomp)
library(lme4)
library(lmerTest)

# UI
ui <- fluidPage(
  
  
  tags$head(
    tags$link(rel = "stylesheet", href = "https://fonts.googleapis.com/css2?family=Roboto:wght@300;400;700&display=swap"),
    tags$style(HTML("
      body { font-family: 'Roboto', sans-serif; font-size: 16px; line-height: 1.7;  color: #333; }
      .intro-container { padding: 20px; }
      h2 { font-size: 30px; font-weight: 700; color: #2c3e50; }
      h3 { font-size: 24px; font-weight: 700; color: #2c3e50; margin-top: 30px; }
      h4 { font-size: 20px; font-weight: 400; }
      .info-box { background-color: #f0f7ff; border-left: 4px solid #2a7df6; padding: 20px; margin: 25px 0; border-radius: 4px; font-size: 15.5px; }
      .info-box2 { background-color: #f0ffff; border: 1px solid #2a7dff; padding: 20px; border-radius: 4px; font-size: 15.5px; }
      .info-box3 { border: 2px solid #fff; padding: 20px; border-radius: 4px; font-size: 15.5px; }
      .info-box strong { color: #2a7df6; }
      .info-box br { margin-bottom: 10px; content: ' '; display: block; }
      .note-title { font-weight: bold; color: #2a7df6; text-transform: uppercase; font-size: 0.85em; }
      .img-container { text-align: center; margin: 35px 0; border: 1px solid #EE3B3B; padding: 15px; background: #fff; border-radius: 8px; }  
      .logo-container { text-align: center; border-radius: 8px; }
      .txt-container { margin: 35px 0; border: 1px solid #ddd; padding: 15px; background: #fff; border-radius: 8px; }
      .fig-caption { font-style: italic; color: #666; margin-top: 12px; display: block; font-size: 14px; }
      .section-divider { border-top: 1px solid #C1CDCD; margin: 20px 0; }
    "))
),
  
  
  tags$script(HTML("
    $(document).on('keyup', function(e) {
      if (e.which == 13) {
        var activeTab = $('#func_tabs .active a').attr('data-value');
        var btnId = '';
        if (activeTab === 'meanTech') btnId = '#run_mt';
        if (activeTab === 'ANOVA_DCt') btnId = '#run_dc';
        if (activeTab === 'ANOVA_DDCt') btnId = '#run_ddct';
        if (activeTab === 'Efficiency') btnId = '#run_eff';
        if (activeTab === 'TTEST_DDCt') btnId = '#run_tt';
        if (activeTab === 'WILCOX_DDCt') btnId = '#run_wx';
        if (activeTab === 'plotFactor') btnId = '#run_pf';
        
        if (btnId !== '') {
          $(btnId).click();
        }
      }
    });
  ")),

  
  titlePanel('rtpcr: qPCR Data Analysis & Plotting'),
  fluidRow(column(12,
           div(style = "margin-bottom: 20px; color: #333; line-height: 1.6;", 
               p("Welcome! This is the shiny version of the rtpcr package, a web application developed using R/Shiny for comparative (dCt and ddCt) analysis of qPCR data.")))),
  br(),
  
  sidebarLayout(
    sidebarPanel(
      tabsetPanel(id = "func_tabs",
                  tabPanel("Home", value = "Home",
                           br(),
                           div(class = "info-box3", h4("Citation"), "Ghader Mirzaghaderi (2025). rtpcr: a package for statistical analysis 
                               and graphical presentation of qPCR data in R. PeerJ 13:e20185.", 
                               a(href="https://doi.org/10.7717/peerj.20185", "doi.org/10.7717/peerj.20185")),
                           
                           ),

                  
                  
                  
                  
                  
                  tabPanel("ANOVA_DCt", value = "ANOVA_DCt",
                           radioButtons("src_dc", "Data Source", 
                                        choices = c("Upload CSV" = "user", "Sample Data" = "sample", "Result: meanTech" = "res_mt"), 
                                        selected = "user"),
                           conditionalPanel("input.src_dc == 'user'", fileInput("file_dc", "Upload CSV", accept = ".csv")),
                           numericInput("numFactors_dc", "Number of factors", 1, min = 1),
                           numericInput("numRefGenes_dc", "Number of reference genes", 1, min = 1),
                           textInput("block_dc", "Block column name (if present)", ""),
                           textInput("model_dc", "model (optional; e.g., model = wDCt ~ factorA * factorB)", ""),
                           checkboxInput("set40_dc", "Set missing Ct to 40", FALSE),
                           checkboxInput("setModelse_dc", "Model_based se", TRUE),
                           selectInput("pAdj_dc", "p-value adjustment", choices = c("none","holm","bonferroni","fdr")),
                           actionButton("run_dc", "Run ANOVA_DCt"),
                  ),
                  
                  tabPanel("ANOVA_DDCt", value = "ANOVA_DDCt",
                           radioButtons("src_ddct", "Data Source", 
                                        choices = c("Upload CSV" = "user", "Sample Data" = "sample", "Result: meanTech" = "res_mt"), 
                                        selected = "user"),
                           conditionalPanel("input.src_ddct == 'user'", fileInput("file_ddct", "Upload CSV", accept = ".csv")),
                           numericInput("numFactors_ddct", "Number of factors", 1, min = 1),
                           numericInput("numRefGenes_ddct", "Number of reference genes", 1, min = 1),
                           textInput("specs_ddct", "specs (e.g., A or A | B or A | B * C)", ""),
                           textInput("block_ddct", "Block column name (if present)", ""),
                           textInput("calibrator_ddct", "Calibrator Level (Optional; default is the first level)", ""),
                           selectInput("seType_ddct", "SE type", choices = c("single.group",  "two.group", "paired.group")),
                           textInput("model_ddct", "model (optional; e.g., model = wDCt ~ factorA * factorB)", ""),
                           selectInput("pAdj_ddct", "p-value adjustment", choices = c("none","holm","bonferroni","fdr")),
                           checkboxInput("set40_ddct", "Set missing Ct to 40", FALSE),
                           checkboxInput("setModelse_ddct", "Model_based se", TRUE),
                           actionButton("run_ddct", "Run ANOVA_DDCt")
                  ),
                  
                  tabPanel("Efficiency", value = "Efficiency",
                           radioButtons("src_eff", "Data Source",choices = c("Upload CSV" = "user", "Sample Data" = "sample"), selected = "user"),
                           conditionalPanel("input.src_eff == 'user'", fileInput("file_eff", "Upload CSV", accept = ".csv")),
                           numericInput("baseSize_eff", "Base font size", 12, min = 6),
                           textInput("extra_eff", "Further ggplot layers (Optional; e.g., + ylab('Ct value'))", ""),
                           hr(),
                           h4("Plot size"),
                           fluidRow(
                             column(6, numericInput("eff_w", "Width (in)", 8, min = 1)),
                             column(6, numericInput("eff_h", "Height (in)", 6, min = 1))
                           ),
                           actionButton("run_eff", "Run Efficiency")
                  ),
                  
                  tabPanel("TTEST_DDCt", value = "TTEST_DDCt",
                           radioButtons("src_tt", "Data Source", 
                                        choices = c("Upload CSV" = "user", "Sample Data" = "sample", "Result: meanTech" = "res_mt"), 
                                        selected = "user"),
                           conditionalPanel("input.src_tt == 'user'", fileInput("file_tt", "Upload CSV", accept = ".csv")),
                           numericInput("numRefGenes_tt", "Number of reference genes", 1),
                           textInput("factorLevels_tt", "Factor levels (comma separated)", ""),
                           checkboxInput("paired_tt", "Paired test", FALSE),
                           checkboxInput("equalVar_tt", "Equal variance", TRUE),
                           selectInput("pAdj_tt", "p-value adjustment", choices = c("none","holm","bonferroni","fdr")),
                           checkboxInput("set40_tt", "Set missing Ct to 40", FALSE),
                           actionButton("run_tt", "Run TTEST_DDCt")
                  ),
                  
                  tabPanel("WILCOX_DDCt", value = "WILCOX_DDCt",
                           radioButtons("src_wx", "Data Source", 
                                        choices = c("Upload CSV" = "user", "Sample Data" = "sample", "Result: meanTech" = "res_mt"), 
                                        selected = "user"),
                           conditionalPanel("input.src_wx == 'user'", fileInput("file_wx", "Upload CSV", accept = ".csv")),
                           numericInput("numRefGenes_wx", "Number of reference genes", 1),
                           textInput("factorLevels_wx", "Factor levels (comma separated)", ""),
                           checkboxInput("paired_wx", "Paired test", FALSE),
                           selectInput("pAdj_wx", "p-value adjustment", choices = c("none","holm","bonferroni","fdr")),
                           checkboxInput("set40_wx", "Set missing Ct to 40", FALSE),
                           actionButton("run_wx", "Run WILCOX_DDCt")
                  ),
                  
                  tabPanel("plotFactor", value = "plotFactor",
                           radioButtons("src_pf", "Data Source", 
                                        choices = c("Upload CSV" = "user", 
                                                    "Sample Data" = "sample",
                                                    "Result: ANOVA_DCt" = "res_dc",
                                                    "Result: ANOVA_DDCt" = "res_ddct",
                                                    "Result: TTEST_DDCt" = "res_tt",
                                                    "Result: WILCOX_DDCt" = "res_wx"), 
                                        selected = "user"),
                           conditionalPanel("input.src_pf == 'user'", fileInput("file_pf", "Upload CSV", accept = ".csv")),
                           selectInput("pf_x", "X Axis Column", choices = NULL),
                           selectInput("pf_y", "Y Axis Column", choices = NULL),
                           selectInput("pf_low", "Lower Error Column", choices = NULL),
                           selectInput("pf_up", "Upper Error Column", choices = NULL),
                           selectInput("pf_group", "Group (Fill) Column", choices = c("None" = "")),
                           selectInput("pf_facet", "Facet Column", choices = c("None" = "")),
                           selectInput("pf_letters", "Letters Column (Optional)", choices = c("None" = "")),
                           numericInput("pf_letters_d", "Letters vertical offset", 0.2, step = 0.05),
                           numericInput("pf_col_width", "Column width", 0.8, step = 0.05),
                           numericInput("pf_err_width", "Error bar width", 0.15, step = 0.05),
                           numericInput("pf_dodge_width", "Dodge width", 0.8, step = 0.05),
                           textInput("pf_fill_colors", "Fill colors (comma separated)", ""),
                           textInput("pf_color", "Outline color", "black"),
                           sliderInput("pf_alpha", "Transparency (alpha)", min = 0, max = 1, value = 0.5, step = 0.05),
                           numericInput("pf_base_size", "Base font size", 12, min = 6),
                           checkboxInput("pf_legend_none", "Hide legend", FALSE),
                           checkboxInput("pf_removeRows", "Remove Calibrator Cols", FALSE),
                           checkboxInput("pf_removeText", "Remove Calibrator Text", FALSE),
                           fluidRow(column(6, numericInput("pf_legend_x", "Legend X", 0.85, min = 0, max = 1, step = 0.05) ),
                                    column(6, numericInput("pf_legend_y", "Legend Y", 0.85, min = 0, max = 1, step = 0.05) ) ),
                           fluidRow(column(6, numericInput("pf_w", "Width (in)", 8, min = 1)),
                                    column(6, numericInput("pf_h", "Height (in)", 9, min = 1))),
                           textInput("extra_pf", "Further ggplot layers (Optional; e.g., + ylab('Fold change'))", ""),
                           actionButton("run_pf", "Run plotFactor")
                  ),
                  tabPanel("meanTech", value = "meanTech",
                           p("Only apply this function if your dataset requires averaging across technical replicates before further analysis.", 
                             style = "color: #555; font-style: italic; margin-top: 10px; margin-bottom: 15px;"),
                           
                           radioButtons("src_mt", "Data Source", choices = c("Upload CSV" = "user", "Sample Data" = "sample"), selected = "user"),
                           conditionalPanel("input.src_mt == 'user'", fileInput("file_mt", "Upload CSV", accept = ".csv")),
                           numericInput("groups_mt", "Groups (Number of grouping columns)", 3, min = 1),
                           numericInput("numFactors_mt", "Number of factors", 2, min = 1),
                           numericInput("numRefGenes_mt", "Number of reference genes", 1, min = 1),
                           textInput("block_mt", "Block column (optional)", ""),
                           checkboxInput("set40_mt", "Set missing Ct to 40", FALSE),
                           actionButton("run_mt", "Run meanTech")
                  )
      )
    ),
    mainPanel(
      tabsetPanel(id = "output_tabs",
                  tabPanel("Introduction", value = "Introduction",
                           div(style = "padding: 20px; line-height: 1.6;",
                               h2("Getting Started with rtpcr"),
                               p("The rtpcr package facilitates relative expression analysis using delta Ct (dCt) and delta delta Ct (ddCt) methods. It supports t-test, ANOVA, and publication-ready visualizations. The package implements a general calculation method adopted from Ganger et al. (2017) and Taylor et al. (2019), covering both the Livak and Pfaffl methods."),
                               
                               div(class = "section-divider"),
                               
            
                              h2("rtpcr Funtions"),
                              p("In the rtpcr package, functions with _DDCt at the end of their name (ANOVA_DDCt, TTEST_DDCt, WILCOX_DDCt) perform expression analysis based on the delta delta Ct (ddCt)
                               method, while ANOVA_DCt function analyze gene expression using the delta Ct (dCt) method. The ANOVA prefix in the function name means that the function uses analysis of variance (using a default full factorial
                               model or a user defined model) for statistical analysis, and mean comparisons. Mean comparisons is performed by the emmeans function using the model resulting from the ANOVA analysis."),

                              helpText(div(class = "info-box2",
                                           tags$strong("ANOVA_DCt: "), " dCt expression analysis for all the level combinations of factor(s).", 
                                           br(),
                                             tags$strong("ANOVA_DDCt: "), " ddCt expression analysis for levels of a factor (geneally or per levels of another factors(s)), specified by the `specs` argument. ", tags$br(),
                                             tags$strong("TTEST_DDCt: "), " ddCt method t.test analysis for paired or unpaired samples.", tags$br(),
                                             tags$strong("WILCOX_DDCt: "), " ddCt method wilcox.test analysis for paired or unpaired samples.", tags$br(),
                                             tags$strong("plotFactor: "), " Bar plot of gene expression", tags$br(),
                                             tags$strong("efficiency: "), " Amplification efficiency statistics and standard curves", tags$br(),
                                             tags$strong("meanTech: "), " Calculate mean of technical replicates. This is used if your data needs averaging over biological replicates. ")),

                              

                               h3("1. Input Data Structure"),
                               p("For relative expression analysis, use TTEST_DDCt, WILCOX_DDCt, ANOVA_DCt, and ANOVA_DDCt functions from the left side panel. 
                               You need to prepare your input data based on the experimenta design, amplification efficiency 
                                 (E) and Ct values (the mean of technical replicates) which should include the following columns from left to 
                                wright:"),
                               tags$ol(
                                 tags$li("Experimental condition columns (and one block column if you want the variation from different plates be concidered by the nodel)."),
                                 tags$li("Replicate information (biological replicates or subjects)."),
                                 tags$li("Target genes efficiency (E) and Ct values (paired columns)."),
                                 tags$li("Reference genes efficiency (E) and Ct values (paired columns).")
                               ),
                               p("If the E values are not available you should use 
                                 ‘2’ instead representing the complete primer amplification efficiency. The package supports one or more target or reference gene(s),
                               supplied as efficiency–Ct column pairs. Reference gene columns must
                               always appear last. Two sample input data sets are presented below."),
                               
                               div(class = "img-container",
                                   img(src = "sampleData1.png", style = "width: 70%; max-width: 600px;"),
                                   span(class = "fig-caption", tags$strong("Figure 1: "), "A sample input data with one experimetal factor, replicate column and E/Ct information of target and reference genes.")),
                               
                               p("If there is no blocking factor, omit that column. However, a replicate column (e.g., 'Rep' or 'id') is always required."),
                               
                               div(class = "img-container",
                                   img(src = "dataStructure1.png", style = "width: 100%; max-width: 800px;"),
                                   span(class = "fig-caption", tags$strong("Figure 2: "), "A sample input data with two experimental factors, blocking factor, replicate column and E/Ct information of target and reference genes.")),

                               
                               h3("2. Notes"),
                               div(class = "info-box",
                                   tags$strong("Note 1: Blocking:"), " When a qPCR experiment is done in multiple qPCR plates, variation resulting from the plates may interfere with the actual amount of gene 
                                                                       expression. One solution is to conduct each plate as a randomized block so that at least one replicate of each treatment and control is present on a plate. Block effect is usually considered as random and its
                                                                       interaction with any main effect is not considered.", tags$br(),
                                   br(),
                                   tags$strong("Note 2: Biological Replicates:"), " For TTEST_DDCt and WILCOX_DDCt (independent groups), ANOVA_DCt,
                                                                       and ANOVA_DDCt each row is from a separate and unique biological replicate. For example, a data frame with 12 rows has come from an experiment with 12 individuals. The repeated measure models are intended for experiments with repeated observations (e.g. time-course data). In repeated measure experiments the Replicate column contains identifiers for each individual (id or subject). For example, all rows with a `1` at Rep column correspond to a single individual, all rows with a `2` correspond to another individual, and so on, which have been sampled at
                                                                       specific time points.", tags$br(),
                                   br(),
                                   tags$strong("Note 3: Technical Replicates:"), " Your data table may also include a column of technical replicates (For example, using one target and one reference genes, if you want to have 4 
                                                biological and 3 technical replicates under Control and Treatment conditions, there would be a table of 24 rows containing both biological 
                                                replicates and technical replicate columns in the data). In this case, the meanTech function should be applied first to calculate the mean of 
                                                the technical replicates. The resulting collapsed table is then used as the input for expression analysis. To use the meanTech function correctly, the technical replicate column must appear immediately after the biological replicate column (see the sample data of the meanTech function).", tags$br(),
                                   br(),
                                   tags$strong("Note 4: Efficiency (E):"), " Complete amplification efficiency (E) in the input data is denoted by
                                                2. This means that 2 indicates 100%, and 1.85 and 1.70 indicate 0.85%
                                                and 0.70% amplification efficiencies."),

                               
                               
                               
                               h3("3. Handling Missing Data"),
                               p("Missing Ct values for target genes is Handled using the `set_missing_target_Ct_to_40` function. If TRUE, missing target gene
                                  Ct values become 40; if FALSE (default), they become NA. missing Ct values of reference genes are always converted to NA. If there are more
                                  than one reference gene, NA in the place of the E or the Ct value cause skipping that gene and remaining references are geometrically averaged in that replicate."),
                               
                               div(class = "section-divider"),
                               
                               p("Contact email: gh.mirzaghaderi at uok.ac.ir"),
                               p("For further details please visit:",
                                 a(href="https://github.com/mirzaghaderi/rtpcr", "https://github.com/mirzaghaderi/rtpcr")),
                               )
                  ),
                  tabPanel("ANOVA_DCt", value = "ANOVA_DCt",
                           br(),
                           tabsetPanel(id = "sub_dc",
                                       tabPanel("Input data", tableOutput("preview_dc")),
                                       tabPanel("Relative Expression", br(), downloadButton("download_dc","Download"), br(), tableOutput("table_dc")),
                                       tabPanel("Per-Gene Results",
                                                selectInput("gene_dc", "Select gene:", choices = character(0)),
                                                tabsetPanel(
                                                  tabPanel("ANOVA, Normality & Singularity", verbatimTextOutput("singular_dc"),
                                                           p("NOTE ", style = "color: #EE3B3B;"),
                                                           p("If the singularity of the model is TRUE, the statistical results including ANOVA and significance is not reliable !",
                                                           style = "color: #9C9C9C; font-style: italic; margin-top: 10px; margin-bottom: 15px;")),
                                                  tabPanel("LM Object", br(), downloadButton("download_lm_dc", "Download LM (.rds)")),
                                                  tabPanel("Final table", br(), downloadButton("download_final_dc", "Download Final Table"), br(), tableOutput("final_data_dc"))
                                                ))
                           )
                  ),
                  tabPanel("ANOVA_DDCt", value = "ANOVA_DDCt",
                           br(),
                           tabsetPanel(id = "sub_ddct",
                                       tabPanel("Input data", tableOutput("preview_ddct")),
                                       tabPanel("Relative Expression", br(), downloadButton("download_ddct", "Download"), br(), tableOutput("table_ddct")),
                                       tabPanel("Per-Gene Results",
                                                selectInput("gene_ddct", "Select gene:", choices = character(0)),
                                                tabsetPanel(
                                                  tabPanel("ANOVA, Normality & Singularity", verbatimTextOutput("singular_ddct"), 
                                                           p("NOTE ", style = "color: #EE3B3B;"),
                                                           p("If the singularity of the model is TRUE, the statistical results including ANOVA and significance is not reliable !",
                                                           style = "color: #9C9C9C; font-style: italic; margin-top: 10px; margin-bottom: 15px;")),
                                                  tabPanel("LM Object", br(), downloadButton("download_lm_ddct", "Download LM (.rds)")),
                                                  tabPanel("Final table", br(), downloadButton("download_final_ddct", "Download Final Table"), br(), tableOutput("final_data_ddct"))
                                                ))
                           )
                  ),
                  tabPanel("Efficiency", value = "Efficiency",
                           br(),
                           tabsetPanel(id = "sub_eff",
                                       tabPanel("Input data", tableOutput("preview_eff")),
                                       tabPanel("Values & Statistics", br(), downloadButton("download_eff", "Download CSV"), br(), verbatimTextOutput("table_eff")),
                                       tabPanel("Plot", br(), downloadButton("download_eff_png", "PNG"), downloadButton("download_eff_pdf", "PDF"), hr(), plotOutput("plot_eff"))
                           )
                  ),
                  tabPanel("TTEST_DDCt", value = "TTEST_DDCt",
                           br(),
                           tabsetPanel(id = "sub_tt",
                                       tabPanel("Input Data", tableOutput("preview_tt")),
                                       tabPanel("Results", br(), downloadButton("download_tt", "Download CSV"), br(), tableOutput("table_tt"))
                           )
                  ),
                  tabPanel("WILCOX_DDCt", value = "WILCOX_DDCt",
                           br(),
                           tabsetPanel(id = "sub_wx",
                                       tabPanel("Input Data", tableOutput("preview_wx") ),
                                       tabPanel("Results", br(), downloadButton("download_wx", "Download CSV"), br(), tableOutput("table_wx"))
                           )
                  ),
                  tabPanel("plotFactor", value = "plotFactor",
                           br(),
                           tabsetPanel(id = "sub_pf",
                                       tabPanel("Input Data", tableOutput("preview_pf")),
                                       tabPanel("Plot", br(), downloadButton("download_pf_plot", "PNG"), downloadButton("download_pf_pdf", "PDF"), hr(), plotOutput("plot_pf_main"))
                           )
                  ),
                  tabPanel("meanTech", value = "meanTech",
                           br(),
                           tabsetPanel(id = "sub_mt",
                                       tabPanel("Input data", tableOutput("preview_mt")),
                                       tabPanel("Results", br(), downloadButton("download_mt", "Download Results"), br(), tableOutput("table_mt"))
                           )
                  )
      )
    )
  )
)

# SERVER
server <- function(input, output, session) {
  
  # Helper for data sources
  get_data <- function(src_input, file_input, sample_path, mean_res = NULL) {
    if (src_input == "sample") {
      return(read.csv(sample_path))
    } else if (src_input == "res_mt") {
      req(mean_res)
      return(mean_res)
    } else {
      req(file_input)
      return(read.csv(file_input$datapath))
    }
  }
  
  # meanTech Logic
  data_mt <- reactive({ get_data(input$src_mt, input$file_mt, "mean/data_withTechRep.csv") })
  
  res_mt <- eventReactive(input$run_mt, {
    tryCatch({
      meanTech(data_mt(), groups = input$groups_mt, numOfFactors = input$numFactors_mt, 
               numberOfrefGenes = input$numRefGenes_mt, block = if(input$block_mt=="") NULL else input$block_mt,
               set_missing_target_Ct_to_40 = input$set40_mt)
    }, error = function(e) { showNotification(e$message, type="error"); NULL })
  })
  
  output$preview_mt <- renderTable({ req(data_mt()); head(data_mt(), 50) })
  output$table_mt   <- renderTable({ req(res_mt()); head(res_mt(), 50) }, digits = 5)
  output$download_mt <- downloadHandler(filename = "meanTech_results.csv", content = function(f) write.csv(res_mt(), f, row.names=F))
  
  # Analysis Logic
  data_dc   <- reactive({ get_data(input$src_dc, input$file_dc, "exp/data_3factor.csv", res_mt()) })
  data_ddct <- reactive({ get_data(input$src_ddct, input$file_ddct, "exp/data_2factorBlock3ref.csv", res_mt()) })
  data_eff  <- reactive({ get_data(input$src_eff, input$file_eff, "eff/data_efficiency1.csv") })
  data_tt   <- reactive({ get_data(input$src_tt, input$file_tt, "exp/data_ttest18genes.csv", res_mt()) })
  data_wx   <- reactive({ get_data(input$src_wx, input$file_wx, "exp/data_ttest18genes.csv", res_mt()) })
  
  output$preview_dc <- renderTable({ req(data_dc()); head(data_dc(), 50) })
  output$preview_ddct <- renderTable({ req(data_ddct()); head(data_ddct(), 50) })
  output$preview_eff <- renderTable({ req(data_eff()); head(data_eff(), 50) })
  output$preview_tt <- renderTable({ req(data_tt()); head(data_tt(), 50) })
  output$preview_wx <- renderTable({ req(data_wx()); head(data_wx(), 50) })
  
  # Sample updates
  observeEvent(input$src_mt, {
    if (input$src_mt == "sample") {
      updateNumericInput(session, "groups_mt", value = 2)
      updateNumericInput(session, "numFactors_mt", value = 1)
      updateNumericInput(session, "numRefGenes_mt", value = 1)
    }
  })
  observeEvent(input$src_dc, {
    if (input$src_dc == "sample") { updateNumericInput(session, "numFactors_dc", value = 3); updateNumericInput(session, "numRefGenes_dc", value = 1) }
  })
  observeEvent(input$src_ddct, {
    if (input$src_ddct == "sample") { updateNumericInput(session, "numFactors_ddct", value = 2); updateNumericInput(session, "numRefGenes_ddct", value = 3); updateTextInput(session, "specs_ddct", value = "Concentration"); updateTextInput(session, "block_ddct", value = "block") }
  })
  
  # Event Handlers
  observeEvent(input$run_mt,   { updateTabsetPanel(session, "sub_mt", selected = "Results") })
  observeEvent(input$run_dc,   { updateTabsetPanel(session, "sub_dc", selected = "Relative Expression") })
  observeEvent(input$run_ddct, { updateTabsetPanel(session, "sub_ddct", selected = "Relative Expression") })
  observeEvent(input$run_eff,  { updateTabsetPanel(session, "sub_eff", selected = "Table") })
  observeEvent(input$run_tt,   { updateTabsetPanel(session, "sub_tt", selected = "Results") })
  observeEvent(input$run_wx,   { updateTabsetPanel(session, "sub_wx", selected = "Results") })
  observeEvent(input$run_pf,   { updateTabsetPanel(session, "sub_pf", selected = "Plot") })
  
  res_dc <- eventReactive(input$run_dc, {
    tryCatch({ ANOVA_DCt(data_dc(), numOfFactors = input$numFactors_dc, numberOfrefGenes = input$numRefGenes_dc,
                         block = if(input$block_dc == "") NULL else input$block_dc, 
                         model = if(input$model_dc == "") NULL else input$model_dc,
                         p.adj = input$pAdj_dc, set_missing_target_Ct_to_40 = input$set40_dc,
                         modelBased_se = input$setModelse_dc)
    }, error = function(e) { showNotification(e$message, type="error"); NULL })
  })
  
  res_ddct <- eventReactive(input$run_ddct, {
    tryCatch({ ANOVA_DDCt(data_ddct(), numOfFactors = input$numFactors_ddct, numberOfrefGenes = input$numRefGenes_ddct,
                          specs = input$specs_ddct, block = if(input$block_ddct == "") NULL else input$block_ddct,
                          calibratorLevel = if(input$calibrator_ddct == "") NULL else input$calibrator_ddct,
                          p.adj = input$pAdj_ddct, set_missing_target_Ct_to_40 = input$set40_ddct, 
                          model = if(input$model_ddct == "") NULL else input$model_ddct,
                          se.type = input$seType_ddct, modelBased_se = input$setModelse_ddct)
    }, error = function(e) { showNotification(e$message, type="error"); NULL })
  })
  
  res_eff <- eventReactive(input$run_eff, {
    res <- efficiency(data_eff(), base_size = input$baseSize_eff)
    if (nzchar(trimws(input$extra_eff))) {
      tryCatch({ res$plot <- eval(parse(text = paste0("res$plot", input$extra_eff))) }, error = function(e) { showNotification(e$message, type = "error") })
    }
    res
  })
  
  res_tt <- eventReactive(input$run_tt, {
    TTEST_DDCt(data_tt(), numberOfrefGenes = input$numRefGenes_tt,
               Factor.level.order = if(input$factorLevels_tt == "") NULL else trimws(unlist(strsplit(input$factorLevels_tt, ","))),
               paired = input$paired_tt, var.equal = input$equalVar_tt, p.adj = input$pAdj_tt, 
               set_missing_target_Ct_to_40 = input$set40_tt)
  })
  
  res_wx <- eventReactive(input$run_wx, {
    tryCatch({ WILCOX_DDCt(data_wx(), numberOfrefGenes = input$numRefGenes_wx,
                           Factor.level.order = if(input$factorLevels_wx == "") NULL else trimws(unlist(strsplit(input$factorLevels_wx, ","))),
                           paired = input$paired_wx, p.adj = input$pAdj_wx, set_missing_target_Ct_to_40 = input$set40_wx)
    }, error = function(e) { showNotification(e$message, type = "error"); NULL })
  })
  
  # Gene selection updates
  observe({
    req(res_dc())
    gene_list_dc <- names(res_dc()$perGene)
    if(!is.null(gene_list_dc)) {
      updateSelectInput(session, "gene_dc", choices = gene_list_dc, selected = gene_list_dc[1])
    }
  })
  
  observe({
    req(res_ddct())
    gene_list_ddct <- names(res_ddct()$perGene)
    if(!is.null(gene_list_ddct)) {
      updateSelectInput(session, "gene_ddct", choices = gene_list_ddct, selected = gene_list_ddct[1])
    }
  })
  
  output$table_dc <- renderTable({ res_dc()$relativeExpression }, digits = 5)
  output$singular_dc <- renderPrint({ req(res_dc(), input$gene_dc); list(
    Linear_model = res_dc()$perGene[[input$gene_dc]]$lm_formula,
    ANOVA_table = res_dc()$perGene[[input$gene_dc]]$ANOVA_table,
    Normality_test_of_residuals = shapiro.test(res_dc()$perGene[[input$gene_dc]]$Final_data$residual),
    Singularity_check_of_model = res_dc()$perGene[[input$gene_dc]]$is_singular)})
  output$final_data_dc <- renderTable({ req(res_dc(), input$gene_dc); res_dc()$perGene[[input$gene_dc]]$Final_data }, digits = 5)
  
  output$table_ddct <- renderTable({ res_ddct()$relativeExpression }, digits = 5)
  output$singular_ddct <- renderPrint({ req(res_ddct(), input$gene_ddct); list(
    Linear_model = res_ddct()$perGene[[input$gene_ddct]]$lm_formula,
    res_ddct()$perGene[[input$gene_ddct]]$ANOVA_table,
    Normality_test_of_residuals = shapiro.test(res_ddct()$perGene[[input$gene_ddct]]$Final_data$residual),
    Singularity_check_of_model = res_ddct()$perGene[[input$gene_ddct]]$singular)})
  output$final_data_ddct <- renderTable({ req(res_ddct(), input$gene_ddct); res_ddct()$perGene[[input$gene_ddct]]$Final_data }, digits = 5)
  
  output$table_eff <- renderPrint({ list(Efficiency = res_eff()$Efficiency, 
                                         Slope_compare = res_eff()$Slope_compare) })
  output$plot_eff <- renderPlot({ res_eff()$plot }, width = function() input$eff_w * 72, height = function() input$eff_h * 72)
  output$table_tt <- renderTable({ res_tt() }, digits = 5)
  output$table_wx <- renderTable({ res_wx() }, digits = 5)
  
  # plotFactor Logic
  df_pf <- reactive({ 
    if (input$src_pf == "user") { req(input$file_pf); read.csv(input$file_pf$datapath) }
    else if (input$src_pf == "sample") { read.csv("plot/out_ANOVA_DDCt.csv") }
    else if (input$src_pf == "res_dc") { req(res_dc()); res_dc()$relativeExpression }
    else if (input$src_pf == "res_ddct") { req(res_ddct()); res_ddct()$relativeExpression }
    else if (input$src_pf == "res_tt") { req(res_tt()); res_tt() }
    else if (input$src_pf == "res_wx") { req(res_wx()); res_wx() }
  })
  
  observeEvent(df_pf(), {
    cols <- colnames(df_pf()); pick_col <- function(i) cols[min(i, length(cols))]
    updateSelectInput(session, "pf_x", choices = cols, selected = pick_col(2))
    updateSelectInput(session, "pf_y", choices = cols, selected = pick_col(4))
    updateSelectInput(session, "pf_low", choices = cols, selected = pick_col(9))
    updateSelectInput(session, "pf_up", choices = cols, selected = pick_col(10))
    updateSelectInput(session, "pf_group", choices = c("None" = "", cols), selected = if(input$src_pf=="sample") "gene" else "")
    updateSelectInput(session, "pf_facet", choices = c("None" = "", cols), selected = if(input$src_pf=="sample") "gene" else "")
    updateSelectInput(session, "pf_letters", choices = c("None" = "", cols), selected = if(input$src_pf=="sample") "sig" else "")
  })
  
  pf_plot_obj <- eventReactive(input$run_pf, {
    p <- plotFactor(data = df_pf(), x_col = input$pf_x, y_col = input$pf_y, Lower.se_col = input$pf_low,
                    Upper.se_col = input$pf_up, group_col = if(input$pf_group == "") NULL else input$pf_group,
                    facet_col = if(input$pf_facet == "") NULL else input$pf_facet,
                    letters_col = if(input$pf_letters == "") NULL else input$pf_letters,
                    letters_d = input$pf_letters_d, col_width = input$pf_col_width, err_width = input$pf_err_width,
                    dodge_width = input$pf_dodge_width, fill_colors = if (input$pf_fill_colors == "") NULL else trimws(unlist(strsplit(input$pf_fill_colors, ","))),
                    color = if (input$pf_color == "") NA else input$pf_color, 
                    alpha = input$pf_alpha, 
                    base_size = input$pf_base_size,
                    legend_position = if (input$pf_legend_none) "none" else c(input$pf_legend_x, input$pf_legend_y),
                    removeCalibratorCols = input$pf_removeRows, removeCalibratorText = input$pf_removeText)
    if (nzchar(trimws(input$extra_pf))) {
      tryCatch({ p <- eval(parse(text = paste0("p", input$extra_pf))) }, error = function(e) { showNotification(paste("Plot adjustment error:", e$message), type = "error") })
    }
    p
  })
  output$preview_pf <- renderTable({ req(df_pf()); head(df_pf(), 50) })
  output$plot_pf_main <- renderPlot({ req(pf_plot_obj()); pf_plot_obj() }, width = function() input$pf_w * 72, height = function() input$pf_h * 72)
  
  # Navigation logic
  observeEvent(input$func_tabs, {
    if (input$func_tabs == "Home") updateTabsetPanel(session, "output_tabs", selected = "Introduction")
    else updateTabsetPanel(session, "output_tabs", selected = input$func_tabs)
  })
  
  # Downloads
  output$download_dc <- downloadHandler(filename = "ANOVA_DCt.csv", content = function(f) write.csv(res_dc()$relativeExpression, f, row.names = FALSE))
  output$download_ddct <- downloadHandler(filename = "ANOVA_DDCt.csv", content = function(f) write.csv(res_ddct()$relativeExpression, f, row.names=F))
  output$download_eff <- downloadHandler(filename = "Efficiency.csv", content = function(f) write.csv(res_eff()$Efficiency, f, row.names=F))
  output$download_tt <- downloadHandler(filename = "TTEST.csv", content = function(f) write.csv(res_tt(), f, row.names=F))
  output$download_wx <- downloadHandler(filename = "WILCOX_DDCt.csv", content = function(f) write.csv(res_wx(), f, row.names = FALSE))
  output$download_final_dc <- downloadHandler(filename = function() paste0("Final_Table_DCt_", input$gene_dc, ".csv"), content = function(f) write.csv(res_dc()$perGene[[input$gene_dc]]$Final_data, f, row.names = FALSE))
  output$download_final_ddct <- downloadHandler(filename = function() paste0("Final_Table_DDCt_", input$gene_ddct, ".csv"), content = function(f) write.csv(res_ddct()$perGene[[input$gene_ddct]]$Final_data, f, row.names = FALSE))
  output$download_eff_png <- downloadHandler(filename = "Eff.png", content = function(f) ggsave(f, plot = res_eff()$plot, device = "png", width = input$eff_w, height = input$eff_h))
  output$download_eff_pdf <- downloadHandler(filename = "Eff.pdf", content = function(f) ggsave(f, plot = res_eff()$plot, device = "pdf", width = input$eff_w, height = input$eff_h))
  output$download_pf_plot <- downloadHandler(filename = "qPCR_Plot.png", content = function(f) ggsave(f, plot = pf_plot_obj(), device = "png", width = input$pf_w, height = input$pf_h))
  output$download_pf_pdf <- downloadHandler(filename = "qPCR_Plot.pdf", content = function(f) ggsave(f, plot = pf_plot_obj(), device = "pdf", width = input$pf_w, height = input$pf_h))
  output$download_lm_dc <- downloadHandler(filename = function() paste0("LM_", input$gene_dc, ".rds"), content = function(f) saveRDS(res_dc()$perGene[[input$gene_dc]]$lm, f))
  output$download_lm_ddct <- downloadHandler(filename = function() paste0("LM_", input$gene_ddct, ".rds"), content = function(f) saveRDS(res_ddct()$perGene[[input$gene_ddct]]$lm, f))
}

shinyApp(ui, server)