PreProcessing.R

#################### UI ##################
preProcessingUI <- function(id, prefix="") {
    ns <- NS(id)
    tagList(
        useShinyjs(),  # Include shinyjs to enable/disable UI elements
        fluidRow(
            hidden(column(width = 6, id = ns("pr_c1"),
                          h4("Add or delete data columns"),
                          fluidRow(column(10, pickerInput(ns("remove_reps"), "Pick the samples you want to remove", 
                                                          choices = NULL, multiple = T,
                                                          options = list(
                                                              `live-search` = TRUE,
                                                              `actions-box` = FALSE,
                                                              `max-options` = -2, # Prevent selecting more than n-2 columns
                                                              `max-options-text` = "Cannot select more than n-2 columns"
                                                          ))),
                                   column(2, actionBttn(ns("h_balancing"),
                                                        icon = icon("info-circle"),
                                                        style = "pill", 
                                                        color = "royal", size = "xs")
                                   )),
                          fluidRow(column(6,
                                          p(htmlOutput(ns("res_num_reps")), style = "text-color:#AA2222")
                          ),
                          column(6,
                                 checkboxInput(ns("add_na_columns"), "Fill with empty columns")
                          )
                          ),
                          h5("Summary:"),
                          htmlOutput(ns("ptable_summary"), style = "border:solid;border-width:1px;padding:10px;")
            )),
            hidden(column(width = 3, id = ns("pr_c2"),
                          h4("Data manipulation and adjustment"),
                          fluidRow(column(10, checkboxInput(ns("logtrafo"), "Is the data already log-transformed?", value = F)),
                                   column(2, actionBttn(ns("h_logtrafo"),
                                                        icon = icon("info-circle"),
                                                        style = "pill", 
                                                        color = "royal", size = "xs")
                                   )),
                          fluidRow(column(10, numericInput(ns("max_na"), label = "Maximum number of missing values per feature", 
                                                           min = 0, max = 0, step = 1, value = 100)),
                                   column(2, actionBttn(ns("h_max_na"),
                                                        icon = icon("info-circle"),
                                                        style = "pill", 
                                                        color = "royal", size = "xs")
                                   )),
                          fluidRow(column(10, selectInput(ns("normalization"), label = "Normalization method",
                                                          choices = c(None = "none", Median = "colMedians", "Mean" = "colMeans", "Cyclic LOESS (LIMMA)" = "cyclicloess"),
                                                          selected = "none")),
                                   column(2, actionBttn(ns("h_normalization"),
                                                        icon = icon("info-circle"),
                                                        style = "pill", 
                                                        color = "royal", size = "xs")
                                   )),
                          fluidRow(column(10, selectInput(ns("summarize"), label = "Summarize to id features", 
                                                          choices = c(None = "none", "By sum" = "colSums",
                                                                      "By mean" = "colMeans", "By median" = "colMedians", 
                                                                      "Robust median (medpolish)" = "medianPolish"
                                                          ))),
                                   column(2, actionBttn(ns("h_summarize"),
                                                        icon = icon("info-circle"),
                                                        style = "pill", 
                                                        color = "royal", size = "xs")
                                   )),
                          
                          fluidRow(column(10, 
                                          selectInput(ns("batch_correction_method"), "Correction method", 
                                                      choices = c("None", "limma", "Combat"), selected = "none") # New input for batch correction method
                          ),
                          column(2, actionBttn(ns("h_batch_effect"),
                                               icon = icon("info-circle"),
                                               style = "pill", 
                                               color = "royal", size = "xs")
                          )),
                          style = 'border-left: 1px solid; margin-bottom: 20px;'  # Add bottom margin to separate from next section
                          
            )),
            hidden(column(3, id = ns("pr_c3"),
                          h4("Proceed to interaction with apps"),
                          actionButton(ns("proceed_to_apps"), "Proceed"),
                          textOutput(ns("txt_proceed_apps")),
                          fluidRow(
                              column(10, 
                                     style = "text-align: left;",
                                     checkboxInput(ns("map2uniprot"), "Map identifiers to UniProt accession numbers?", value = FALSE),
                              ),
                              column(2, actionBttn(ns("h_mapping"),
                                                   icon = icon("info-circle"),
                                                   style = "pill", 
                                                   color = "royal", size = "xs")
                              )
                          ),
                          selectInput(ns("idtype"), 
                                      "Select identifier type to map to UniProt",
                                      choices = c("Generic gene names (will be mapped to human)" = "Gene_Name", 
                                                  "Ensembl gene ids" = "Ensembl", 
                                                  "GeneIDs (Entrez IDs)" = "GeneID"), 
                                      selected = "GeneID"),
                          style = 'border-left: 1px solid; margin-bottom: 20px;') 
            )),    
        hr(style="border:solid;border-width:1px;"),
        hidden(fluidRow(id = ns("pr_plots"),
                        fluidRow(
                            column(4, fluidRow(
                                column(10,
                                       style = "text-align: left;",
                                       radioButtons(ns("show_pca_labels"), 
                                                    "How to show labels", 
                                                    choices = c("Add labels" = TRUE, 
                                                                "No labels" = FALSE), 
                                                    selected = TRUE)),
                                column(2, actionBttn(ns("h_pcaplot"),
                                                     icon = icon("info-circle"),
                                                     style = "pill", 
                                                     color = "royal", size = "xs")
                                )
                            )
                            ),
                            column(4, fluidRow(
                                column(10, 
                                       style = "text-align: left;",
                                       checkboxInput(ns("scale_corrplot"), 
                                                     "Adjust color gradient to available/calculated correlation values",
                                       ),                                        
                                ), 
                                column(2, actionBttn(ns("h_corrplot"),
                                                     icon = icon("info-circle"),
                                                     style = "pill", 
                                                     color = "royal", size = "xs")
                                ))
                            ),
                            
                            column(4, fluidRow(
                                column(10, 
                                       style = "text-align: left;",
                                       radioButtons(ns("show_missing_rows"), 
                                                    "Configure plot for missing values", 
                                                    choices = c("Distribution of missing values over features" = TRUE, 
                                                                "Missing values per sample" = FALSE), 
                                                    selected = TRUE)),
                                column(2, actionBttn(ns("h_missingplot"),
                                                     icon = icon("info-circle"),
                                                     style = "pill", 
                                                     color = "royal", size = "xs")
                                ))
                            )
                        ),
                        fluidRow(column(4, 
                                        girafeOutput(ns("pca_plot"), height = "500px")
                        ),
                        column(4, 
                               plotOutput(ns("corrplot"), height = "500px")
                        ),
                        column(4, 
                               plotOutput(ns("missingplot"), height = "500px")
                        ),  style = 'margin-bottom: 20px;'
                        ),
                        fluidRow(column(4,
                                        downloadBttn(ns("download_pca"), label = "Download figure (pdf)")),
                                 column(4,
                                        downloadBttn(ns("download_corrplot"), label = "Download figure (pdf)")),
                                 column(4,
                                        downloadBttn(ns("download_missingplot"), label = "Download figure (pdf)")),
                                 
                        )
        )
        )
    )
}





#################### Server ##################
preProcessingServer <- function(id, parent, expDesign, log_operations, SM) {
    moduleServer(
        id,
        function(input, output, session) {
            
            process_table <- reactiveVal(NULL) 
            processed_table <- reactiveVal(NULL)
            other_cols <- reactiveVal(NULL)
            pexp_design <- reactiveVal(NULL)
            exp_design <- reactiveVal(NULL)
            next_tab <- reactiveVal(NULL)
            removed_cols <- reactiveVal(NULL)
            
            # Reactive values to store group, batch and replicate information
            group_info <- reactiveVal(NULL)
            batch_info <- reactiveVal(NULL)
            replicate_info <- reactiveVal(NULL)
            
            # Reactive values to added columns
            added_columns <- reactiveVal(c())
            
            # Register ALL of them under this module's namespace
            ns_id <- session$ns("preProcessing")  
            SM$register_vals(ns_id, list(
                # processed_table = processed_table,
                other_cols = other_cols,
                pexp_design = pexp_design,
                exp_design = exp_design,
                process_table = process_table,
                # next_tab = next_tab,
                removed_cols = removed_cols,
                group_info = group_info,
                batch_info = batch_info,
                replicate_info = replicate_info,
                added_columns = added_columns
            ))
            
            SM$set_input_reassert(session$ns("remove_reps"))
            SM$set_input_reassert(session$ns("batch_correction_method"))
            SM$set_input_reassert(session$ns("max_na"))
            
            
            init_data <- function() {
                if(!is.null(process_table())) {
                    # Initialize processed_table with 'id' and 'quant' columns from process_table
                    initial_data <- process_table()
                    id_column <- initial_data[, grep("id", sapply(initial_data, class)), drop = FALSE]
                    # Substitute NA values in id column
                    id_column[is.na(id_column)] <- "No_ID_Given"
                    id_column[id_column == ""] <- "No_ID_Given"
                    quant_columns <- initial_data[, grep("quant", sapply(initial_data, class)), drop = FALSE]
                    processed_table(cbind(id_column, quant_columns)) # Initialize processed_table
                    # Keep all other columns separate, only to merge when summarizing
                    # This is why we need to reload them
                    other_cols(initial_data[, !(colnames(initial_data) %in% c(colnames(id_column), colnames(quant_columns))), drop = FALSE])
                }
            }
            

            ##########################################################################
            ##########################################################################
            observeEvent(expDesign$next_tab(), {
                if (!is.null(expDesign$next_tab()))  {
                    
                    print("init preprocessing")
                    process_table(expDesign$process_table())
                    pexp_design(expDesign$pexp_design())
                    exp_design(expDesign$pexp_design())
                    
                    init_data()
                    
                    shinyjs::show("pr_c1")
                    shinyjs::show("pr_c2")
                    shinyjs::show("pr_c3")
                    shinyjs::show("pr_plots")
                    
                    id_column <- processed_table()[, 1]
                    quant_columns <- processed_table()[, -1]
                    
                    # Substitute NA values in id column
                    id_column[is.na(id_column)] <- "No_ID_Given"
                    id_column[id_column == ""] <- "No_ID_Given"
                    
                    
                    # Check if there are duplicate IDs and hide/show the summarize input accordingly
                    #print(id_column)
                    if (sum(duplicated(id_column)) > 0) {
                        shinyjs::show("summarize")  # Show the 'Summarize to id features' selector if duplicates are present
                    } else {
                        shinyjs::disable("summarize")  # Hide the 'Summarize to id features' selector if no duplicates
                        #shinyjs::disable("h_summarize") 
                    }
                    
                    
                    # Update pickerInput with new column names and set max-options dynamically
                    total_columns <- ncol(expDesign$pexp_design())
                    if (!SM$restoring()) 
                        updatePickerInput(session, "remove_reps", 
                                          choices = colnames(exp_design()),
                                          options = list(
                                              `max-options` = total_columns - 2, # Set max-options dynamically
                                              `max-options-text` = "Cannot select more than n-2 columns"
                                          ))
                    
                    # Initialize group, batch and replicate information
                    group_info(pexp_design()[1,]) # Assuming Group info is in the first row
                    replicate_info(pexp_design()[2,]) # Adjust index based on data structure
                    batch_info(pexp_design()[3,]) # Adjust index based on data structure
                    
                    # Check for different batches and batch sample number >= 2
                    if (length(unique(batch_info())) > 1 && min(table(batch_info()) >= 2)) {
                        shinyjs::enable("batch_correction_method")
                        updateSelectInput(session, "batch_correction_method")
                    } else {
                        shinyjs::disable("batch_correction_method")
                    }
                    
                    # Determine if the data has been log-transformed
                    if (isTRUE(SM$restoring())) return()   # don't clobber restored inputs
                    tlog <- log_operations()
                    if (max(quant_columns, na.rm = T) / min(quant_columns, na.rm = T) < 100 || min(quant_columns, na.rm = T) < 0) {
                        updateCheckboxInput(session, "logtrafo", value = TRUE)
                        tlog[["preprocess_take_log2"]] <- FALSE
                    } else {
                        tlog[["preprocess_take_log2"]] <- TRUE
                    }
                    updateNumericInput(session, "max_na",
                                       min = 0, max = ncol(quant_columns),
                                       value = ncol(quant_columns)
                    )
                    log_operations(tlog)
                }
            })
            
            # choices for remove_reps come from exp_design()
            observeEvent(exp_design(), {
                print("Update remove_reps list")
                cn <- colnames(exp_design()); req(length(cn) > 0)
                sel_remove <- isolate(input$remove_reps)
                sel_remove <- intersect(sel_remove %||% character(0), cn)
                updatePickerInput(session, "remove_reps",
                                  choices = cn,
                                  selected = sel_remove,
                                  options = list(
                                      `max-options`      = ncol(exp_design()) - 2,
                                      `max-options-text` = "Cannot select more than n-2 columns"
                                  )
                )
                updateNumericInput(session, "max_na",
                                   min = 0, max = ncol(exp_design())
                )
            })
            
            # enable/disable batch control button
            observe({
                if (length(unique(batch_info())) > 1 && min(table(batch_info()) >= 2)) {
                    shinyjs::enable("batch_correction_method")
                } else {
                    shinyjs::disable("batch_correction_method")
                }
            })
            
            observeEvent(input$remove_reps, {
                req(processed_table())
                removed_cols(input$remove_reps)
            }, ignoreNULL = FALSE)
            
            ##########################################################################
            ##########################################################################
            # Add or delete data columns
            remove_cols <- function(tdata, remove_reps, texpdesign) {
                print("Removing columns...")
                withProgress(message = "Processing...", value = 0, min = 0, max = 1, {
                    cnames <- colnames(texpdesign)
                
                    # Removing selected columns
                    rem <- -which(names(tdata) %in% remove_reps)
                    if (length(rem) > 0) {
                        incProgress(0.1, detail = "Removing replicates")
                        tdata <- tdata[, rem]
                        rem2 <- -which(cnames %in% remove_reps)
                        pexp_design(texpdesign[, rem2])
                        cnames <- cnames[rem2]
                        tlog <- log_operations()
                        tlog[["preprocess_removed_replicates"]] <- remove_reps
                        log_operations(tlog)
                    } else {
                        pexp_design(texpdesign)
                    }
                    
                    # Update processed_table to keep only 'id' and 'quant' columns
                    id_column <- tdata[, grep("id", sapply(tdata, class)), drop = FALSE]
                    quant_columns <- tdata[, grep("quant", sapply(tdata, class)), drop = FALSE]
                    processed_table(cbind(id_column, quant_columns))
                })
            }
            
            
            ##########################################################################
            ##########################################################################
            ## Check for balanced exp. design
            check_balance <- function(texpdesign) {
                print("check for balancing")

                # Calculate the number of replicates per experimental condition
                ed_stats <- as.vector(table(texpdesign[1, ]))
                
                # Check if all experimental conditions have the same number of replicates
                if (length(unique(ed_stats)) > 1) {
                    # If the design is unbalanced, disable the Proceed button
                    disable("proceed_to_apps")
                    enable("add_na_columns")  # Enable fill with empty columns when unbalanced
                    tout <- paste(
                        "This unbalanced design has between ",
                        min(ed_stats),
                        " and maximally", max(ed_stats),
                        "replicates for each experimental condition (sample type). <br/><b>Please click checkbox to make data balanced.</b>"
                    )
                    tout2 <- paste("Unbalanced design!")
                } else {
                    # If the design is balanced, enable the Proceed button
                    enable("proceed_to_apps")
                    #disable("add_na_columns")  # Disable fill with empty columns when balanced
                    tout <- paste("Experimental design is balanced.")
                    tout2 <- NULL
                }
                
                output$res_num_reps <- renderText(tout)
                output$txt_proceed_apps <- renderText(tout2)
            }
            
            
            
            ##########################################################################
            ##########################################################################
            # "Fill with empty columns" button
            add_na_columns <- function(tdata, expdesign) {
                print("Adding NA columns...")
                
                reps <- table(expdesign[1, ])
                max_reps <- max(reps)
                tedes <- expdesign
                
                added_columns <- c()  # Track added columns for the summary
                
                # Add the new columns
                for (cond in unique(tedes[1, ])) {
                    tt <- tedes[, tedes[1, ] == cond, drop = F]
                    if (!is.na(reps[as.character(cond)])) {
                        for (i in seq_len(max_reps - reps[as.character(cond)])) {
                            # Create a new column that matches the number of rows in tedes
                            new_col <- c(cond, max(tt[2, ]) + i, 1)
                            
                            # Name the new column as new_oldname_numbering
                            oldname <- colnames(tt)[1]  # Assuming the first column is the old name
                            new_col_name <- paste0("new_", oldname, "_", i)  # Add numbering for each new column
                            
                            # Add the new column to tedes
                            tedes <- cbind(tedes, new_col)
                            
                            # Rename the new column in tedes
                            colnames(tedes)[ncol(tedes)] <- new_col_name
                            
                            # Add an NA column to tdata (the actual data)
                            tdata[new_col_name] <- NA
                            tdata[, new_col_name] <- as.numeric(tdata[, new_col_name])
                            class(tdata[, new_col_name]) <- "quant"
                            
                            # Track the added column name for the summary
                            added_columns <- c(added_columns, new_col_name)
                        }
                    }
                }
                
                # Store the added columns in the reactive value
                added_columns(added_columns)
                
                # log
                tlog <- log_operations()
                tlog[["added_columns"]] <- added_columns
                log_operations(tlog)
                
                
                # Reorder columns according to experimental design
                tedes <- tedes[, order(tedes[1, ], tedes[2, ])]
                pexp_design(tedes)
                
                # Update the processed tables
                id_column <- tdata[, grep("id", sapply(tdata, class)), drop = FALSE]
                quant_columns <- tdata[, grep("quant", sapply(tdata, class)), drop = FALSE]
                quant_columns <- quant_columns[, colnames(pexp_design())]
                
                # Update the processed table
                #print(head(quant_columns))
                processed_table(cbind(id_column, quant_columns))  # Make sure new columns are added here
                
            }
            
            ##########################################################################
            ##########################################################################
            # Remove any previously added na columns
            remove_na_columns <- function(tdata, expdesign) {
                print("Removing NA columns...")
                
                added_cols <- added_columns()
                
                if (length(added_cols) > 0) {
                    # Remove the added NA columns from tdata
                    tdata <- tdata[, !(colnames(tdata) %in% added_cols)]
                    
                    # Remove the added NA columns from expdesign
                    tedes <- expdesign[, !(colnames(expdesign) %in% added_cols)]
                    pexp_design(tedes)
                    
                    # Clear the added_columns reactive value
                    added_columns(c())
                    
                    # Update the processed table
                    id_column <- tdata[, grep("id", sapply(tdata, class)), drop = FALSE]
                    quant_columns <- tdata[, grep("quant", sapply(tdata, class)), drop = FALSE]
                    processed_table(cbind(id_column, quant_columns))
                }
            }
            
            
            
            ##########################################################################
            ##########################################################################
            # Apply maximum number of missing values per feature
            filter_nas <- function(tdata,  max_na) {
                print("Filtering NAs...")
                if (!is.null(max_na)) {
                    to_delete <- rowSums(is.na(tdata[, -1])) <= max_na
                    tdata <- tdata[to_delete, , drop=F]
                    oc <- other_cols()
                    if (!is.null(oc)) {
                        oc <- oc[to_delete, ,drop=F]
                    }
                    other_cols(oc)
                    tlog <- log_operations()
                    tlog[["max_na"]] <- max_na
                    log_operations(tlog)
                    
                }
                processed_table(tdata)
            }
            
            ##########################################################################
            ##########################################################################
            # Log transformation
            log_transformation <- function(tdata, logtrafo) {
                print("Log transformation...")
                if (!logtrafo) {
                    ttt <- as.matrix(tdata[, -1, drop=F])
                    ttt <- log2(ttt)
                    ttt[!is.finite(ttt)] <- NA
                    tdata[, -1] <- ttt
                }
                processed_table(tdata)
            }
            
            ##########################################################################
            ##########################################################################
            # Normalization
            normalize_data <- function(tdata, method) {
                print("Normalizing data...")
                # Apply normalization method
                if (method != "none") {
                    if (method == "colMedians") {
                        tdata[, -1] <- t(t(tdata[, -1]) - colMedians(as.matrix(tdata[, -1]), na.rm = TRUE))
                    } else if (method == "colMeans") {
                        tdata[, -1] <- t(t(tdata[, -1]) - colMeans(as.matrix(tdata[, -1]), na.rm = TRUE))
                    } else if (method == "cyclicloess") {
                        tdata[, -1] <- limma::normalizeBetweenArrays(as.matrix(tdata[, -1]), method = "cyclicloess")
                    }
                    tlog <- log_operations()
                    tlog[["normalization"]] <- method
                    log_operations(tlog)
                    processed_table(tdata)
                }
            }
            
            
            ##########################################################################
            ##########################################################################
            # Summarization
            summarize_data <- function(tdata, method) {
                if (sum(duplicated(tdata[,1])) > 0) {
                    print("Summarizing data...")
                    
                    
                    # update other_cols by merging values using bars
                    o_cols <- other_cols()
                    if (!is.null(o_cols)) {
                        to_rbind <- by(o_cols, tdata[, 1], function(x) {
                            apply(x, 2, function(y) paste(as.character(y), collapse="|"))
                        })
                        o_cols <- do.call(rbind, as.list(to_rbind))
                    }
                    
                    # Apply summarization method
                    if (method != "none") {
                        if (method == "colSums") {
                            tdata <- aggregate(. ~ tdata[, 1], data = tdata[, -1], FUN = sum, na.rm = TRUE)
                        } else if (method == "colMeans") {
                            tdata <- aggregate(. ~ tdata[, 1], data = tdata[, -1], FUN = mean, na.rm = TRUE)
                        } else if (method == "colMedians") {
                            tdata <- aggregate(. ~ tdata[, 1], data = tdata[, -1], FUN = median, na.rm = TRUE)
                        } else if (method == "medianPolish") {
                            tdata <- MsCoreUtils::aggregate_by_vector(as.matrix(tdata[, -1]), tdata[, 1], FUN = medianPolish, na.rm = TRUE)
                            tdata <- data.frame(ids=rownames(tdata), tdata)
                        }
                        tlog <- log_operations()
                        tlog[["summarization"]] <- method
                        log_operations(tlog)
                        enable("proceed_to_apps")
                        
                    } 
                    
                    if (!is.null(o_cols))
                        other_cols(o_cols[as.character(tdata[, 1]), ])
                    
                    # Update processed table after adjustments
                    processed_table(tdata)
                }
            }
            
            ##########################################################################
            ##########################################################################
            # Batch Effect Detection
            check_batch_effect <- function(tdata, expdesign) {
                print("Checking for batch effects...")
                # Ensure that there is enough data
                if (is.null(tdata) || ncol(tdata) < 2 || nrow(tdata) < 10) {
                    return()
                }
                
                # Store total number of features before processing
                total_features <- nrow(tdata)
                
                # Prepare the data
                texp_design <- expdesign
                
                # Select only quantitative columns and remove columns with all NAs
                tdata <- tdata[, -1, drop = FALSE]
                # Filter NA: we go full in to avoid strange effect from potential
                # imputation in the method
                tdata <- tdata[, colSums(!is.na(tdata)) > 0, drop=F]
                tdata <- tdata[complete.cases(tdata), , drop=F]
                
                texp_design <- texp_design[, colnames(tdata), drop = FALSE]
                
                # # Store number of features after removing NAs
                # features_without_na <- nrow(tdata[complete.cases(tdata), ])
                
                # # Filter out rows with missing values
                # tdata <- tdata[complete.cases(tdata), ]
                
                # Get batch labels
                batch_labels <- batch_info()
                
                # Ensure batch_labels has at least two levels
                if (length(unique(batch_labels)) < 2) {
                    return()
                }
                
                # Align `batch_labels` to match the data columns
                batch_labels <- batch_labels[match(colnames(tdata), colnames(texp_design))]
                
                # Filter rows for having at least one value per batch
                for (b in unique(batch_labels)) {
                    tdata <- tdata[rowSums(tdata[, batch_labels == b, drop=F]) > 0, , drop=F]
                }
                
                # Create a samples data frame required for BEclear
                sample_ids <- colnames(tdata)  # Assuming column names of tdata are the sample IDs
                samples <- data.frame(sample_id = sample_ids, batch_id = batch_labels)
                
                
                # Use BEclear to calculate batch effects
                # batch_effect_results <- tryCatch({
                batch_effect_results <-  BEclear::calcBatchEffects(
                    data = tdata, 
                    samples = samples,
                    adjusted = TRUE, 
                    method = "fdr"
                )
                
                if (is.null(batch_effect_results)) return() # Exit if error occurred
                
                # Extract median differences and p-values from the results
                mdifs <- batch_effect_results$med
                pvals <- batch_effect_results$pval
                summary <- calcSummary(medians = mdifs, pvalues = pvals, pvaluesTreshold = 0.01, mediansTreshold = 0.5)
                return(list(summary=summary, outdata=tdata))
            }
            
            
            #########################################################################
            #################### Batch Effect Correction ###########################
            batch_correction <- function(tdata, expdesign, batch_labels, method) {
                req(!SM$restoring())
                print("Correcting batch effects...")
                
                withProgress(message = 'Running batch correction', value = 0, {
                    # Align `batch_labels` to match the data columns
                    batch_labels <- batch_info()
                    batch_labels <- batch_labels[names(batch_labels) %in% colnames(tdata)]
                    
                    id_col <- NULL
                    # Determine the number of features with significant batch effects (p < 0.01)
                    if (method == "None") {
                        return(NULL)
                    } 
                    # Filter rows for having at least one value per batch
                    rownames(tdata) <- paste("r", 1:nrow(tdata))
                    red_for_correction <- tdata
                    # Reduced data to batch_labeled set
                    red_for_correction <- cbind(red_for_correction[, 1], red_for_correction[, names(batch_labels)])                    
                    for (b in unique(batch_labels)) {
                        red_for_correction <- red_for_correction[rowSums(!is.na(red_for_correction[, names(batch_labels[batch_labels == b]), drop=F])) > 1, , drop=F]
                    }
                    if (nrow(red_for_correction) < 10) {
                        sendSweetAlert(session,
                                       title = "Batch Effect Correction Error",
                                       text = "Not enough data after filtering for missing values to run batch effect correction.",
                                       type = "error")
                        return(NULL)
                    }
                    
                    id_col <- tdata[, 1]
                    if (method == "limma") {
                        
                        # Use limma to calculate and correct batch effects
                        batch_effect_corrected <- tryCatch({
                            removeBatchEffect(red_for_correction[, -1], batch = as.factor(batch_labels))
                        }, error = function(e) {
                            sendSweetAlert(session,
                                           title = "Batch Effect Correction Error",
                                           text = paste("An error occurred during batch effect correction with limma:", e$message),
                                           type = "error")
                            return(NULL)
                        })
                    } else if (method == "Combat") {
                        
                        # Use Combat to correct batch effects
                        batch_effect_corrected <- tryCatch({
                            sva::ComBat(dat = red_for_correction[, -1], batch = as.factor(batch_labels), mod = NULL, par.prior = TRUE, prior.plots = FALSE)
                        }, error = function(e) {
                            sendSweetAlert(session,
                                           title = "Batch Effect Correction Error",
                                           text = paste("An error occurred during batch effect correction with Combat:", e$message),
                                           type = "error")
                            return(NULL)
                        })
                    } else {
                        sendSweetAlert(session,
                                       title = "Batch Effect Correction Error",
                                       text = "Unknown batch correction method selected.",
                                       type = "error")
                        return()
                    }
                    
                    if (is.null(batch_effect_corrected)) return() # Exit if error occurred
                    
                    tlog <- log_operations()
                    tlog[["batch_correction"]] <- method
                    log_operations(tlog)
                    
                    # Increment progress to 90%
                    incProgress(0.9, detail = "Updating processed table with corrected data")
                    
                    # # Get the ID column (which is non-numeric column in processed_table)
                    # origdata <- processed_table()
                    # # Get data after batch correction
                    # origdata[, colnames(batch_effect_corrected)] <- batch_effect_corrected
                    
                    # write data back
                    tdata[rownames(batch_effect_corrected), names(batch_labels)] <- batch_effect_corrected
                    processed_table(tdata)
                    
                    # Increment progress to 100%
                    incProgress(1, detail = "Batch correction completed")
                    
                })
            }
            
            
            ## Delay reaction to selecting rows in data table
            triggerUpdate <- debounce(reactive(list( input$normalization,
                                                     input$summarize,
                                                     input$max_na,
                                                     input$add_na_columns,
                                                     removed_cols(),
                                                     input$logtrafo,
                                                     expDesign$next_tab(),
                                                     input$batch_correction_method)
                                               
            ),1000)
            
            
            
            ##########################################################################
            ##########################################################################
            # Wrapper for data manipulation and adjustment
            
            # Separate observer for normalization to avoid multiple reactivity loops
            observe({
                triggerUpdate()
                
                isolate({
                    init_data()
                    
                    tdata <- processed_table()
                    
                    if (is.null(tdata)) return() # Exit if tdata is NULL
                    remove_cols(tdata, removed_cols(), exp_design())
                    
                    if (input$add_na_columns) {
                        add_na_columns(processed_table(), pexp_design())
                    } else {
                        remove_na_columns(processed_table(), pexp_design())
                    }
                    
                    check_balance(pexp_design())
                    
                    filter_nas(processed_table(), input$max_na)
                    
                    log_transformation(processed_table(), input$logtrafo)
                    
                    normalize_data(processed_table(), input$normalization)
                    
                    summarize_data(processed_table(), input$summarize)
                    
                    # set classes of processed_table to id and quant
                    id_column <- processed_table()[,1, drop = FALSE]
                    class(id_column[,1]) <- "id"
                    quant_columns <- processed_table()[, -1]
                    for (i in 1:ncol(quant_columns)) {
                        class(quant_columns[,i]) <- "quant"
                    }
                    processed_table(cbind(id_column, quant_columns))
                    
                    batch_correction(processed_table(), pexp_design(), batch_info(), input$batch_correction_method)
                    
                    #print(colnames(processed_table()))
                    #print(pexp_design())
                    
                })
                
            })
            
            
            ##########################################################################
            ##########################################################################
            # "Summary:"
            output$ptable_summary <- renderText({
                print("Generating summary...")
                
                # Get the processed table, which reflects any changes made in data treatment
                tdata <- processed_table()
                isolate({
                    
                    # Ensure processed_table is available and valid
                    if (is.null(tdata) || ncol(tdata) < 2 || nrow(tdata) < 10) { 
                        shiny::validate(
                            need(FALSE, "Data matrix too small for generating summary")
                        )
                        return()
                    }
                    
                    # Identify "id" and "quant" columns
                    id_column <- tdata[,1]
                    quant_columns <- tdata[, -1, drop = FALSE]
                    
                    
                    # Ensure there is still enough data after filtering quantitative columns
                    if (ncol(quant_columns) < 2 || nrow(quant_columns) < 10) {
                        shiny::validate(
                            need(FALSE, "Not enough quantitative columns or rows for summary generation")
                        )
                        return()
                    }
                    
                    added_columns <- added_columns()  # Store the result in a variable
                    
                    # check for batch effect
                    tout <- check_batch_effect(processed_table(), pexp_design())
                    summary <- tout$summary
                    
                    # Determine the number of features with significant batch effects (p < 0.01)
                    significant_batches <-  ifelse(is.null(summary), 0, nrow(summary)) 
                    
                    
                    # Prepare the summary text
                    paste(
                        "<p><b>Before Filling:</b>
          <br/><b>Size:</b> The current data table contains ", 
                        ncol(quant_columns), " samples and ", nrow(quant_columns), " features.
          <br/><b>Missingness: </b>The proportion of missing values is ",
                        round(sum(is.na(quant_columns)) / length(as.matrix(quant_columns)), 3), 
                        ", with the number of missing values ranging from ", min(colSums(is.na(quant_columns))), 
                        " to ", max(colSums(is.na(quant_columns))), " per sample.
          <br/><b>Range:</b> The dynamic range (on log-scale) is from ",
                        round(min(quant_columns, na.rm = TRUE), 2), " to ", round(max(quant_columns, na.rm = TRUE), 2), 
                        ".<br/><b>Summarization:</b> The ID column contains ", 
                        ifelse(sum(duplicated(id_column)) > 0, "non-unique IDs, and thus needs summarization.", "unique IDs, so summarization is not required."), 
                        "<br/><br/><b>Batch Effect:</b> ", ifelse(
                            length(unique(batch_info())) > 1, 
                            ifelse(min(table(batch_info())) >= 2, 
                                   ifelse(
                                       significant_batches > 0, 
                                       paste("There are currently", significant_batches, " features potentially affected by batch effects (p < 0.01)."), 
                                       "No significant batch effects detected."),
                                   "There need to be at least 2 samples in each batch for running batch effect correction."),
                            "Only one batch detected, no batch effect correction needed."),
                        "<br/>Used batch effect correction method: ", input$batch_correction_method, ".",
                        "<br/><br/><b>After Filling:</b>", 
                        ifelse(length(added_columns()) == 0, " No columns added", paste(" New empty columns added: ", paste(added_columns(), collapse = ", "))),  "</p>"
                    )
                })
            })

            ##########################################################################
            ##########################################################################
            ## detect id type
            observeEvent(input$map2uniprot, {
                ids <- unique(processed_table()[, 1])
                regex_uniprot <- "^[A-NR-Z][0-9]{5}(-[0-9]+)?$|^[OPQ][0-9][A-Z0-9]{3}[0-9](-[0-9]+)?$"
                detect_key_type <- function(ids) {
                    if (all(grepl("^ENSG[0-9]+", ids))) {
                        return("Ensembl")
                    } else if (all(grepl("^[0-9]+$", ids))) {
                        return("GeneID")
                    } else if (all(grepl(regex_uniprot, ids))) {
                        # Simple UniProtKB pattern (e.g., P12345, Q8NBP7)
                        return("UniProtKB")            
                    } else {
                        return("Gene_Name")
                    }
                }
                
                keytype <- detect_key_type(ids)
                message("Detected key type: ", keytype)
                if (isTRUE(SM$restoring())) return()   # don't clobber restored inputs
                updatePickerInput(session, "idtype", 
                                  selected = keytype,
                                  choices = c("Ensembl", "GeneID", "UniProtKB", "Gene_Name"))
            })
            
            
            
            ##########################################################################
            ##########################################################################
            ## Send further to next tab
            observeEvent(input$proceed_to_apps, {
                # Ensure the processed_table is up to date
                tdata <- processed_table()  # Fetch the current processed_table data
                if (!is.null(tdata)) {
                    processed_table(tdata)  # Reassign the same data to ensure it is up-to-date
                    print("processed_table has been updated before proceeding.")
                }
                
                # add Uniprot IDs
                map_gene_to_main_uniprot <- function(ids) {
                    
                    keytype <- input$idtype
                    
                    # Create UniProt.ws object with species restriction
                    # try connecting and return NULL if it fails
                    up <- NULL
                    tryCatch({
                        up <- UniProt.ws()
                    }, error = function(e) {
                        warning("Failed to connect to UniProt.ws: ", conditionMessage(e))
                        showModal(modalDialog(
                            title = "Warning", 
                            "Failed to connect to UniProt.ws. Please check your internet connection or try again later.",
                            size = "s",
                            easyClose = T
                        ))
                        
                        return(NULL)
                    })
                    
                    available_keytypes <- keytypes(up)
                    
                    if (!(keytype %in% available_keytypes)) {
                        message("Detected key type", keytype, "is not available in UniProt.ws.\nAvailable key types include:\n",
                                paste(available_keytypes, collapse = ", "))
                        showModal(modalDialog(
                            title = "Warning", 
                            paste("Detected key type", keytype, "is not available in UniProt.ws. Please check your input."),
                            size = "s",
                            easyClose = T
                        ))
                        return(NULL)
                    }    
                    
                    
                    # Exit early if input is UniProt accessions
                    if (keytype == "UniProtKB") {
                        message("Input appears to be UniProtKB accessions. No mapping performed.")
                        showModal(modalDialog(
                            title = "Warning", 
                            "Input appears to be UniProtKB accessions. No mapping performed.",
                            size = "s",
                            easyClose = T
                        ))
                        return(NULL)
                    }
                    
                    
                    # Columns we want (filter down later)
                    cols <- c("accession")
                    cols <- intersect(cols, columns(up))
                    
                    
                    # Fetch mappings
                    # Helper function to perform chunked select() calls
                    chunked_uniprot_select <- function(up, ids, keytype, cols, chunk_size = 1000) {
                        # Split IDs into chunks
                        id_chunks <- split(ids, ceiling(seq_along(ids) / chunk_size))
                        
                        # For collecting results
                        all_results <- vector("list", length(id_chunks))
                        
                        for (i in seq_along(id_chunks)) {
                            chunk <- id_chunks[[i]]
                            if (length(chunk) == 0) next
                            message(sprintf("Querying chunk %d of %d...", i, length(id_chunks)))
                            incProgress(1 / length(id_chunks), detail = sprintf("Processing chunk %d of %d", i, length(id_chunks)))                            
                            q <- if (keytype == "Gene_Name")
                                list(taxId = 9606, ids = chunk)
                            else
                                list(ids = chunk)
                            res <- tryCatch({
                                
                                tmp <- mapUniProt(
                                    from      = keytype,
                                    to        = "UniProtKB-Swiss-Prot",
                                    query     = q,
                                    verbose = TRUE,
                                    #debug = TRUE,
                                    columns   = cols,
                                    paginate = F # Pagination is not needed here, we handle chunks manually
                                    # pageSize = chunk_size      # up to 10 000 is accepted
                                )
                                tmp[!duplicated(tmp$From),,drop=F]
                            }, error = function(e) {
                                warning(sprintf("Mapping failed for chunk %d: %s", i, conditionMessage(e)))
                                return(data.frame())
                            })
                            all_results[[i]] <- res
                        }
                        
                        # Combine all results into one data.frame
                        do.call(rbind, all_results)
                    }
                    
                    res <- tryCatch({
                        chunked_uniprot_select(up, ids = ids, keytype = keytype, cols = cols)
                    }, error = function(e) {
                        warning("Mapping failed: ", conditionMessage(e))
                        showModal(modalDialog(
                            title = "Warning", 
                            "Mapping failed. Please check your input.",
                            size = "s",
                            easyClose = T
                        ))
                        return(data.frame())
                    })
                    
                    # Deduplicate: keep one UniProt per gene (based on keytype)
                    res_unique <- res[!duplicated(res$From), ]
                    
                    # Report unmapped
                    mapped <- unique(res_unique$From)
                    unmapped <- setdiff(ids, mapped)
                    
                    list(
                        mapping = res_unique,
                        unmapped = unmapped,
                        keytype = keytype
                    )
                }
                # progressbar
                withProgress(message = 'Mapping identifiers to UniProt accession numbers via Uniprot API (this can take a while)', value = 0, {
                    
                    if(input$map2uniprot) {
                        id_column <- processed_table()[, grep("id", sapply(processed_table(), class)), drop = TRUE]
                        uniprots <- map_gene_to_main_uniprot(id_column)$mapping
                        uniprots <- as.data.frame(uniprots)
                        new_table <- other_cols()
                        if (!is.null(uniprots) && nrow(uniprots) > 0) {
                            # Safe mapping: build a named vector for lookup
                            uniprot_lookup <- setNames(uniprots$Entry, uniprots$From)
                            # Assign UniProt values using lookup by id_column
                            matched_uniprots <- uniprot_lookup[as.character(id_column)]  # names(id_column) will match names(uniprot_lookup)
                            matched_uniprots <- matched_uniprots[!is.na(matched_uniprots)]
                            new_table <- cbind(Uniprot = NA, new_table)
                            new_table[id_column %in% names(matched_uniprots), "Uniprot"] <- matched_uniprots
                        }
                        other_cols(new_table)
                        
                    }
                })
                
                
                # Proceed to the next tab
                updateTabsetPanel(parent, "mainpage", selected = "apps")
                if (!is.null(next_tab())) {
                    next_tab(paste0(next_tab(), "_new"))
                } else {
                    next_tab("ready")
                }
                
                
                print("Proceeding to interaction with apps, processed_table has been updated.")
            })
            
            
            
            
            
            ##########################################################################
            ##########################################################################
            ### PCA Plot
            
            output$pca_plot <- renderGirafe({
                print("PCA Plot")
                
                # Get the processed table, which reflects any changes made in data treatment
                tdata <- processed_table()
                
                # Check if processed_table is available and valid
                if (is.null(tdata) || ncol(tdata) < 3 || nrow(tdata) < 10) { 
                    shiny::validate(
                        need(FALSE, "Data matrix too small to perform PCA after batch correction")
                    )
                    return()
                }
                
                # Convert the processed_table to the format before batch correction
                # Identify "id" and "quant" columns
                id_column <- tdata[, 1]
                quant_columns <- tdata[,-1, drop = FALSE]
                
                
                # Ensure there is still enough data after filtering quantitative columns
                if (ncol(quant_columns) < 2) {
                    shiny::validate(
                        need(FALSE, "Not enough quantitative columns for PCA calculation")
                    )
                    return()
                }
                
                # Remove columns with all NAs
                quant_columns <- quant_columns[, colSums(!is.na(quant_columns)) > 0, drop = FALSE]
                
                # Remove constant columns to avoid PCA errors
                constant_columns <- apply(quant_columns, 2, function(col) var(col, na.rm = TRUE) == 0)
                quant_columns <- quant_columns[, !constant_columns]
                
                # Remove rows with missing or infinite values
                quant_columns <- quant_columns[complete.cases(quant_columns),]
                
                # Check if there is still enough data for PCA
                shiny::validate(
                    need(ncol(quant_columns) > 2, "Not enough columns with variance for PCA calculation"),
                    need(nrow(quant_columns) > 20, "Not enough samples for PCA calculation")
                )
                
                # Filter for non-zero sds
                row_sd <- apply(quant_columns, 1, sd, na.rm = TRUE)
                quant_columns <- quant_columns[row_sd > 0, ]
                
                # Perform PCA
                pca <- prcomp(t(quant_columns), scale = TRUE, retx = TRUE)
                loadings <- pca$x
                
                # Get the current experimental design to align replicate and batch information
                texp_design <- pexp_design()
                
                # Ensure replicate_info() and batch_info() match the data after processing
                # Find the common indices based on column names that are retained after processing
                valid_indices <- match(colnames(quant_columns), colnames(texp_design))
                replicate_info_filtered <- replicate_info()[valid_indices]
                batch_info_filtered <- batch_info()[valid_indices]
                group_info_filtered <- group_info()[valid_indices]
                
                # Check if the valid_indices has NAs which mean the alignment was not successful
                if (any(is.na(valid_indices))) {
                    shiny::validate(
                        need(FALSE, "Mismatch between data and experimental design. Ensure alignment of data columns.")
                    )
                    return()
                }
                
                # Create data frame for plotting
                pca_df <- data.frame(
                    PC1 = loadings[, 1], 
                    PC2 = loadings[, 2], 
                    Group = as.factor(group_info_filtered),
                    Replicate = as.factor(replicate_info_filtered), 
                    Batch = as.factor(batch_info_filtered),
                    Name = rownames(loadings)
                )
                
                pca_df$tooltip_info <- paste0(
                    "Sample: ", pca_df$Name, "\n",
                    "Group: ", pca_df$Group, "\n",
                    "Replicate: ", pca_df$Replicate, "\n",
                    "Batch: ", pca_df$Batch
                )
                
                # Use ggplot2 for combined PCA plot
                p <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Group, shape = Batch, label = Name, text = tooltip_info)) +
                    geom_point(size = 3) +
                    labs(title = paste0("PCA plot with first two components"),
                         x = "PC1", y = "PC2") +
                    theme_minimal() +
                    
                    geom_point_interactive(aes(color = Group, shape = Batch, tooltip = tooltip_info), size = 3) +
                    
                    theme(legend.position = "right",
                          plot.title = element_text(size = 12, face = "bold", hjust = 0.5))
                
                if (input$show_pca_labels) {
                    p <- p + geom_text(data = pca_df, aes(label = Name), , nudge_y = 3, size = 3)
                }
                output$download_pca <- downloadHandler(
                    filename = function() {
                        paste("OmicsQ_pca_plot_", Sys.Date(), ".pdf", sep = "")
                    },
                    content = function(file) {
                        pdf(file)
                        print(p)
                        dev.off()
                    }
                )
                girafe(ggobj = p)
            })
            
            
            
            
            ##########################################################################
            ##########################################################################
            ##### Correlation Plot
            output$corrplot <- renderPlot({
                print("corrplot")
                tdata <- processed_table()
                
                # Check if processed_table is available and valid
                if (is.null(tdata) || ncol(tdata) < 3 || nrow(tdata) < 10 || ncol(tdata) < 2) { 
                    shiny::validate(
                        need(FALSE, "Data matrix too small to perform correlation analysis")
                    )
                    return()
                }
                
                # Identify "quant" columns
                quant_columns <- tdata[, -1, drop = FALSE]
                
                # Remove columns with all NAs
                quant_columns <- quant_columns[, colSums(!is.na(quant_columns)) > 0, drop = FALSE]
                
                # Remove constant columns to avoid errors in correlation calculation
                constant_columns <- apply(quant_columns, 2, function(col) var(col, na.rm = TRUE) == 0)
                quant_columns <- quant_columns[, !constant_columns]
                
                # Check for missing or infinite values
                quant_columns <- quant_columns[complete.cases(quant_columns),]
                
                # Ensure there is still enough data for correlation plot
                shiny::validate(
                    need(ncol(quant_columns) > 2, "Not enough columns with variance for correlation calculation"),
                    need(nrow(quant_columns) > 20, "Not enough samples for correlation calculation")
                )
                
                # Compute correlation matrix
                correlation_matrix <- cor(quant_columns, use = "pairwise.complete.obs")
                
                # Scale the correlation range
                corr_range <- seq(-1, 1, 0.01)
                corr_col <- gplots::redblue(200)
                if (input$scale_corrplot) {
                    # Restrict red to negative values
                    if(min(correlation_matrix, na.rm = TRUE) >= 0) {
                        corr_range <- 201
                        corr_col <- gplots::colorpanel(200, "white","blue")  # Reverse the color scale
                    } 
                }
                
                # Plot the correlation matrix with customized color scale and legend position
                par(mar = c(5.1, 4.1, 5.1, 3))
                gplots::heatmap.2(correlation_matrix,
                                  main = "Pairwise sample correlations",
                                  symm = TRUE,        # Symmetrical plot
                                  scale = "none",     # No scaling
                                  col = corr_col,  # Blue to red color scale
                                  breaks = corr_range,  # Breaks from -1 to 1 for colors
                                  trace = "none",     # No trace lines
                                  cex.main = 1.5,     # Size of the main title
                                  
                                  # Add cell borders
                                  sepwidth = c(0.01, 0.01),  # Width of the separation between cells
                                  sepcolor = "white",        # Color of the separation lines (white)
                                  colsep = 1:ncol(correlation_matrix),  # Add separation for all columns
                                  rowsep = 1:nrow(correlation_matrix),  # Add separation for all rows
                                  # adapt margins to max label length
                                  margins = c(12, 12),  # Margins for the plot
                                  dendrogram = "both" # Show dendrograms on both axes
                )
                output$download_corrplot <- downloadHandler(
                    filename = function() {
                        paste("OmicsQ_correlation_plot_", Sys.Date(), ".pdf", sep = "")
                    },
                    content = function(file) {
                        pdf(file, )
                        # Plot the correlation matrix with customized color scale and legend position
                        
                        gplots::heatmap.2(correlation_matrix,
                                          main = "Pairwise sample correlations",
                                          symm = TRUE,        # Symmetrical plot
                                          scale = "none",     # No scaling
                                          col = corr_col,  # Blue to red color scale
                                          breaks = corr_range,  # Breaks from -1 to 1 for colors
                                          trace = "none",     # No trace lines
                                          cex.main = 1.5,     # Size of the main title
                                          
                                          # Add cell borders
                                          sepwidth = c(0.01, 0.01),  # Width of the separation between cells
                                          sepcolor = "white",        # Color of the separation lines (white)
                                          colsep = 1:ncol(correlation_matrix),  # Add separation for all columns
                                          rowsep = 1:nrow(correlation_matrix),  # Add separation for all rows
                                          
                                          margins = c(12, 12),  # Margins for the plot
                                          dendrogram = "both" # Show dendrograms on both axes
                        )
                        dev.off()
                    }
                )
                
                
            })
            
            
            ##########################################################################
            ##########################################################################
            ############### Plot Missing Value Distributions
            output$missingplot <- renderPlot({
                print("Missing Value Plot")
                
                # Get the processed table, which reflects any changes made in data treatment
                tdata <- processed_table()
                
                # Check if processed_table is available and valid
                if (is.null(tdata) || ncol(tdata) < 2 || nrow(tdata) < 10) { 
                    shiny::validate(
                        need(FALSE, "Data matrix too small to plot missing value distributions")
                    )
                    return()
                }
                
                
                # Get the full dataset
                # Initialize processed_table with 'id' and 'quant' columns from process_table
                full_data <- process_table()
                full_data <- full_data[, grep("quant", sapply(full_data, class)), drop = FALSE]
                
                
                # Identify "quant" columns
                quant_columns <- tdata[, -1, drop = FALSE]
                
                # Calculate the proportion of missing values per sample
                missing_dist <- missing_dist_orig <- NULL
                if (input$show_missing_rows) {
                    missing_dist <- table(rowSums(is.na(quant_columns)))
                    missing_dist_orig <- table(rowSums(is.na(full_data)))
                } else {
                    missing_dist <- colSums(is.na(quant_columns))
                    missing_dist_orig <- colSums(is.na(full_data))
                }
                
                # set colors
                if (input$show_missing_rows){ 
                    diff_len <- length(missing_dist_orig) - length(missing_dist)
                    colors <- c(rep("lightblue", length(missing_dist)), rep("darkblue", diff_len))
                    missing_dist <- missing_dist_orig
                } else {
                    colors <- c("lightgreen", "darkgreen")
                    missing_dist <- rbind(missing_dist, missing_dist_orig)
                }
                
                # plot with ifelse to distinguish rows/columns
                par(mar = c(8.1, 4.1, 4.1, 2.1))
                barplot(missing_dist, 
                        main = ifelse(input$show_missing_rows, "Missing Values per Feature", "Missing Values per Sample"),
                        xlab = ifelse(input$show_missing_rows, "Number of missing values", ""),
                        ylab = "Counts",
                        # leave space for legend
                        ylim = c(0, max(missing_dist) * 1.2),
                        col = colors,
                        border = "black",
                        cex.names = 0.7,  # Adjust label size
                        beside = TRUE,
                        # rotate labels
                        las = 2
                )
                if (input$show_missing_rows)
                    legend("topright", legend = c("Filtered", "Original"), fill = c("lightblue","darkblue"), bty = "n")
                else {
                    legend("topright", legend = c("Filtered", "Original"), fill = c("lightgreen", "darkgreen"), bty = "n")
                }
                output$download_missingplot <- downloadHandler(
                    filename = function() {
                        paste("OmicsQ_missing_value_plot_", Sys.Date(), ".pdf", sep = "")
                    },
                    content = function(file) {
                        pdf(file)
                        par(mar = c(8.1, 4.1, 4.1, 2.1))
                        barplot(rbind(missing_dist, missing_dist_orig), 
                                main = ifelse(input$show_missing_rows, "Missing Values per Feature", "Missing Values per Sample"),
                                xlab = ifelse(input$show_missing_rows, "Number of missing values", ""),
                                ylab = "Counts",
                                col = ifelse(input$show_missing_rows, c("lightblue","darkblue"), c("lightgreen", "darkgreen")),
                                border = "black",
                                cex.names = 0.7,  # Adjust label size
                                beside = TRUE,
                                legend = c("Filtered", "Original"),
                                # rotate labels
                                las = 2
                        )
                        dev.off()
                    }
                )
            })
            
            
            ##########################################################################
            ##########################################################################
            ############### Help messages
            observeEvent(input$h_balancing, sendSweetAlert(session,
                                                           title = "Create a balanced design",
                                                           text = HTML("<p align='justify'>The different experimental conditions
    (sample types) need to be at least nearly balanced. For the
              further analysis, this means that the number of replicates
              per sample should be the same for each of them. You can balance the data
              by adding empty columns and/or removing excess replicates.
              Replicates with the same number have already been summarized.<br/>
              <i>Select below</i> the samples you want to exclude.<br/>
              <i>Fill with empty columns: </i>Use this switch to reach
              the same number of replicates per condition. This does
              <b>not make a data set more balanced</b>.
              We recommend removing samples from groups that have a
              much larger number of replicates.</p>"),
                                                           type = "info", html = T
            ))
            
            observeEvent(input$h_logtrafo, sendSweetAlert(session,
                                                          title = "Log transformation",
                                                          text = HTML("<p align='justify'>In general, quantitative omics data is log-transformed.
              OmicsQ estimates whether the data was transformed from the values (data range and
              negative values). Deselecting this box will transform the data by taking the logarithm
              with base 2.</p>"),
                                                          type = "info", html = T
            ))
            
            observeEvent(input$h_max_na, sendSweetAlert(session,
                                                        title = "Deleting features with low coverage",
                                                        text = HTML("<p align='justify'>Despite of the capability of PolySTest and VSClust to include features
              with missing values, we recommend removing features that have very low coverage as their measurements
              are usually very noisy. Filter for the maximum number of missing values here. The default is
              to take all features.</p>"),
                                                        type = "info", html = T
            ))
            
            observeEvent(input$h_normalization, sendSweetAlert(session,
                                                               title = "Normalization",
                                                               text = HTML("<p align='justify'>The total amount of molecules per sample can vary when loading them
              into the instrument. This systematic error can be corrected by assuming that most
              of the features do not change between samples. We offer the main normalization methods
              used in the field.<br/>
              <i>None:</i> No normalization (e.g. when the data has already been normalized)<br/>
              <i>Median: </i>Subtract the median of each sample from the log-transformed values (recommended
              as more outlier insensitive)<br/>
              <i>Mean: </i>Subtract the mean of the sample from the log-transformed values<br/>
              <i>Cyclic Loess: </i>Apply a non-linear transformation to accommodate for non-linear effects and
              large changes and batch effects in the data (see <a href='http://web.mit.edu/~r/current/arch/i386_linux26/lib/R/library/limma/html/normalizeCyclicLoess.html'>cyclicLoess</a>. <br/></p>"),
                                                               type = "info", html = T
            ))
            
            observeEvent(input$h_summarize, sendSweetAlert(session,
                                                           title = "Summarization of feature values",
                                                           text = HTML("<p align='justify'><i>This option is only needed for id columns with duplicated values.</i>
              This can be multiple features measured over all samples like peptides of the same protein.
              Check the data summary on the right side whether there are any duplicates.
              The here given options summarize the rows with the same feature name using one of the following
              methods:<br/>
              <i>None: </i>Do not summarize. In the case of duplicated id values, you won't be able to proceed (default)<br/>
              <i>Sum: </i>Here we take the sum of the non-log-transformed values. This can have impact on the normalization.<br/>
              <i>Mean: </i>Take the mean of log-transformed values<br/>
              <i>Median: </i>Take the median of log-transformed values<br/>
              <i>Robust median: </i>Take the median of log-transformed values but remove outliers
              (see <a href='https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/medpolish'>medpolish</a>)<br/>
              <i>Note: </i>Summarization happens after normalization.</p>"),
                                                           type = "info", html = T
            ))
            
            observeEvent(input$h_batch_effect, sendSweetAlert(session,
                                                              title = "Batch Effect Detection and Correction",
                                                              text = HTML("<p align='justify'>
    Batch effects occur when systematic technical variations (e.g., due to different experimental conditions, days, or instrument runs) 
    introduce biases into your data, making it difficult to discern true biological differences. 
    Detecting and correcting for these effects is essential for reliable downstream analysis.<br/>
    <b>Important:</b>Apply batch correction only when you have a good reason to assume that there are batch effects!<br/><br/>
    
    <b>Detection:</b> The <i>BEclear</i> package is used to detect batch effects by calculating the median differences and 
    p-values for each feature across different batches. Features with significant differences between batches (p < 0.01 and median difference > 0.5) 
    are flagged as being affected by batch effects.<br/><br/>
    
    <b>Correction:</b> You can select one of the following batch correction methods:<br/>
    <i>None: </i>No batch correction is applied.<br/>
    <i>limma: </i>Batch effects are corrected using the linear modeling approach from the limma package 
    (<a href='https://bioconductor.org/packages/release/bioc/html/limma.html'>limma</a>), which estimates and removes the batch effect.<br/>
    <i>ComBat: </i>This method uses an empirical Bayes framework to adjust for batch effects by shrinking the batch means and variances 
    towards the overall mean (<a href='https://bioconductor.org/packages/release/bioc/html/sva.html'>sva: ComBat</a>).<br/>
    <b>Note:</b> Batch correction is only applied when there are at least two samples per batch.<br/><br/>

    Both methods are widely used for batch correction in omics data and help reduce false positives 
    and ensure that biological signal is retained. Select the method that best fits your experimental setup.
    </p>"),
                                                              type = "info", html = T
            ))
            
            observeEvent(input$h_pcaplot, sendSweetAlert(session,
                                                         title = "PCA Plot",
                                                         text = HTML("<p align='justify'>
                                                              The PCA plot visualizes the variance in the data by projecting it onto the first two principal components.
              This helps to identify patterns, clusters, and potential outliers in the data. The points are colored by group and shaped by batch.
              You can also choose to label the points with feature names for better identification of samples with different behavior.</p>"),
                                                         type = "info", html = T
            ))
            
            observeEvent(input$h_corrplot, sendSweetAlert(session,
                                                          title = "Correlation Plot",
                                                          text = HTML("<p align='justify'>The correlation plot shows the pairwise correlations between samples.
              This helps to identify relationships and potential batch effects in the data. The color scale indicates the strength of the correlation,
              with blue representing positive correlations and red representing negative correlations. You can scale the plot to better see smaller differences 
                                                                      when the correlations are close to one.</p>"),
                                                          type = "info", html = T
            ))
            
            observeEvent(input$h_missingplot, sendSweetAlert(session,
                                                             title = "Missing Value Distribution",
                                                             text = HTML("<p align='justify'>The missing value distribution plot shows the number of missing values per feature or sample.
              This helps to identify features or samples with high missingness, which can affect the analysis. When choosing the distribution per feature, the plot counts the number of 
                                                                        features with a given number of missing values.</p>"),
                                                             type = "info", html = T
            ))
            
            observeEvent(input$h_mapping, sendSweetAlert(session,
                                                         title = "UniProt Mapping",
                                                         text = HTML("<p align='justify'>This step maps the identifiers in your data to UniProt accessions.
              This is useful for standardizing identifiers and ensuring compatibility with protein complex analysis in ComplexBrowser.
              The mapping is performed using the UniProt.ws package, which connects to the UniProt database.
              You can select the key type of your identifiers (e.g., Ensembl, GeneID, UniProtKB, Gene_Name) to ensure accurate mapping. The app
                                                                      tries to estimate the type automatically. <br/>Note that gene names are too generic 
                                                                      to be mapped to a particular species, and thus are mapped to homo sapiens.</p>"),
                                                         type = "info", html = T
            ))
            
            
            
            return(list(
                next_tab = next_tab,
                processed_table = processed_table,
                processed_table = processed_table,
                other_cols = other_cols,
                pexp_design = pexp_design
            ))
            
        }
    )
}