server.R

server <- function(input, output, session) { 
  #use_python("/opt/conda39/envs/pyscenic/bin/python") #environment containing phate
  source("global.R", local=TRUE)
  source("addPeak2GeneLinks_shiny.R")
  
  options(shiny.maxRequestSize=100*1024^3) # 2GB limit in the online version
  session$sendCustomMessage("handler_disableTabs", "sidebarMenu") # disable all tab panels (except Data Input) until files are uploaded
  
  session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " GENE REGULATORY NETWORK\nANALYSIS", " UTILITY OPTIONS")) #enable GRN tab only for RNA
  hideTab(inputId="grnTabPanel", target="scATAC-seq", session = session)
  
  hideAllLoaders() # helper function (in global.R) that initially hides all loaders
  metaD <- reactiveValues(my_project_name="-", all_lin="0")
  
  #------------------Upload tab--------------------------------
  observeEvent(input$uploadCountMatrixConfirm, {
    if(!is.null(proj_default) | !is.null(seurat_object))
    {
      showModal(modal_confirm)
    }
    else
    {
      session$sendCustomMessage("handler_disableTabs", "sidebarMenu") 
      session$sendCustomMessage("handler_disableAllButtons", "uploadCountMatrixConfirm") 
      tryCatch({
        showModal(modalDialog(div('Data upload in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        # Create the user directory for the input and output of the analysis
        metaD$my_project_name <- input$uploadCountMatrixprojectID
        minimum_cells <<- input$uploadCountMatrixminCells
        minimum_features <<- input$uploadCountMatrixminFeatures
        organism <<- input$uploadCountMatrixRadioSpecies
        userId <- session$token
        user_dir <<- paste0("./usr_temp/", userId, metaD$my_project_name, gsub(pattern = "[ ]|[:]", replacement = "_", x = paste0("_", Sys.time()))) 
        dir.create(user_dir)
        file.copy(from = input$countMatrix$datapath, to = paste0(user_dir, "/countMatrix.txt"), overwrite = TRUE)
        
        testMatrix <- read.table(paste0(user_dir, "/countMatrix.txt"))
        seurat_object <<- CreateSeuratObject(counts = testMatrix, project = metaD$my_project_name, min.cells = as.numeric(minimum_cells), min.features = as.numeric(minimum_features))
        
        init_seurat_object <<- CreateSeuratObject(counts = testMatrix, project = metaD$my_project_name, min.cells = as.numeric(minimum_cells), min.features = as.numeric(minimum_features))
        
        if(organism == "mouse")
        {
          seurat_object[["percent.mt"]] <<- PercentageFeatureSet(seurat_object, pattern = "^mt-")
          init_seurat_object[["percent.mt"]] <<- PercentageFeatureSet(init_seurat_object, pattern = "^mt-")
        }
        else
        {
          seurat_object[["percent.mt"]] <<- PercentageFeatureSet(seurat_object, pattern = "^MT-")
          init_seurat_object[["percent.mt"]] <<- PercentageFeatureSet(init_seurat_object, pattern = "^MT-")
        }
        
        updateMetadata()
        
        updateSelInpColor()
        updateGeneSearchFP()
        updateQC_choices()
        updateMaxHVGs()
        updatePCs_selection(10)
        updateUtilitiesAssays()
        disableTabsATAC()
        #update signature text area
        if(organism == "human")
          updateTextAreaInput(session, inputId="findMarkersSignatureMembers", placeholder = "PRG4\nTSPAN15\nCOL22A1\nHTRA4")
        
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " QUALITY CONTROL", " DATA NORMALIZATION\n& SCALING"))
      }, error = function(e) {
        print(paste("Error :  ", e))
        session$sendCustomMessage("handler_alert", "Data Input error. Please, refer to the help pages for input format.")
      }, finally = { # with or without error
        removeModal()
        Sys.sleep(1) # giving some time for rendering, for smoother transition
        session$sendCustomMessage("handler_enableAllButtons", "uploadCountMatrixConfirm")
      })
    }
  })
  
  observeEvent(input$upload10xExampleRNACountMatrixConfirm, { 
    if(!is.null(proj_default) | !is.null(seurat_object))
    {
      showModal(modal_confirm)
    }
    else
    {
      session$sendCustomMessage("handler_disableTabs", "sidebarMenu")
      session$sendCustomMessage("handler_disableAllButtons", "upload10xExampleRNACountMatrixConfirm") 
      tryCatch({
        showModal(modalDialog(div('Data upload in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        # Create the user directory for the input and output of the analysis
        metaD$my_project_name <- input$uploadCountMatrixprojectID
        minimum_cells <<- input$uploadCountMatrixminCells
        minimum_features <<- input$uploadCountMatrixminFeatures
        organism <<- "human"
        userId <- session$token
        user_dir <<- paste0("./usr_temp/", userId, metaD$my_project_name, gsub(pattern = "[ ]|[:]", replacement = "_", x = paste0("_", Sys.time())))
        dir.create(user_dir)
        
        testMatrix <- read.table("exampleRNA_matrix/exampleMatrix.txt")
        seurat_object <<- CreateSeuratObject(counts = testMatrix, project = metaD$my_project_name, min.cells = as.numeric(minimum_cells), min.features = as.numeric(minimum_features))
        
        init_seurat_object <<- CreateSeuratObject(counts = testMatrix, project = metaD$my_project_name, min.cells = as.numeric(minimum_cells), min.features = as.numeric(minimum_features))
        
        seurat_object[["percent.mt"]] <<- PercentageFeatureSet(seurat_object, pattern = "^MT-")
        init_seurat_object[["percent.mt"]] <<- PercentageFeatureSet(init_seurat_object, pattern = "^MT-")
        
        updateMetadata()
        
        updateSelInpColor()
        updateGeneSearchFP()
        updateQC_choices()
        updateMaxHVGs()
        updatePCs_selection(10)
        updateUtilitiesAssays()
        disableTabsATAC()
        #update signature text area
        if(organism == "human")
          updateTextAreaInput(session, inputId="findMarkersSignatureMembers", placeholder = "PRG4\nTSPAN15\nCOL22A1\nHTRA4")
        
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " QUALITY CONTROL", " DATA NORMALIZATION\n& SCALING"))
      }, error = function(e) {
        print(paste("Error :  ", e))
        session$sendCustomMessage("handler_alert", "Data Input error. Please, refer to the help pages for input format.")
      }, finally = { # with or without error
        removeModal()
        Sys.sleep(1) # giving some time for rendering for smoother transition
        session$sendCustomMessage("handler_enableAllButtons", "upload10xExampleRNACountMatrixConfirm")
      })
    }
  })
  
  observeEvent(input$upload10xRNAConfirm, {
    if(!is.null(proj_default) | !is.null(seurat_object))
    {
      showModal(modal_confirm)
    }
    else
    {
      session$sendCustomMessage("handler_disableTabs", "sidebarMenu") 
      session$sendCustomMessage("handler_disableAllButtons", "upload10xRNAConfirm")
      tryCatch({
        showModal(modalDialog(div('Data upload in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        metaD$my_project_name <- input$upload10xRNAprojectID
        minimum_cells <<- input$upload10xRNAminCells
        minimum_features <<- input$upload10xRNAminFeatures
        organism <<- input$upload10xRNARadioSpecies
        
        userId <- session$token
        user_dir <<- paste0("./usr_temp/", userId, metaD$my_project_name, gsub(pattern = "[ ]|[:]", replacement = "_", x = paste0("_", Sys.time()))) 
        dir.create(user_dir)
        file.copy(from = input$matrix$datapath, to = paste0(user_dir, "/matrix.mtx.gz"), overwrite = TRUE) #,input$matrix$name
        file.copy(from = input$barcodes$datapath, to = paste0(user_dir, "/barcodes.tsv.gz"), overwrite = TRUE) #, input$barcodes$name
        file.copy(from = input$genes$datapath, to = paste0(user_dir, "/features.tsv.gz"), overwrite = TRUE) #, input$genes$name
        
        seurat_data <- Read10X(user_dir)
        seurat_object <<- CreateSeuratObject(counts = seurat_data, project = metaD$my_project_name, min.cells = as.numeric(minimum_cells), min.features = as.numeric(minimum_features))
        
        init_seurat_object <<- CreateSeuratObject(counts = seurat_data, project = metaD$my_project_name, min.cells = as.numeric(minimum_cells), min.features = as.numeric(minimum_features))
        
        if(organism == "mouse")
        {
          seurat_object[["percent.mt"]] <<- PercentageFeatureSet(seurat_object, pattern = "^mt-")
          init_seurat_object[["percent.mt"]] <<- PercentageFeatureSet(init_seurat_object, pattern = "^mt-")
        }
        else
        {
          seurat_object[["percent.mt"]] <<- PercentageFeatureSet(seurat_object, pattern = "^MT-")
          init_seurat_object[["percent.mt"]] <<- PercentageFeatureSet(init_seurat_object, pattern = "^MT-")
        }
        
        updateMetadata()
        
        updateSelInpColor()
        updateGeneSearchFP()
        updateQC_choices()
        updateMaxHVGs()
        updatePCs_selection(10)
        updateUtilitiesAssays()
        disableTabsATAC()
        #update signature text area
        if(organism == "human")
          updateTextAreaInput(session, inputId="findMarkersSignatureMembers", placeholder = "PRG4\nTSPAN15\nCOL22A1\nHTRA4")
        
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " QUALITY CONTROL", " DATA NORMALIZATION\n& SCALING"))
      }, error = function(e) {
        print(paste("Error :  ", e))
        session$sendCustomMessage("handler_alert", "Data Input error. Please, refer to the help pages for input format.")
      }, finally = { # with or without error
        removeModal()
        Sys.sleep(1) # giving some time for rendering for smoother transition
        session$sendCustomMessage("handler_enableAllButtons", "upload10xRNAConfirm")
      })
    }
  })
  
  observeEvent(input$upload10xExampleRNA10xFilesConfirm, {
    if(!is.null(proj_default) | !is.null(seurat_object))
    {
      showModal(modal_confirm)
    }
    else
    {
      session$sendCustomMessage("handler_disableTabs", "sidebarMenu") # disable all tab panels (except Data Input) until files are uploaded
      session$sendCustomMessage("handler_disableAllButtons", "upload10xExampleRNA10xFilesConfirm")
      tryCatch({
        showModal(modalDialog(div('Data upload in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        metaD$my_project_name <- "PBMC3k"
        minimum_cells <<- 3
        minimum_features <<- 200
        organism <<- "human"
        
        
        userId <- session$token
        user_dir <<- paste0("./usr_temp/", userId, metaD$my_project_name, gsub(pattern = "[ ]|[:]", replacement = "_", x = paste0("_", Sys.time()))) 
        dir.create(user_dir)
        
        seurat_data <- Read10X("exampleRNA_10xFiles/")
        seurat_object <<- CreateSeuratObject(counts = seurat_data, project = metaD$my_project_name, min.cells = as.numeric(minimum_cells), min.features = as.numeric(minimum_features))
        
        init_seurat_object <<- CreateSeuratObject(counts = seurat_data, project = metaD$my_project_name, min.cells = as.numeric(minimum_cells), min.features = as.numeric(minimum_features))
        
        seurat_object[["percent.mt"]] <<- PercentageFeatureSet(seurat_object, pattern = "^MT-")
        init_seurat_object[["percent.mt"]] <<- PercentageFeatureSet(init_seurat_object, pattern = "^MT-") #as.double(sprintf("%.3f", as.double(seurat_object@meta.data$percent.mt)))
        seurat_object@meta.data$percent.mt <<- as.double(seurat_object@meta.data$percent.mt)
        init_seurat_object@meta.data$percent.mt <<- as.double(init_seurat_object@meta.data$percent.mt)
        
        updateMetadata()
        
        updateSelInpColor()
        updateGeneSearchFP()
        updateQC_choices()
        updateMaxHVGs()
        updatePCs_selection(10)
        updateUtilitiesAssays()
        disableTabsATAC()
        #update signature text area
        if(organism == "human")
          updateTextAreaInput(session, inputId="findMarkersSignatureMembers", placeholder = "PRG4\nTSPAN15\nCOL22A1\nHTRA4")
        
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " QUALITY CONTROL", " DATA NORMALIZATION\n& SCALING"))
      }, error = function(e) {
        print(paste("Error :  ", e))
        session$sendCustomMessage("handler_alert", "Data Input error. Please, refer to the help pages for input format.")
      }, finally = { # with or without error
        removeModal()
        Sys.sleep(1) # giving some time for rendering for smoother transition
        session$sendCustomMessage("handler_enableAllButtons", "upload10xExampleRNA10xFilesConfirm")
      })
    }
  })
  
  observeEvent(input$uploadSeuratRdsConfirm, {
    if(!is.null(proj_default) | !is.null(seurat_object))
    {
      showModal(modal_confirm)
    }
    else
    {
      session$sendCustomMessage("handler_disableAllButtons", "uploadSeuratRdsConfirm")
      
      tryCatch({
        showModal(modalDialog(div('Data upload in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        #seurat_object <<- readRDS(input$uploadRdsFile$datapath)
        #init_seurat_object <<- readRDS(input$uploadRdsFile$datapath)
        userId <- session$token
        user_dir <<- paste0("./usr_temp/", userId, gsub(pattern = "[ ]|[:]", replacement = "_", x = paste0("_", Sys.time()))) 
        dir.create(user_dir)
        file.copy(from = input$uploadRdsFile$datapath, to = paste0(user_dir, "/init_seurat_object.RDS"), overwrite = TRUE)
        
        seurat_object <<- readRDS(paste0(user_dir, "/init_seurat_object.RDS"))
        init_seurat_object <<- readRDS(paste0(user_dir, "/init_seurat_object.RDS"))
        organism <<- input$uploadRdsRadioSpecies
        
        print(seurat_object)
        print(init_seurat_object)
        
        #checking slots and enabling tabs
        if(is.null(seurat_object@meta.data))
          print("Metadata table is missing")
        else
        {
          disableTabsATAC()
          updateMetadata()
          updateSelInpColor()
          updateGeneSearchFP()
          updateQC_choices()
          updateMaxHVGs()
          updatePCs_selection(10)
          updateUtilitiesAssays()
          updateRegressOut()
          #CHECK percent.mt
          if(is.null(seurat_object@meta.data$percent.mt))
          {
            #calculate percent.mt and enable 
            if(organism == "mouse")
            {
              seurat_object[["percent.mt"]] <<- PercentageFeatureSet(seurat_object, pattern = "^mt-")
              init_seurat_object[["percent.mt"]] <<- PercentageFeatureSet(init_seurat_object, pattern = "^mt-")
            }
            else
            {
              seurat_object[["percent.mt"]] <<- PercentageFeatureSet(seurat_object, pattern = "^MT-")
              init_seurat_object[["percent.mt"]] <<- PercentageFeatureSet(init_seurat_object, pattern = "^MT-")
            }
          }
          
          session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " UTILITY OPTIONS", " QUALITY CONTROL", " GENE REGULATORY NETWORK\nANALYSIS"))
          
          #CHECK @data, VariableFeatures, @scale.data
          if( !is.null(seurat_object@assays$RNA@data) & length(VariableFeatures(seurat_object)) != 0 & !is.null(seurat_object@assays$RNA@scale.data) )
          {
            session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " DATA NORMALIZATION\n& SCALING"))
            shinyjs::show("hvgScatter_loader")
            #plot highly variable genes
            #--Rendering
            output$hvgScatter <- renderPlotly({ # tooltip
                plot1 <- VariableFeaturePlot(seurat_object)
                varplot <- plot1$data
                varplot$gene <- rownames(varplot)
                varplot$colors[varplot$colors == "yes"] <- paste0("Highly Variable genes(", length(VariableFeatures(seurat_object)), ")")
                varplot$colors[varplot$colors == "no"] <- paste0("Not Variable genes(", length(rownames(seurat_object)) - length(VariableFeatures(seurat_object)), ")")
                print("before method")
                print(seurat_object@commands$FindVariableFeatures.RNA$selection.method)
                if(seurat_object@commands$FindVariableFeatures.RNA$selection.method == "vst"){
                  p <- ggplot(varplot, aes(x=log10(mean), y=variance.standardized, color=colors, text = paste0("log10(mean expression): ", log10(mean),
                                                                                                               "\nstandardized variance: ", variance.standardized,
                                                                                                               "\ngene: ", gene)
                  )) +
                    geom_point() +
                    theme_bw() +
                    scale_color_manual(
                      values = c("red", "black")
                    )+
                    labs(x="Average Expression", y="Standardized Variance", color="")
                } else {
                  p <- ggplot(varplot, aes(x=mean, y=dispersion, color=colors, text = paste0("mean expression: ", mean,
                                                                                             "\ndispersion: ", dispersion,
                                                                                             "\ngene: ", gene)
                  )) +
                    geom_point() +
                    theme_bw() +
                    scale_color_manual(
                      values = c("red", "black")
                    )+ 
                    labs(x="Average Expression", y="Dispersion", color="")
                }
                gp <- plotly::ggplotly(p, tooltip = "text")
                gp
            }) 
            
            output$hvgTop10Stats <- renderPrint(
              {
                top10hvgs <- head(VariableFeatures(seurat_object), 10)
                cat("\nTop10 varibale features:", c(top10hvgs[1:10]))
              })
            #--End Rendering
          }
          
          #CHECK @reductions$pca and other dimensionality reductions
          if( !is.null(seurat_object@reductions$pca) )
          {
            #print 2 plots
            session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " PCA/LSI", " ADDITIONAL DIMENSIONALITY\nREDUCTION METHODS", " FEATURE INSPECTION"))
            shinyjs::show("elbowPlotPCA_loader")
            shinyjs::show("PCAscatter_loader")
            #--Rendering
            output$elbowPlotPCA <- renderPlotly(
              {
                plot1 <- ElbowPlot(seurat_object, ndims = ncol(seurat_object@reductions$pca))
                plot1_data <- plot1$data
                colnames(plot1_data)[1] <- "PC"
                colnames(plot1_data)[2] <- "SD"
                
                p <- ggplot(plot1_data, aes(x=PC, y=SD))+
                    geom_point() +
                    theme_bw() +
                    labs(x="PC", y="Standard Deviation")
                
                
                gp <- plotly::ggplotly(p, tooltip = c("x", "y"))
                print(gp)
              }
            )
            
            output$PCAscatter <- renderPlotly(
            {
                #prepare metadata
                meta <- seurat_object@meta.data
                meta$Cell_id <- rownames(meta)
                meta <- meta[, ]#meta[, c('Cell_id', 'seurat_clusters', 'orig.ident')]
                reduc_data <- data.frame()
                
                #prepare colors
                cols = colorRampPalette(brewer.pal(12, "Paired"))(1)
                
                #umap data frame
                seurat_object_reduc <- as.data.frame(seurat_object@reductions$pca@cell.embeddings)
                seurat_object_reduc <- seurat_object_reduc[, c(1:ncol(seurat_object_reduc))]
                seurat_object_reduc$Cell_id <- rownames(seurat_object_reduc)
                reduc_data <- left_join(seurat_object_reduc, meta)
                
                p <- ggplot(data=reduc_data, aes_string(x="PC_1", y="PC_2", fill="orig.ident")) +
                  geom_point(size=2, shape=19, stroke=0)+
                  scale_fill_manual(values = cols)+
                  scale_size()+
                  theme_bw() +
                  theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                        axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                        axis.title.y = element_text(face = "bold", color = "black", size = 25),
                        axis.title.x = element_text(face = "bold", color = "black", size = 25),
                        legend.text = element_text(face = "bold", color = "black", size = 9),
                        legend.title = element_text(face = "bold", color = "black", size = 9),
                        legend.position="right",
                        title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
                  labs(x="PC 1", y="PC 2", title = "", fill="Color")
                print(plotly::ggplotly(p))
              }
            )
            
            #--End rendering
            #CHECK $seurat@reductions
            updateUmapTypeChoices(names(seurat_object@reductions))
          }
          
          #CHECK $seurat_clusters
          if( !is.null(seurat_object@meta.data$seurat_clusters) )
          {
            updateInputLRclusters()
            session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " CLUSTERING", " MARKERS' IDENTIFICATION", " CELL CYCLE PHASE ANALYSIS", " DOUBLETS' DETECTION",
                                                              " LIGAND - RECEPTOR\nANALYSIS"))
            #--print cluster table
            updateClusterTab()
            #--
            
            if(!is.null(seurat_object@reductions$umap))
              session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " TRAJECTORY ANALYSIS"))
          }
          
          #CHECK Seurat@misc$markers
          if( !is.null(seurat_object@misc$markers) )
          {
            seurat_object@misc$markers <- seurat_object@misc$markers[, c(7, 6, 1:5)]
            output$findMarkersTable <- renderDataTable(seurat_object@misc$markers, options = list(pageLength = 20), filter = 'top', rownames = FALSE)
            markers_logFCBase <<- colnames(seurat_object@misc$markers)[grep(pattern = "avg_log", x = colnames(seurat_object@misc$markers))]
            export_markerGenes_RNA <<- seurat_object@misc$markers
            session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " FUNCTIONAL/MOTIF\nENRICHMENT ANALYSIS"))
          }
        }
      }, error = function(e) {
          print(paste("Error :  ", e))
          session$sendCustomMessage("handler_alert", "Data Input error. Please, refer to the help pages for input format.")
      }, finally = { # with or without error
          removeModal()
          Sys.sleep(1) # giving some time for rendering for smoother transition
          session$sendCustomMessage("handler_enableAllButtons", "uploadSeuratRdsConfirm")
        }
      )
    }
  })
  
  observeEvent(input$upload10xATACConfirm, {
    if(!is.null(proj_default) | !is.null(seurat_object))
    {
      showModal(modal_confirm)
    }
    else
    {
      session$sendCustomMessage("handler_disableTabs", "sidebarMenu") # disable all tab panels (except Data Input) until files are uploaded
      session$sendCustomMessage("handler_disableAllButtons", "upload10xATACConfirm")
      tryCatch({
        showModal(modalDialog(div('Data upload in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        
        #user dir creation
        projectNameATAC <<- input$uploadATACprojectID
        userId <- session$token
        user_dir <<- paste0("./usr_temp/", userId, projectNameATAC, gsub(pattern = "[ ]|[:]", replacement = "_", x = paste0("_", Sys.time()))) 
        dir.create(user_dir)
        file.copy(from = input$uploadATACArrow$datapath, to = paste0(user_dir, "/arrowFile.arrow"), overwrite = TRUE)
        dir.create(paste0(user_dir, "/default"))
        #select genome version and organism
        addArchRGenome(input$upload10xATACRadioSpecies)
        if(length(grep("mm", input$upload10xATACRadioSpecies)) == 1)
        {
          organism <<- "mouse"
        }
        else
        {
          organism <<- "human"
        }
        
        addArchRThreads(threads = as.numeric(input$upload10xATACThreads)) 
        
        ########################################
        ######### Create Arch Project ##########
        ########################################
        proj_default <<- ArchRProject(
          ArrowFiles = paste0(user_dir, "/arrowFile.arrow"),
          outputDirectory = paste0(user_dir, "/default"),
          copyArrows = FALSE #This is recommended so that if you modify the Arrow files you have an original copy for later usage.
        )
        
        saveArchRProject(proj_default)
        print("Project saved")
        
        updateMetadataATAC()
        
        addArchRThreads(threads = as.numeric(as.numeric(input$upload10xATACThreads)))
        updateSelectizeInput(session, "visualizeTracksGene", choices = unique(proj_default@geneAnnotation$genes$symbol), server = T)
        
        disableTabsRNA()
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " QUALITY CONTROL"))
        }, error = function(e) {
        print(paste("Error :  ", e))
        session$sendCustomMessage("handler_alert", "Data Input error. Please, refer to the help pages for input format.")
      }, finally = { # with or without error
        removeModal()
        Sys.sleep(1) # giving some time for rendering for smoother transition
        session$sendCustomMessage("handler_enableAllButtons", "upload10xATACConfirm")
      })
    }
  })
  
  observeEvent(input$upload10xExampleATACConfirm, {
    if(!is.null(proj_default) | !is.null(seurat_object))
    {
      showModal(modal_confirm)
    }
    else
    {
      session$sendCustomMessage("handler_disableTabs", "sidebarMenu") # disable all tab panels (except Data Input) until files are uploaded
      session$sendCustomMessage("handler_disableAllButtons", "upload10xExampleATACConfirm")
      tryCatch({
        showModal(modalDialog(div('Data upload in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        
        #user dir creation
        projectNameATAC <<- "example"
        userId <- session$token
        user_dir <<- paste0("./usr_temp/", userId, projectNameATAC, gsub(pattern = "[ ]|[:]", replacement = "_", x = paste0("_", Sys.time()))) 
        dir.create(user_dir)
        dir.create(paste0(user_dir, "/default"))
        
        #select genome version and organism
        addArchRGenome("hg19")
        
        
        #set number of threads
        addArchRThreads(threads = as.numeric(input$upload10xATACThreads)) 
        
        ########################################
        ######### Create Arch Project ##########
        ########################################
        proj_default <<- ArchRProject(
         ArrowFiles = "exampleATAC/arrowFile.arrow",
         outputDirectory = paste0(user_dir, "/default"),
         copyArrows = FALSE #This is recommended so that if you modify the Arrow files you have an original copy for later usage.
        )
        
        saveArchRProject(proj_default)
        
        print("Project saved")
        
        updateMetadataATAC()
        
        addArchRThreads(threads = as.numeric(input$upload10xATACThreads))
        updateSelectizeInput(session, "visualizeTracksGene", choices = unique(proj_default@geneAnnotation$genes$symbol), server = T)
        
        disableTabsRNA()
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " QUALITY CONTROL"))
      }, error = function(e) {
        print(paste("Error :  ", e))
        session$sendCustomMessage("handler_alert", "Data Input error. Please, refer to the help pages for input format.")
      }, finally = { # with or without error
        removeModal()
        Sys.sleep(1) # giving some time for rendering for smoother transition
        session$sendCustomMessage("handler_enableAllButtons", "upload10xExampleATACConfirm")
      })
    }
  })
  
  output$uploadMetadataExport <- downloadHandler(
    filename = function() { 
      paste("metadataTableATAC-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_metadata_ATAC, file, sep = "\t", quote = F, row.names = F)
    })
  
  output$uploadMetadataExportRNA <- downloadHandler(
    filename = function() { 
      paste("metadataTableRNA-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_metadata_RNA, file, sep = "\t", quote = F, row.names = F)
    })
    
  
  #------------------Utilities------------------------------------------
  #save rna object as rds file
  output$utilitiesConfirmExport1 <- downloadHandler(
    filename = function() { 
      paste("processed_seurat_object-", Sys.Date(), ".RDS", sep="")
    },
    content = function(file) {
      saveRDS(seurat_object, file)
    })
  
  output$utilitiesConfirmExport2 <- downloadHandler(
    filename = function() { 
      paste("processed_seurat_object-", Sys.Date(), ".RDS", sep="")
    },
    content = function(file) {
      saveRDS(seurat_object, file)
    })
  
  output$utilitiesConfirmExport <- downloadHandler(
    filename = function() { 
      paste("processed_seurat_object-", Sys.Date(), ".RDS", sep="")
    },
    content = function(file) {
      saveRDS(seurat_object, file)
    })
  
  #rename-merge-delete clusters
  observeEvent(input$utilitiesConfirmRename, {
    if(!identical(seurat_object, NULL) & input$utilitiesRenameOldName != "-" & input$utilitiesRenameNewName != "")
    {
      showModal(modal_confirm_utilities_rename) 
    }
  })
  
  observeEvent(input$utilitiesConfirmDelete, {
    if(!identical(seurat_object, NULL) & input$utilitiesDeleteCluster != "-")
    {
      showModal(modal_confirm_utilities_delete)
    }
  })
  
  #change active clusters
  observeEvent(input$utilitiesConfirmChangeCluster, {
    tryCatch({
    seurat_object$seurat_clusters <<- seurat_object@meta.data[, input$utilitiesActiveClusters]
    Idents(seurat_object) <<- seurat_object$seurat_clusters
    
    updateInputLRclusters()
    updateMetadata()
    updateClusterTab()
    shinyalert(title = "Active clustering column changed succesfully!", type = "success")
    }, error = function(e) {
      print(paste("Error :  ", e))
      shinyalert(title = "There was an error with changing the clustering column.", type = "error")
    }, finally = {
      
    })
  })
  
  #change active assay 
  observeEvent(input$utilitiesConfirmChangeAssay, {
    tryCatch({
      DefaultAssay(seurat_object) <<- input$utilitiesActiveAssay
      updateGeneSearchFP()
      shinyalert(title = "Default assay was changed succesfully!", type = "success")
    }, error = function(e) {
      print(paste("Error :  ", e))
      shinyalert(title = "There was an error with changing the default assay.", type = "error")
    }, finally = {
      
    })
    #print("def assay")
    #print(DefaultAssay(seurat_object))
  })
  
  observeEvent(input$utilitiesPlotCommands, {
    if(!is.null(seurat_object))
    {
      if(length(seurat_object@commands) > 0)
      {
        output$history <- renderPrint(
          {
            for(i in 1:length(seurat_object@commands))
            {
              print(seurat_object@commands[[i]])
              print("-------------------------------------------------------------")
            }
          })
      }
      else
      {
        Sys.sleep(1)
        shinyalert("Error. No command history was found.", type = "error")
      }
    }
    else
    {
      Sys.sleep(1)
      shinyalert("Error. No project is loaded.", type = "error")
    }
  })
  
  #------------------Quality Control tab--------------------------------
  observeEvent(input$qcDisplay, {
    start.time <- Sys.time()
    print(paste0("!!!!-----DisplayQC_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "qcDisplay")
    tryCatch({
      if (identical(init_seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else{
        shinyjs::show("nFeatureViolin_loader")
        shinyjs::show("totalCountsViolin_loader")
        shinyjs::show("mitoViolin_loader")
        shinyjs::show("genesCounts_loader")
        shinyjs::show("mtCounts_loader")
        showModal(modalDialog(div('Quality control in progress. This operation may take a few minutes for large datasets. Please wait...', style='color:#ffffff; background-color:#222d32;'), 
                              footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        
        output$nFeatureViolin <- renderPlotly(
          {
            p <- VlnPlot(init_seurat_object, features = c("nFeature_RNA"), pt.size = 0, group.by = "orig.ident",
                         cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(init_seurat_object@meta.data[, 'orig.ident'])))) + 
              theme_bw() + 
              geom_hline(yintercept=c(as.numeric(input$minUniqueGenes), as.numeric(input$maxUniqueGenes)), linetype="dashed", color = "red", size=1) + 
              theme(
                plot.title = element_blank(),
                axis.title.x = element_blank(),
                legend.position = "none") +
              labs(title = "", y="Genes detected/cell")
            
            p_temp <- VlnPlot(init_seurat_object, features = c("nFeature_RNA"), pt.size = 0, group.by = "orig.ident",
                              cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(init_seurat_object@meta.data[, 'orig.ident'])))) + 
              theme_bw() + 
              geom_hline(yintercept=c(as.numeric(input$minUniqueGenes), as.numeric(input$maxUniqueGenes)), linetype="dashed", color = "red", size=1) + 
              theme(
                plot.title = element_blank(),
                axis.title.x = element_blank(),
                legend.position = "none") +
              labs(title = "", y="Genes detected/cell")
            
            plotly::ggplotly(p + ylim(min(init_seurat_object$nFeature_RNA)-1, max(init_seurat_object$nFeature_RNA)+2), tooltip = c("x", "y")) 
          }
        )
        
        output$totalCountsViolin <- renderPlotly(
          {
            p <- VlnPlot(init_seurat_object, features = c("nCount_RNA"), pt.size = 0, group.by = "orig.ident",
                         cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(init_seurat_object@meta.data[, 'orig.ident'])))) + 
              theme_bw() + 
              theme(
                plot.title = element_blank(),
                axis.title.x = element_blank(),
                legend.position = "none")+
              labs(title = "", y="Total counts/cell")
            plotly::ggplotly(p, tooltip = c("x", "y")) 
          }
        )
        
        output$mitoViolin <- renderPlotly(
          {
            p <- VlnPlot(init_seurat_object, features = c("percent.mt"), pt.size = 0, group.by = "orig.ident",
                         cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(init_seurat_object@meta.data[, 'orig.ident'])))) + 
              theme_bw() + 
              #scale_y_continuous(labels = number_format(accuracy = 0.001)) +
              geom_hline(yintercept= as.numeric(input$maxMtReads), linetype="dashed", color = "red", size=1) +
              theme(
                plot.title = element_blank(),
                axis.title.x = element_blank(),
                legend.position = "none")+
              labs(title = "", y="% of reads mapped to mitochonrial genome/cell")
            plotly::ggplotly(p, tooltip = c("x", "y")) #%>%
              #layout(yaxis = list(hoverformat = ".2f"))
          }
        )
        
        output$mtCounts <- renderPlotly(
          {
            p <- FeatureScatter(init_seurat_object, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "orig.ident", raster = F,
                                cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(init_seurat_object@meta.data[, 'orig.ident']))))
            gp <- plotly::ggplotly(p)
            print(gp)
          }
        )
        
        output$genesCounts <- renderPlotly(
          {
            p <- FeatureScatter(init_seurat_object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "orig.ident", raster = F,
                                cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(init_seurat_object@meta.data[, 'orig.ident']))))
            gp <- plotly::ggplotly(p)
            print(gp)
          }
        )
        
        output$cellStats <- renderPrint(
          {
            cat(paste0("Total number of cells: ", nrow(init_seurat_object@meta.data)))
          }
        )
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "The selected Quality Control arguments cannot produce meaningful visualizations.")
    }, finally = {
      removeModal()
      
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----DisplayQC_end_Mins_taken-----", time.taken))
      
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "qcDisplay")
    })
  })
  
  observeEvent(input$qcConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----QCFiltering_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "qcConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else{
        updateTabsetPanel(session, inputId = "qc_tabs_rna", selected = "Post-filtering plots")
        showModal(modalDialog(div('Analysis in Progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("filteredNFeatureViolin_loader")
        shinyjs::show("filteredTotalCountsViolin_loader")
        shinyjs::show("filteredMitoViolin_loader")
        shinyjs::show("filteredGenesCounts_loader")
        shinyjs::show("filteredMtCounts_loader")
        
        qc_minFeatures <<- input$minUniqueGenes
        qc_maxFeatures <<- input$maxUniqueGenes
        qc_maxMtPercent <<- input$maxMtReads
        
        seurat_object <<- subset(init_seurat_object, subset = nFeature_RNA > as.numeric(qc_minFeatures) & nFeature_RNA < as.numeric(qc_maxFeatures) & percent.mt < as.double(qc_maxMtPercent)) #filter object
        
        output$filteredNFeatureViolin <- renderPlotly(
          {
            p <- VlnPlot(seurat_object, features = c("nFeature_RNA"), pt.size = 0, group.by = input$qcColorBy,
                         cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, input$qcColorBy])))) + 
              theme_bw() + 
              theme(
                plot.title = element_blank(),
                axis.title.x = element_blank(),
                legend.position = "none") +
              labs(title = "", y="Genes detected/cell")
            plotly::ggplotly(p, tooltip = c("x", "y")) 
          }
        )
        output$filteredTotalCountsViolin <- renderPlotly(
          {
            p <- VlnPlot(seurat_object, features = c("nCount_RNA"), pt.size = 0, group.by = input$qcColorBy,
                         cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, input$qcColorBy])))) + 
              theme_bw() + 
              theme(
                plot.title = element_blank(),
                axis.title.x = element_blank(),
                legend.position = "none")+
              labs(title = "", y="Total counts/cell")
            plotly::ggplotly(p, tooltip = c("x", "y")) 
          }
        )
        output$filteredMitoViolin <- renderPlotly(
          {
            p <- VlnPlot(seurat_object, features = c("percent.mt"), pt.size = 0, group.by = input$qcColorBy,
                         cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, input$qcColorBy])))) + 
              theme_bw() + 
              theme(
                plot.title = element_blank(),
                axis.title.x = element_blank(),
                legend.position = "none")+
              labs(title = "", y="% of reads mapped to mitochonrial genome/cell")
            plotly::ggplotly(p, tooltip = c("x", "y"))
          }
        )
        output$filteredMtCounts <- renderPlotly(
          {
            p <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = input$qcColorBy, raster = F, 
                                cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, input$qcColorBy]))))
            gp <- plotly::ggplotly(p)
            print(gp)
          }
        )
        output$filteredGenesCounts <- renderPlotly(
          {
            p <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = input$qcColorBy, raster = F,
                                cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, input$qcColorBy]))))
            gp <- plotly::ggplotly(p)
            print(gp)
          }
        )
        output$filteredCellStats <- renderPrint(
          {
            cat(paste0("\nTotal number of cells after filtering: ", nrow(seurat_object@meta.data)))
          }
        )
        updateSelInpColor()
        session$sendCustomMessage("handler_disableTabs", "sidebarMenu")
        updateMetadata()
        cleanAllPlots(F) # fromDataInput -> FALSE
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " QUALITY CONTROL", " DATA NORMALIZATION\n& SCALING"))
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "The selected Quality Control arguments cannot produce meaningful visualizations.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----QCfilter_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "qcConfirm")
    })
  })
  
  #ATAC qc
  observeEvent(input$qcDisplayATAC, {
    start.time <- Sys.time()
    print(paste0("!!!!-----DisplayQC_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "qcDisplayATAC")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else{
        showModal(modalDialog(div('Quality control metrics are prepared. Rendering mighr take several minutes for large datasets. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("TSS_plot_loader")
        shinyjs::show("nFrag_plot_loader")
        shinyjs::show("TSS_nFrag_plot_loader")
    #TSS plot
    p1 <- plotGroups(
      ArchRProj = proj_default,
      colorBy = "cellColData",
      name = "TSSEnrichment",
      plotAs = "violin",
      alpha = 0.4,
      addBoxPlot = TRUE
    )
    output$TSS_plot <- renderPlotly( expr =  ggplot(p1$data, aes(x=x, y=y, fill=x)) + geom_violin() + theme_bw() + labs(y="TSS Enrichment"))
    
    #nFrags plot
    p2 <- plotGroups(
      ArchRProj = proj_default,
      colorBy = "cellColData",
      name = "log10(nFrags)",
      plotAs = "ridges"
    )
    output$nFrag_plot <- renderPlot( expr = ggplot(p2$data, aes(y=x, x=y, fill=x)) + geom_density_ridges() + theme_bw() + scale_y_discrete(expand = c(0, 0)) + labs(x="log10(nFrags)"), 
                                             width = 500, height = 500
                                          )
    
    #nFrags-TSS plot
    df <- getCellColData(proj_default, select = c("log10(nFrags)", "TSSEnrichment"))
    p4 <- ggPoint(
      x = df[,1],
      y = df[,2],
      #size = 3,
      #baseSize=30,
      #legendSize=5,
      #dpi=1200,
      colorDensity = TRUE,
      continuousSet = "sambaNight",
      xlabel = "Log10 Unique Fragments",
      ylabel = "TSS Enrichment",
      xlim = c(log10(500), quantile(df[,1], probs = 0.99)),
      ylim = c(0, quantile(df[,2], probs = 0.99))
    ) + geom_hline(yintercept = 4, lty = "dashed") + geom_vline(xintercept = 3, lty = "dashed")
    output$TSS_nFrag_plot <- renderPlotly( expr =  ggplotly(p4))
    
    #cells 
    output$CellStatsATAC <- renderPrint(
      {
        cat(paste0("\nTotal number of cells after soft filtering: ", nrow(getCellColData(proj_default))))
      }
    )
    
    session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " QUALITY CONTROL", " PCA/LSI"))
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "The selected Quality Control arguments cannot produce meaningful visualizations.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----DisplayQC_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "qcDisplayATAC")
    })
  })
  
  #------------------Normalization tab--------------------------------
  observeEvent(input$normalizeConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----Scaling_start-----", start.time))
    
    session$sendCustomMessage("handler_log", " ### Starting normalization procedure ###")
    session$sendCustomMessage("handler_disableAllButtons", "normalizeConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else{
        shinyjs::show("hvgScatter_loader")
        showModal(modalDialog(div('Normalization and scaling in progress. This operation will be dramatically slowed down if all genes are selected.
                                  Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 
                                'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        normalize_normMethod <- "LogNormalize"
        normalize_normScaleFactor <- input$normScaleFactor
        seurat_object <<- NormalizeData(seurat_object, normalization.method = normalize_normMethod, scale.factor = as.numeric(normalize_normScaleFactor))
        session$sendCustomMessage("handler_log", "Finished performing log-normalization.")
        
        normalize_hvgMethod <<- input$radioHVG
        normalize_hvgNGenes <<- input$nHVGs

        seurat_object <<- FindVariableFeatures(seurat_object, selection.method = normalize_hvgMethod, nfeatures = as.numeric(normalize_hvgNGenes))
        session$sendCustomMessage("handler_log", "Finished calculating gene and feature variances of standardized and clipped values.")
        
        normalize_scaleRegressOut <- input$normalizeRegressColumns
        
        all.genes <- VariableFeatures(seurat_object)
        if(input$normalizeScaleGenes == "all_genes")
        {
          all.genes <- rownames(seurat_object)
        }
        
        
        if(is.null(normalize_scaleRegressOut)) seurat_object <<- ScaleData(seurat_object, features = all.genes) 
        else seurat_object <<- ScaleData(seurat_object, vars.to.regress=normalize_scaleRegressOut, features = all.genes)
        
        session$sendCustomMessage("handler_log", "Finished centering and scaling data matrix.")
        
        # Rendering
        output$hvgScatter <- renderPlotly({ # tooltip
          if(length(VariableFeatures(seurat_object)) != 0){
            plot1 <- VariableFeaturePlot(seurat_object)
            varplot <- plot1$data
            varplot$gene <- rownames(varplot)
            varplot$colors[varplot$colors == "yes"] <- paste0("Highly Variable genes(", length(VariableFeatures(seurat_object)), ")")
            varplot$colors[varplot$colors == "no"] <- paste0("Not Variable genes(", length(rownames(seurat_object)) - length(VariableFeatures(seurat_object)), ")")
            
            if(normalize_hvgMethod == "vst"){
              p <- ggplot(varplot, aes(x=log10(mean), y=variance.standardized, color=colors, text = paste0("log10(mean expression): ", log10(mean),
                                                                                                           "\nstandardized variance: ", variance.standardized,
                                                                                                           "\ngene: ", gene)
              )) +
                geom_point() +
                theme_bw() +
                scale_color_manual(
                  values = c("red", "black")
                )+
                labs(x="Average Expression", y="Standardized Variance", color="")
            } else {
              p <- ggplot(varplot, aes(x=mean, y=dispersion, color=colors, text = paste0("mean expression: ", mean,
                                                                                                "\ndispersion: ", dispersion,
                                                                                                "\ngene: ", gene)
              )) +
                geom_point() +
                theme_bw() +
                scale_color_manual(
                  values = c("red", "black")
                )+ 
                labs(x="Average Expression", y="Dispersion", color="")
            }
            gp <- plotly::ggplotly(p, tooltip = "text")
            return(gp)
          } else session$sendCustomMessage("handler_alert", "There are no variable features in the Seaurat object.")
        }) # End Rendering
        
        output$hvgTop10Stats <- renderPrint(
          {
            top10hvgs <- head(VariableFeatures(seurat_object), 10)
            cat("\nTop10 varibale features:", c(top10hvgs[1:10]))
          })
        
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " PCA/LSI"))
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "The selected Normalization arguments cannot produce meaningful visualizations.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Scaling_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "normalizeConfirm")
      session$sendCustomMessage("handler_log", " ### Finished normalization procedure ###")
    })
  }) 
  
  #------------------PCA tab------------------------------------------
  observeEvent(input$PCrunPCA, {
    start.time <- Sys.time()
    print(paste0("!!!!-----PCA_start-----", start.time))
    
    session$sendCustomMessage("handler_log", " ### Starting PCA Analysis ###")
    session$sendCustomMessage("handler_disableAllButtons", "PCrunPCA")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (!"ScaleData.RNA" %in% names(seurat_object@commands)) session$sendCustomMessage("handler_alert", "Data need to be scaled first. Please, execute the previous step in DATA NORMALIZATION & SCALING.")
      else {
        showModal(modalDialog(div('PCA analysis in progress. Automatic optimum value calculation will increase the runtime of this step
                                  dramatically. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("elbowPlotPCA_loader")
        shinyjs::show("PCAscatter_loader")
        
        seurat_object <<- RunPCA(seurat_object, features = VariableFeatures(object = seurat_object), npcs = 100)
        
        optimal_nPCs <- 0 #only for visualization purposes
        if(input$pcaRadio == "yes")
        {
          
          data.use <- PrepDR(seurat_object, genes.use = VariableFeatures(seurat_object), use.imputed = F, assay.type = "RNA")
          initial_optimal_nPCs <- PCA_estimate_nPC(data.use, from.nPC = 1, to.nPC=100, by.nPC= 10, k = 10)
          
          optimal_nPCs <- PCA_estimate_nPC(data.use, from.nPC = 1, to.nPC=initial_optimal_nPCs+10, by.nPC= 1, k = 10)
          
          #Update selected PCs 
          updatePCs_selection(optimal_nPCs)
        }
        
        updateUmapTypeChoices("pca")
        
        output$elbowPlotPCA <- renderPlotly(
          {
            plot1 <- ElbowPlot(seurat_object, ndims = 100)
            plot1_data <- plot1$data
            colnames(plot1_data)[1] <- "PC"
            colnames(plot1_data)[2] <- "SD"
            
            if(input$pcaRadio == "yes" & optimal_nPCs != 0)
            {
              p <- ggplot(plot1_data, aes(x=PC, y=SD))+
                geom_point() +
                theme_bw() +
                geom_vline(xintercept= optimal_nPCs, linetype="dashed", color = "red", size=1) +
                labs(x="PC", y="Standard Deviation")
            }
            else
            {
              p <- ggplot(plot1_data, aes(x=PC, y=SD))+
                geom_point() +
                theme_bw() +
                labs(x="PC", y="Standard Deviation")
            }
            
            gp <- plotly::ggplotly(p, tooltip = c("x", "y"))
            print(gp)
          }
        )
        
        output$PCAscatter <- renderPlotly(
          {
            #prepare metadata
            meta <- seurat_object@meta.data
            meta$Cell_id <- rownames(meta)
            meta <- meta[, ]#meta[, c('Cell_id', 'seurat_clusters', 'orig.ident')]
            reduc_data <- data.frame()
            
            #prepare colors
            cols = colorRampPalette(brewer.pal(12, "Paired"))(1)
            
            #umap data frame
            seurat_object_reduc <- as.data.frame(seurat_object@reductions$pca@cell.embeddings)
            seurat_object_reduc <- seurat_object_reduc[, c(1:ncol(seurat_object_reduc))]
            seurat_object_reduc$Cell_id <- rownames(seurat_object_reduc)
            reduc_data <- left_join(seurat_object_reduc, meta)
            
            p <- ggplot(data=reduc_data, aes_string(x="PC_1", y="PC_2", fill="orig.ident")) +
              geom_point(size=2, shape=19, stroke=0)+
              scale_fill_manual(values = cols)+
              scale_size()+
              theme_bw() +
              theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                    axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                    axis.title.y = element_text(face = "bold", color = "black", size = 25),
                    axis.title.x = element_text(face = "bold", color = "black", size = 25),
                    legend.text = element_text(face = "bold", color = "black", size = 9),
                    legend.title = element_text(face = "bold", color = "black", size = 9),
                    legend.position="right",
                    title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
              labs(x="PC 1", y="PC 2", title = "", fill="Color")
            print(plotly::ggplotly(p))
          }
        )
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " CLUSTERING", " CELL CYCLE PHASE ANALYSIS", " FEATURE INSPECTION"))
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the PCA analysis.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----PCA_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "PCrunPCA")
      session$sendCustomMessage("handler_log", " ### Finished PCA Analysis ###")
    })
  })
  
  observeEvent(input$PCconfirm, {
    session$sendCustomMessage("handler_log", " ### Starting PC Exploration ###")
    session$sendCustomMessage("handler_disableAllButtons", "PCconfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (!"pca" %in% names(seurat_object)) session$sendCustomMessage("handler_alert", "Please, first Run PCA above.")
      else {
        shinyjs::show("PCAloadings_loader")
        shinyjs::show("PCAheatmap_loader")
        showModal(modalDialog(div('Rendering in progress, please wait...', style='color:#ffffff; background-color:#222d32;'), 
                              footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        
        activePC <- as.numeric(input$PCin)
        
        output$PCAloadings <- renderPlotly(
          {
            plot1 <- VizDimLoadings(seurat_object, dims = activePC, reduction = "pca", balanced = TRUE)
            plot1_data <- plot1$data
            plot1_data$orig.ident <- unique(seurat_object@meta.data$orig.ident)
            colnames(plot1_data)[1] <- "LoadingScore"
            activePC <- paste0("PC", "_", activePC)
            
            export_loadingScoresTable_RNA <<- plot1$data
            
            p <- ggplot(plot1_data, aes(x=LoadingScore, y=feature, color=orig.ident))+
              geom_point() +
              theme_bw() +
              scale_color_manual(values="blue")+
              theme(legend.position = "none")+
              labs(x=activePC, y="")
            
            gp <- plotly::ggplotly(p, tooltip = c("x", "y"))
            print(gp)
          }
        )
        
        output$PCAheatmap <- renderPlotly(
          {
            p <- DimHeatmap(seurat_object, dims = activePC, cells = 500, balanced = TRUE, fast = F)
            p
            plotly::ggplotly(p) 
          }
        )
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the visualization of PCA analysis.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "PCconfirm")
      session$sendCustomMessage("handler_log", " ### Finished PC Exploration ###")
    })
  })
  
  output$pcaRNAExport <- downloadHandler(
    filename = function() { 
      paste("loadingScoresTableRNA-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_loadingScoresTable_RNA, file, sep = "\t", quote = F, row.names = F)
    })
  
  #ATAC LSI
  observeEvent(input$lsiConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----LSI_start-----", start.time))
    
    session$sendCustomMessage("handler_log", " ### Starting LSI ###")
    session$sendCustomMessage("handler_disableAllButtons", "lsiConfirm")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else
      {
        showModal(modalDialog(div('LSI Analysis in progress. High number of iterations increase the runtime of the process dramatically. 
                                  Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        proj_default <<- addIterativeLSI(ArchRProj = proj_default, useMatrix = "TileMatrix", name = "IterativeLSI", force = T,
                                         iterations = as.numeric(input$lsiIterations), varFeatures = as.numeric(input$lsiVarFeatures),
                                         clusterParams = list( resolution = as.numeric(input$lsiResolution), sampleCells = 10000, n.start = 10),dimsToUse=1:as.numeric(input$lsiDmensions), 
                                         logFile = paste0(user_dir, "/LSI_file.log"))
        saveArchRProject(proj_default)
        output$lsiOutput <- renderPrint(
          {
            cat(paste0("LSI executed successfully."))
          }
        )
      }
      session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " CLUSTERING"))
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the LSI analysis.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----LSI_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "lsiConfirm")
      session$sendCustomMessage("handler_log", " ### Finished LSI ###")
    })
  })
  
  #------------------Clustering tab------------------------------------------
  observeEvent(input$snnConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----Clustering_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "snnConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (!"pca" %in% names(seurat_object)) session$sendCustomMessage("handler_alert", "Please, first execute PRINCIPAL COMPONENT ANALYSIS.")
      else {
        showModal(modalDialog(div('Clustering cells in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        snn_dims <<- input$snnPCs
        snn_k <<- input$snnK
        cluster_res <<- input$clusterRes
        cluster_dims <<- input$clusterPCs
        
        seurat_object <<- FindNeighbors(seurat_object, k.param = as.numeric(snn_k), dims = 1:as.numeric(snn_dims), reduction = "pca", nn.method="annoy")
        seurat_object <<- FindClusters(seurat_object, resolution = as.numeric(cluster_res))
        
        updateClusterTab()
        updateSelInpColor()
        updateInputLRclusters()
        updateMetadata()
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " ADDITIONAL DIMENSIONALITY\nREDUCTION METHODS", " TRAJECTORY ANALYSIS", " DOUBLETS' DETECTION",
                                                          " MARKERS' IDENTIFICATION", " LIGAND - RECEPTOR\nANALYSIS", " UTILITY OPTIONS"))
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the clustering procedure.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Clustering_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "snnConfirm")
    })
  })
  
  observeEvent(input$clusterBarplotConfirm, {
    session$sendCustomMessage("handler_disableAllButtons", "clusterBarplotConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first.")
      else {
        showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("clusterBarplot_loader")
        
        if(input$clusterGroupBy == "seurat_clusters")
        {
          #barplot for cell distribution per cluster
          clusterTable <- as.data.frame(table(seurat_object$seurat_clusters)) #as.data.frame(table(Idents(seurat_object)))
          totalCells <- sum(clusterTable$Freq)
          
          clusterTable$Perc <- (clusterTable$Freq)/totalCells
          colnames(clusterTable)[1] <- "Cluster"
          colnames(clusterTable)[2] <-  "Cells"
          colnames(clusterTable)[3] <- "Percentage"
          cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(clusterTable$Cluster)))
          
          p <- ggplot(clusterTable) + theme_bw() +
            geom_bar( mapping = aes(x = Cluster, y = Percentage, fill=Cluster), stat = "identity" ) +
            scale_fill_manual(values = cols ) +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 12, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 12, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 12),
                  axis.title.x = element_text(face = "bold", color = "black", size = 12),
                  panel.background = element_blank(),
                  axis.line = element_line(colour = "black")) +
            labs(x="", y="Percent of cells", fill="Clusters")
          gp <- plotly::ggplotly(p, tooltip = c("x", "y"))
          
          output$clusterBarplot <- renderPlotly({print(gp)})
        }
        else
        {
          cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object$seurat_clusters)))
          
          p <- dittoBarPlot(seurat_object, var = "seurat_clusters",group.by = input$clusterGroupBy, scale = "percent") + theme_bw() +
            scale_fill_manual(values = cols ) +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 12, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 12, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 12),
                  axis.title.x = element_text(face = "bold", color = "black", size = 12),
                  panel.background = element_blank(),
                  axis.line = element_line(colour = "black")) 
          
          gp <- plotly::ggplotly(p, tooltip = c("x", "y"))
          
          output$clusterBarplot <- renderPlotly({print(gp)})
        }
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error during plotting the barplot.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "clusterBarplotConfirm")
    })
  })
  
  observeEvent(input$snnDisplayConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----SNN_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "snnDisplayConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else{
        showModal(modalDialog(div('SNN graph visualization in progress. This a slow operation for large datasets. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("snnSNN_loader")
        
        output$snnSNN <- renderVisNetwork(
          {
            if(!is.null(seurat_object@graphs$RNA_snn))
            {
              #set.seed(9)
              
              mygraph <- as.matrix(seurat_object@graphs$RNA_snn)
              graphOut <- graph_from_adjacency_matrix(mygraph, mode = "undirected", weighted = T, diag = F)
              graphSimple <- simplify(graphOut) #, remove.loops=T)
              weights <- E(graphSimple)$weight
              sub_nodes <- V(graphSimple)$name
              
              tableCl <- seurat_object@meta.data[, ]
              tableCl$Cell_id <- rownames(tableCl)
              tableCl <- tableCl[, c('Cell_id', 'seurat_clusters')]
              tableCl <- tableCl[order(tableCl$seurat_clusters), ]
              tableCl <- tableCl[sub_nodes, ]
              colors_cl <- colorRampPalette(brewer.pal(12, "Paired"))(length(unique(tableCl$seurat_clusters)))
              tableCl$color <- colors_cl[as.numeric(tableCl$seurat_clusters)]
              
              V(graphSimple)$color <- tableCl$color
              visIgraph(graphSimple, layout = "layout_with_lgl", randomSeed = 9) %>%
                visInteraction(navigationButtons = TRUE, hover = TRUE)
            }
          }
        )
        
        
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was a problem in drawing the SNN graph.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----SNN_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "snnDisplayConfirm")
    })
  }) 
  
  #ATAC clustering
  observeEvent(input$clusterConfirmATAC, {
    start.time <- Sys.time()
    print(paste0("!!!!-----Clustering_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "clusterConfirmATAC")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        showModal(modalDialog(div('Cell clustering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("clusterBarplotATAC_loader")
    #Run clustering
    proj_default <<- addClusters(input = proj_default, reducedDims = "IterativeLSI", method = "Seurat", name = "Clusters", 
                                 resolution = as.numeric(input$clusterResATAC), dimsToUse = 1:as.numeric(input$clusterDimensionsATAC), force = T, 
                                 logFile =paste0(user_dir, "/Clusters_file.log"))
    saveArchRProject(proj_default)
    #Cluster table
    cluster_df <- as.data.frame(table(proj_default$Clusters))
    colnames(cluster_df)[1] <- "Cluster"
    colnames(cluster_df)[2] <- "Number of cells"
    cluster_df$`% of cells per cluster` <- as.double(sprintf("%.3f", as.double(cluster_df$`Number of cells`/nrow(getCellColData(proj_default))*100)))
    
    output$clusterTableATAC <- renderDataTable(cluster_df, options = list(pageLength = 10), rownames = F)
    export_clustertable_ATAC <<- cluster_df
    
    #Impute genescores
    proj_default <<- addImputeWeights(proj_default,sampleCells = 5000, 
                                      logFile =paste0(user_dir, "/Impute_weights_file.log"))
    updateSelectizeInput(session, "findMarkersGeneSelectATAC", choices = unique(proj_default@geneAnnotation$genes$symbol), server = T)
    
    cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(cluster_df$Cluster)))
    
    p <- ggplot(cluster_df) + theme_bw() +
      geom_bar( mapping = aes(x = Cluster, y = `% of cells per cluster`, fill=Cluster), stat = "identity" ) +
      scale_fill_manual(values = cols ) +
      theme(axis.text.x = element_text(face = "bold", color = "black", size = 12, angle = 0),
            axis.text.y = element_text(face = "bold", color = "black", size = 12, angle = 0),
            axis.title.y = element_text(face = "bold", color = "black", size = 12),
            axis.title.x = element_text(face = "bold", color = "black", size = 12),
            panel.background = element_blank(),
            axis.line = element_line(colour = "black")) +
      labs(x="", y="Percent of cells", fill="Clusters")
    gp <- plotly::ggplotly(p, tooltip = c("x", "y"))
    
    output$clusterBarplotATAC <- renderPlotly({print(gp)})
    updateMetadataATAC()
    updateSelInpColorATAC()
    session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " ADDITIONAL DIMENSIONALITY\nREDUCTION METHODS", " TRAJECTORY ANALYSIS", 
                                                      " MARKERS' IDENTIFICATION"))
    session$sendCustomMessage("handler_disableAllButtons", "umapConfirmATAC")
    }
      }, error = function(e) {
        print(paste("Error :  ", e))
        session$sendCustomMessage("handler_alert", "There was an error with the clustering procedure.")
      }, finally = {
        end.time <- Sys.time()
        time.taken <- as.numeric (end.time - start.time, units = "mins")
        print(paste0("!!!!-----Clustering_end_Mins_taken-----", time.taken))
        
        removeModal()
        Sys.sleep(1)
        session$sendCustomMessage("handler_finishLoader", "clust4_loader")
        session$sendCustomMessage("handler_enableAllButtons", "clusterConfirmATAC")
      })
  })
  
  output$clusterTableRNAExport <- downloadHandler(
    filename = function() { 
      paste("clusterTableRNA-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_clustertable_RNA, file, sep = "\t", quote = F, row.names = F)
    })
  
  output$clusterTableExportATAC <- downloadHandler(
    filename = function() { 
      paste("clusterTableATAC-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_clustertable_ATAC, file, sep = "\t", quote = F, row.names = F)
    })
  
  #------------------Umap/tSNE/DFM tab---------------------------------------
  observeEvent(input$umapRunUmapTsneATAC, {
    start.time <- Sys.time()
    print(paste0("!!!!-----UMAP/tsne_start-----", start.time))
    
    
    session$sendCustomMessage("handler_disableAllButtons", "umapRunUmapTsneATAC")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        showModal(modalDialog(div('UMAP and tSNE calculation. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;'))
        ## Dimensionality Reduction ##
        proj_default <<- addUMAP(ArchRProj = proj_default, reducedDims = "IterativeLSI", name = "umap", nNeighbors = 30, minDist = 0.5, metric = "cosine", seed = 4,
                                 force = T, n_components = as.numeric(input$umapOutComponentsATAC), dimsToUse = 1:as.numeric(input$umapDimensionsATAC))
        
        #UMAP is used in Trajectory tab only
        proj_default <<- addUMAP(ArchRProj = proj_default, reducedDims = "IterativeLSI", name = "UMAP", nNeighbors = 30, minDist = 0.5, metric = "cosine", seed = 4,
                                 force = T, n_components = as.numeric(input$umapOutComponentsATAC), dimsToUse = 1:as.numeric(input$umapDimensionsATAC)) #3D-->2
        
        df_2d <- proj_default@embeddings$UMAP$df[, 1:2]
        proj_default@embeddings$UMAP$df <- df_2d
        
        proj_default <<- addTSNE(ArchRProj = proj_default, reducedDims = "IterativeLSI", name = "tsne", perplexity = 30, force = T, 
                                 n_components = as.numeric(input$umapOutComponentsATAC), dimsToUse = 1:as.numeric(input$umapDimensionsATAC))
        
        saveArchRProject(proj_default)
        updateReductionATAC()
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " DOUBLETS' DETECTION", " FEATURE INSPECTION", " CLUSTERS' ANNOTATION"))
        }
      }, error = function(e) {
        print(paste("Error :  ", e))
        session$sendCustomMessage("handler_alert", "There was an error with the generation of UMAP or tSNE.")
      }, finally = {
        end.time <- Sys.time()
        time.taken <- as.numeric (end.time - start.time, units = "mins")
        print(paste0("!!!!-----Umap/tsne_end_Mins_taken-----", time.taken))
        
        removeModal()
        Sys.sleep(1)
        session$sendCustomMessage("handler_enableAllButtons", "umapRunUmapTsneATAC")
      })
  })
  
  updateReductionATAC <- function()
  {
    session$sendCustomMessage("handler_disableAllButtons", "umapConfirmATAC")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        shinyjs::show("umapPlotATAC_loader")
        #get input
        dims <- as.numeric(input$umapDimensionsPlotATAC)
        type <- input$umapTypeATAC
        
        #prepare metadata
        meta <- as.data.frame(getCellColData(proj_default))
        meta$Cell_id <- rownames(meta)
        reduc_data <- data.frame()
        
        #prepare colors
        cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorByATAC])))
        
        #for all reductions
        archr_object_reduc <- proj_default@embeddings[[input$umapTypeATAC]]$df
        archr_object_reduc <- archr_object_reduc[, c(1:ncol(archr_object_reduc))]
        archr_object_reduc$Cell_id <- rownames(archr_object_reduc)
        reduc_data <- left_join(archr_object_reduc, meta)
        print(head(reduc_data))
        
        colnames(reduc_data) <- gsub("IterativeLSI#", "", colnames(reduc_data))
        colnames(reduc_data) <- gsub("Dimension_", "", colnames(reduc_data))
        
        if(type == "umap" & dims == 2 & dims <= ncol(archr_object_reduc)-1 )
        {
          p <- ggplot(data=reduc_data, aes_string(x="UMAP_1", y="UMAP_2", fill=input$umapColorByATAC)) +
            geom_point(size= as.numeric(input$umapDotSizeATAC), shape=21, alpha= as.numeric(input$umapDotOpacityATAC), stroke=as.numeric(input$umapDotBorderATAC))+
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="UMAP 1", y="UMAP 2", color="Cell type", title = "", fill="Color")
          output$umapPlotATAC <- renderPlotly({plotly::ggplotly(p)})
        }
        else if(type == "umap" & dims == 3 & dims <= ncol(archr_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~UMAP_1, y=~UMAP_2, z=~UMAP_3, type="scatter3d", alpha = as.numeric(input$umapDotOpacityATAC), mode="markers", color=as.formula(paste0('~', input$umapColorByATAC)),
                       marker = list(size = as.numeric(input$umapDotSizeATAC), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorderATAC))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorByATAC]))) )
          
          output$umapPlotATAC <- renderPlotly({print(p)})
        }
        else if(type == "tsne" & dims == 2 & dims <= ncol(archr_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="TSNE_1", y="TSNE_2", fill=input$umapColorByATAC)) +
            geom_point(size= as.numeric(input$umapDotSizeATAC), shape=21, alpha= as.numeric(input$umapDotOpacityATAC), stroke=as.numeric(input$umapDotBorderATAC)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="tSNE 1", y="tSNE 2", color="Cell type", title = "", fill="Color")
          output$umapPlotATAC <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "tsne" & dims == 3 & dims <= ncol(archr_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~TSNE_1, y=~TSNE_2, z=~TSNE_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacityATAC), color=as.formula(paste0('~', input$umapColorByATAC)), 
                       marker = list(size = as.numeric(input$umapDotSizeATAC), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorderATAC))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorByATAC]))) ) 
          output$umapPlotATAC <- renderPlotly({print(p)})
        }
        else
        {
          session$sendCustomMessage("handler_alert", "There was an error with the number of dimensions of the plot. Please set Dimensions to 2D and update the plot.")
        }
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error during of UMAP or tSNE plotting.")
    }, finally = {
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "umapConfirmATAC")
    })
  }
  
  observeEvent(input$umapConfirmATAC, {
    updateReductionATAC()
  })
  
  observeEvent(input$umapRunUmap, {
    start.time <- Sys.time()
    print(paste0("!!!!-----UMAP_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "umapRunUmap")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunTsne")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunDFM")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunPhate")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (!"pca" %in% names(seurat_object)) session$sendCustomMessage("handler_alert", "Please, first execute PRINCIPAL COMPONENT ANALYSIS.")
      else {
        showModal(modalDialog(div('UMAP reduction is calculated! Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        seurat_object <<- RunUMAP(seurat_object, dims = 1:as.numeric(input$umapPCs), seed.use = as.numeric(input$umapSeed), n.components = as.numeric(input$umapOutComponents), reduction = "pca")
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " TRAJECTORY ANALYSIS", " GENE REGULATORY NETWORK\nANALYSIS"))
        updateUmapTypeChoices("umap")
        updateReduction2("umap")
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the UMAP calculation.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----UMAP_end_Mins_taken-----", time.taken))

      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "umapRunUmap")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunTsne")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunDFM")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunPhate")
    })
  })
  
  observeEvent(input$umapRunTsne, {
    session$sendCustomMessage("handler_disableAllButtons", "umapRunUmap")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunTsne")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunDFM")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunPhate")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (!"pca" %in% names(seurat_object)) session$sendCustomMessage("handler_alert", "Please, first execute PRINCIPAL COMPONENT ANALYSIS.")
      else {
        showModal(modalDialog(div('tSNE reduction is calculated! Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        if(as.numeric(input$umapOutComponents) == 2)
          seurat_object <<- RunTSNE(seurat_object, dims = 1:as.numeric(input$umapPCs), seed.use = as.numeric(input$umapSeed), dim.embed = 2, reduction = "pca", verbose = T)
        else
          seurat_object <<- RunTSNE(seurat_object, dims = 1:as.numeric(input$umapPCs), seed.use = as.numeric(input$umapSeed), dim.embed = 3, reduction = "pca", verbose = T)
        
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " TRAJECTORY ANALYSIS"))
        updateUmapTypeChoices("tsne")
        updateReduction2("tsne")
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the tSNE calculation.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "umapRunUmap")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunTsne")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunDFM")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunPhate")
    })
  })
  
  observeEvent(input$umapRunDFM, {
    session$sendCustomMessage("handler_disableAllButtons", "umapRunUmap")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunTsne")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunDFM")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunPhate")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (!"pca" %in% names(seurat_object)) session$sendCustomMessage("handler_alert", "Please, first execute PRINCIPAL COMPONENT ANALYSIS.")
      else {
        showModal(modalDialog(div('Diffusion map reduction is calculated! Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        #prepare input
        dfm_in <- as.data.frame(seurat_object@reductions[['pca']]@cell.embeddings)
        dfm_in$Cell_id <- rownames(dfm_in)
        dfm_in <- dfm_in[, c(51, 1:as.numeric(input$umapPCs))]
    
        #run DFM default
        dfm <- DiffusionMap(dfm_in, n_eigs = as.numeric(input$umapOutComponents))
        dfm_out <- dfm@eigenvectors
        colnames(dfm_out) <- gsub("DC", "DC_", colnames(dfm_out))
    
        #add new reduction in seurat_object
        seurat_object[["dfm"]] <<- CreateDimReducObject(embeddings = dfm_out, key = "DC_", assay = DefaultAssay(seurat_object), global = T)
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " TRAJECTORY ANALYSIS"))
        updateUmapTypeChoices("dfm")
        updateReduction2("dfm")
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the Diffusion Map calculation.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "umapRunUmap")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunTsne")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunDFM")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunPhate")
    })
  })
  
  observeEvent(input$umapRunPhate, {
    start.time <- Sys.time()
    print(paste0("!!!!-----Phate_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "umapRunUmap")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunTsne")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunDFM")
    session$sendCustomMessage("handler_disableAllButtons", "umapRunPhate")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (!"pca" %in% names(seurat_object)) session$sendCustomMessage("handler_alert", "Please, first execute PRINCIPAL COMPONENT ANALYSIS.")
      else {
        showModal(modalDialog(div('PHATE reduction is calculated! This is slow procedure, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        seurat_object <<- RunPHATE(seurat_object, dims = 1:as.numeric(input$umapPCs), n.components = as.numeric(input$umapOutComponents), reduction = "pca", seed.use = as.numeric(input$umapSeed))
        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " TRAJECTORY ANALYSIS"))
        updateUmapTypeChoices("phate")
        updateReduction2("phate")
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the Phate calculation.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Phate_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "umapRunUmap")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunTsne")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunDFM")
      session$sendCustomMessage("handler_enableAllButtons", "umapRunPhate")
    })
  })
  
  observeEvent(input$umapConfirm, { 
    if(input$umapType != "-")
      updateReduction()
  })
  
  #------------------DEA tab-----------------------------------------------
  observeEvent(input$findMarkersConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----DEA_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "findMarkersConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first and then proceed with Marker Gene identification.")
      else {
        showModal(modalDialog(div('Analysis in Progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        markers_logFCBase <<- input$findMarkersLogBase
        
        if(markers_logFCBase == "avg_logFC")
        {
          seurat_object@misc$markers <<- FindAllMarkers(seurat_object, test.use = input$findMarkersTest, min.pct = as.numeric(input$findMarkersMinPct), logfc.threshold = as.numeric(input$findMarkersLogFC), 
                                                        return.thresh = as.numeric(input$findMarkersPval), base = exp(1))
        }
        else
        {
          seurat_object@misc$markers <<- FindAllMarkers(seurat_object, test.use = input$findMarkersTest, min.pct = as.numeric(input$findMarkersMinPct), logfc.threshold = as.numeric(input$findMarkersLogFC), 
                                                        return.thresh = as.numeric(input$findMarkersPval), base = 2)
        }
        updateInputLRclusters()
        
        output$findMarkersTable <- renderDataTable(
          {
            if (!is.null(seurat_object@misc$markers))
            {
              seurat_object@misc$markers <<- seurat_object@misc$markers[, c(7, 6, 1:5)]
              output$findMarkersTable <- renderDataTable(seurat_object@misc$markers, options = list(pageLength = 20), filter = 'top', rownames = FALSE)  
              export_markerGenes_RNA <<- seurat_object@misc$markers
            }
          }
        )

        session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " FUNCTIONAL/MOTIF\nENRICHMENT ANALYSIS", " CLUSTERS' ANNOTATION"))
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the DE Analysis")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----DEA_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "findMarkersConfirm")
    })
  })

  observeEvent(input$findMarkersTop10HeatmapConfirm, {
    #DEA output rendering
    session$sendCustomMessage("handler_disableAllButtons", "findMarkersTop10HeatmapConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else{
        showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("findMarkersHeatmap_loader")
        
        top10 <- seurat_object@misc$markers %>% group_by(cluster) %>% top_n(n = 10, wt = eval(parse(text=markers_logFCBase))) 
        
        downsampled <- seurat_object
        if(ncol(seurat_object) > 1500)
        {
          set.seed(9)
          downsampled <- subset(seurat_object, cells = sample(Cells(seurat_object), 1500))  
        }
        
        scaled_tabe <- as.data.frame(downsampled@assays$RNA@scale.data)
        scaled_tabe$gene <- rownames(scaled_tabe)
        scaled_tabe_order <- as.data.frame(top10$gene)
        colnames(scaled_tabe_order)[1] <- "gene"
        scaled_tabe_final <- left_join(scaled_tabe_order, scaled_tabe)
        scaled_tabe_final <- na.omit(scaled_tabe_final)
        tableCl <- downsampled@meta.data[, ]
        tableCl$Cell_id <- rownames(tableCl)
        tableCl <- tableCl[, c('Cell_id', 'seurat_clusters')]
        tableCl <- tableCl[order(tableCl$seurat_clusters), ]
        
        clip<-function(x, min=-2, max=2) {
          x[x<min]<-min; 
          x[x>max]<-max; 
          x
        }

        final_mat <- scaled_tabe_final[, -1]
        final_mat <- final_mat[, tableCl$Cell_id]
        
        cols <- colorRampPalette(brewer.pal(12, "Paired"))(length(unique(tableCl$seurat_clusters)))
        names(cols) <- unique(tableCl$seurat_clusters)
        
        output$findMarkersHeatmap <- renderPlotly({
          heatmaply(clip(final_mat), Rowv = F, Colv = F, colors = rev(RdBu(256)), showticklabels = c(F, T), labRow  = scaled_tabe_final$gene, 
                    col_side_colors = tableCl$seurat_clusters, col_side_palette =  cols)
        })
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the visualization of the plot.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "findMarkersTop10HeatmapConfirm")
    })  
      
  })
  
  observeEvent(input$findMarkersTop10DotplotConfirm, {
    session$sendCustomMessage("handler_disableAllButtons", "findMarkersTop10DotplotConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else{
        showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("findMarkersDotplot_loader")
        output$findMarkersDotplot <- renderPlotly(
          {
            top10 <- seurat_object@misc$markers %>% group_by(cluster) %>% top_n(n = 10, wt = eval(parse(text=markers_logFCBase))) 
            p <- DotPlot(seurat_object, features = rev(unique(top10$gene)), dot.scale = 6, cols = c("grey", "red")) + RotatedAxis() + labs(fill="Average\nexpression")
            plotly::ggplotly(p)
          }
        )
        
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the visualization of the plot.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "findMarkersTop10DotplotConfirm")
    })
  })
  
  observeEvent(input$findMarkersSignatureAdd, {
    start.time <- Sys.time()
    print(paste0("!!!!-----Signature_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "findMarkersSignatureAdd")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else{
        showModal(modalDialog(div('Gene signature addition in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        markers <- list()
        varTextarea <- input$findMarkersSignatureMembers
        markers2 <- as.vector(unlist(strsplit(varTextarea, "\\n"))) 
        print(markers2)
        print(class(markers2))
        print(length(markers2))
        markers[[1]] <- unlist(strsplit(varTextarea, "\\n"))
        sig_name <- input$findMarkersSignatureName
        print(sig_name)
        names(markers)[1] <- sig_name
        
        genesFound <- markers2[(which(markers2 %in% rownames(seurat_object)))]
        print(genesFound)
        
        if(length(genesFound) == 0)
        {
          print("inside all zero")
          session$sendCustomMessage("handler_alert", "None of the selected genes were found in the dataset.The signature scores could not be calculated.")
        }
        else
        {
          seurat_object <<- AddModuleScore_UCell(seurat_object, features = markers)
          updateGeneSearchFP()
          updateMetadata()
          if(length(genesFound) != length(markers2))
          {
            print("inside some not found")
            session$sendCustomMessage("handler_alert", "Some of the selected genes were not found in the dataset and are excluded from the signature's calculation.")
          }
        }
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "The signature could not be added.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Signature_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "findMarkersSignatureAdd")
    })
  })
    
observeEvent(input$findMarkersFPConfirm, {
  session$sendCustomMessage("handler_disableAllButtons", "findMarkersFPConfirm")
  tryCatch({
    if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else{
      showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      shinyjs::show("findMarkersFeaturePlot_loader")
      
      if(input$findMarkersReductionType != "-")
      {
        geneS <- input$findMarkersGeneSelect
        
        label_x <- ""
        label_y <- ""
        show_label <- as.logical(input$findMarkersLabels)
        order_exp <- as.logical(input$findMarkersOrder)
        minq <- paste0("q", input$findMarkersMinCutoff)
        maxq <- paste0("q", input$findMarkersMaxCutoff)
        
        if(input$findMarkersReductionType == "umap")
        {
          label_x <- "UMAP_1"
          label_y <- "UMAP_2"
        }
        else if(input$findMarkersReductionType == "tsne")
        {
          label_x <- "tSNE_1"
          label_y <- "tSNE_2"
        }
        else if(input$findMarkersReductionType == "dfm")
        {
          label_x <- "DC_1"
          label_y <- "DC_2"
        }
        else if(input$findMarkersReductionType == "pca")
        {
          label_x <- "PC_1"
          label_y <- "PC_2"
        }
        else if(input$findMarkersReductionType == "phate")
        {
          label_x <- "PHATE_1"
          label_y <- "PHATE_2"
        }
        
        plot_temp <- FeaturePlot(seurat_object, features = geneS, pt.size = 1.5, label = show_label, label.size = 5, cols = c("lightgrey", "red"), 
                            order = order_exp, reduction = input$findMarkersReductionType, max.cutoff = maxq, min.cutoff = minq)
        
        plot <- FeaturePlot(seurat_object, features = geneS, pt.size = 1.5, label = show_label, label.size = 5, cols = c("lightgrey", "red"), 
                            order = order_exp, reduction = input$findMarkersReductionType, max.cutoff = maxq, min.cutoff = minq, raster = F) + 
          xlim(min(plot_temp$data[label_x]), max(plot_temp$data[label_x])) + ylim(min(plot_temp$data[label_y]), max(plot_temp$data[label_y])) +
          theme_bw() +
          theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                axis.title.y = element_text(face = "bold", color = "black", size = 25),
                axis.title.x = element_text(face = "bold", color = "black", size = 25),
                legend.text = element_text(face = "bold", color = "black", size = 9),
                legend.title = element_text(face = "bold", color = "black", size = 9),
                legend.position="right",
                title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
          labs(x=label_x, y=label_y, title = geneS, color="")
        gp <- plotly::ggplotly(plot, tooltip = c("x", "y", geneS))
        output$findMarkersFeaturePlot <- renderPlotly({
          gp
        })
      }
      
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was a problem with the generation of the plot.")
  }, finally = {
    removeModal()
    Sys.sleep(1)
    session$sendCustomMessage("handler_enableAllButtons", "findMarkersFPConfirm")
  })
})

observeEvent(input$findMarkersFeaturePairConfirm, {
  session$sendCustomMessage("handler_disableAllButtons", "findMarkersFeaturePairConfirm")
  tryCatch({
    if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else{
      showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      shinyjs::show("findMarkersFPfeature1_loader")
      shinyjs::show("findMarkersFPfeature2_loader")
      shinyjs::show("findMarkersFPfeature1_2_loader")
      shinyjs::show("findMarkersFPcolorbox_loader")
      
      updateFeaturePair()
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was a problem with the generation of the plot.")
  }, finally = {
    removeModal()
    Sys.sleep(1)
    session$sendCustomMessage("handler_enableAllButtons", "findMarkersFeaturePairConfirm")
  })
})

observeEvent(input$findMarkersViolinConfirm, {
  session$sendCustomMessage("handler_disableAllButtons", "findMarkersViolinConfirm")
  tryCatch({
    if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else{
      showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      shinyjs::show("findMarkersViolinPlot_loader")
      if (!identical(input$findMarkersGeneSelect2, NULL)){
        output$findMarkersViolinPlot <- renderPlotly(
          {
            
            geneS <- input$findMarkersGeneSelect2
            print(geneS)
            plot <- VlnPlot(seurat_object, features = geneS, pt.size = 0, raster = F,
                            cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, 'seurat_clusters'])))) +
              theme_bw() +
              theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                    axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                    axis.title.y = element_text(face = "bold", color = "black", size = 25),
                    axis.title.x = element_text(face = "bold", color = "black", size = 25),
                    legend.text = element_text(face = "bold", color = "black", size = 9),
                    legend.title = element_text(face = "bold", color = "black", size = 9),
                    title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
              labs(x="Cluster", y="", title = geneS, fill="Cluster")
            gp <- plotly::ggplotly(plot, tooltip = c("x", "y", geneS))
            gp
          }
        )
      }
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was a problem with the generation of the plot.")
  }, finally = {
    removeModal()
    Sys.sleep(1)
    session$sendCustomMessage("handler_enableAllButtons", "findMarkersViolinConfirm")
  })
})

observeEvent(input$findMarkersVolcanoConfirm, {
  session$sendCustomMessage("handler_disableAllButtons", "findMarkersVolcanoConfirm")
  tryCatch({
    if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else{
      showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      shinyjs::show("findMarkersVolcanoPlot_loader")
      diff_exp_genes <- seurat_object@misc$markers
      cluster_degs <- diff_exp_genes[which(diff_exp_genes$cluster == input$findMarkersClusterSelect), ]
      cluster_degs$status <- "Down regulated"
      for(i in 1:length(cluster_degs$gene))
      {
        if(cluster_degs[i, markers_logFCBase] > 0) 
        {
          cluster_degs$status[i] <- "Up regulated"
        }
      }
      cluster_degs$log10Pval <- -log10(cluster_degs$p_val)
      
      p <- ggplot(data=cluster_degs, aes_string(x=markers_logFCBase, y="log10Pval", fill="status", label="gene", color="status")) + 
        geom_point(size=1, shape=16)+
        scale_fill_manual(values = c("cyan3", "orange"))+
        scale_color_manual(values = c("cyan3", "orange"))+
        scale_size()+
        theme_bw() +
        theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
              axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
              axis.title.y = element_text(face = "bold", color = "black", size = 25),
              axis.title.x = element_text(face = "bold", color = "black", size = 25),
              legend.text = element_text(face = "bold", color = "black", size = 9),
              legend.title = element_text(face = "bold", color = "black", size = 9),
              legend.position="right",
              title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
        labs(x=paste0("", markers_logFCBase), y="-log10(Pvalue)", fill="Color", color="")
      
      output$findMarkersVolcanoPlot <- renderPlotly(
        {
          plotly::ggplotly(p, tooltip = c("x", "y", "label"))
        })
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was a problem with the generation of the plot.")
  }, finally = {
    removeModal()
    Sys.sleep(1)
    session$sendCustomMessage("handler_enableAllButtons", "findMarkersVolcanoConfirm")
  })
})

#ATAC   
observeEvent(input$findMarkersConfirmATAC, {
  start.time <- Sys.time()
  print(paste0("!!!!-----MarkerGenes_start-----", start.time))
  
  session$sendCustomMessage("handler_disableAllButtons", "findMarkersConfirmATAC")
  tryCatch({
    if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else {
      showModal(modalDialog(div('Marker genes identification in progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      shinyjs::show("findMarkersGenesHeatmapATAC_loader")
      shinyjs::show("findMarkersGenesATACTable_loader")
      markers_cluster <- getMarkerFeatures(
        ArchRProj = proj_default,
        useMatrix = "GeneScoreMatrix",
        groupBy = input$findMarkersGroupByATAC, #"predictedGroup_Co", #"Clusters",
        bias = c("TSSEnrichment", "log10(nFrags)"),
        testMethod = input$findMarkersTestATAC, 
        logFile =paste0(user_dir, "/Marker_Genes_file.log")
      )
      
      markers_clusterList <- getMarkers(markers_cluster, cutOff = paste0("FDR <= ",input$findMarkersFDRATAC," & Log2FC >= ",input$findMarkersLogFCATAC))
      markers_clusters_all <- as.data.frame(unlist(markers_clusterList))
      markers_clusters_all$Clusters <- rownames(markers_clusters_all)
      markers_clusters_all$Clusters <- gsub(pattern = "[.][0-9]+", replacement = "", x = markers_clusters_all$Clusters)
      rownames(markers_clusters_all) <- NULL
      
      output$findMarkersGenesTableATAC <- renderDataTable(expr = markers_clusters_all, filter = 'top', rownames = FALSE, options = list(pageLength = 10))
      export_markerGenes_ATAC <<- markers_clusters_all
      
      #heatmap top10
      markers_clusterList10 <- unlist(getMarkers(markers_cluster, cutOff = paste0("FDR <= ",input$findMarkersFDRATAC," & Log2FC >= ",input$findMarkersLogFCATAC), n = 10))
      heatmap_matrix <- plotMarkerHeatmap(
        seMarker = markers_cluster,
        cutOff = paste0("FDR <= ",input$findMarkersFDRATAC," & Log2FC >= ",input$findMarkersLogFCATAC),
        labelMarkers = NULL,
        transpose = TRUE,
        returnMatrix = TRUE,
        logFile =paste0(user_dir, "/Marker_Heatmap_file.log")
      )
      output$findMarkersGenesHeatmapATAC <- renderPlotly(expr = heatmaply(heatmap_matrix[, markers_clusterList10$name], colors=paletteContinuous(set = "blueYellow", n = 100))) 
      saveArchRProject(proj_default)
      session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " FUNCTIONAL/MOTIF\nENRICHMENT ANALYSIS", " TRACKS"))
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was an error with the the detection of marker genes.")
  }, finally = {
    end.time <- Sys.time()
    time.taken <- as.numeric (end.time - start.time, units = "mins")
    print(paste0("!!!!-----MarkerGenes_end_Mins_taken-----", time.taken))
    
    removeModal()
    Sys.sleep(1)
    session$sendCustomMessage("handler_enableAllButtons", "findMarkersConfirmATAC")
  })	
})

observeEvent(input$findMarkersPeaksConfirmATAC, {
  start.time <- Sys.time()
  print(paste0("!!!!-----MarkerPeaks_start-----", start.time))
  
  
  session$sendCustomMessage("handler_disableAllButtons", "findMarkersPeaksConfirmATAC")
  tryCatch({
    if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else {
      showModal(modalDialog(div('Marker peaks identification in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      updateTabsetPanel(session, inputId = "ATAC_markers_tabs", selected = "Marker peaks (ATAC)")
      shinyjs::show("findMarkersPeaksHeatmapATAC_loader")
      shinyjs::show("findMarkersPeaksATACTable_loader")
      source("Peaks_ArchR_windows.R")
      
      #addArchRThreads(threads = as.numeric(input$upload10xATACThreads))
      
      # if(.Platform$OS.type == "windows")
      # {
      #   pathToMacs2_a<-system("bash -c 'find /home -name macs2'", intern = TRUE)
      #   
      #   proj_default <<- addGroupCoverages(ArchRProj = proj_default, groupBy = "Clusters",force = TRUE, 
      #                                      logFile =paste0(user_dir, "/Group_Coverages_file.log"))
      #   proj_default <<- addReproduciblePeakSet_win(
      #     ArchRProj = proj_default,
      #     groupBy = "Clusters",
      #     pathToMacs2 = pathToMacs2_a, force = T,
      #     logFile =paste0(user_dir, "/Rep_Peakset_file.log")
      #   )
      # }
      # else
      # {
      #   proj_default <<- addGroupCoverages(ArchRProj = proj_default, groupBy = "Clusters", force = T, 
      #                                      logFile =paste0(user_dir, "/Group_Coverages_file.log"))
      #   
      #   proj_default <<- addReproduciblePeakSet(
      #     ArchRProj = proj_default,
      #     groupBy = "Clusters",
      #     pathToMacs2 = input$pathToMacs2, #"/home/user/anaconda3/bin/macs2",
      #     force = T, 
      #     logFile =paste0(user_dir, "/Rep_Peakset_file.log")
      #   )
      #   print("Peakset finished L")
      # }
      
      #----custom peaks from BED
      tile_matrix<-getMatrixFromProject(proj_default,"TileMatrix",binarize = TRUE, 
                                        logFile =paste0(user_dir, "/get_Matrix_From_Project_file.log"))
      valid_chromosomes<-unique(rowData(tile_matrix)$seqnames)
      tile_matrix<-c()
      file.copy(from = input$findMarkersPeaksCustomPeaks$datapath, to = paste0(user_dir, "/bedFile.bed"), overwrite = TRUE)#load bed
      proj_default_peaks<-read.delim(paste0(user_dir, "/bedFile.bed") ,header=F)
      proj_default_peaks<-proj_default_peaks[,c(1:3)]
      if(substr(proj_default_peaks[1,1], 1, 3) != "chr" || !is.numeric(proj_default_peaks[1,2]) || !is.numeric(proj_default_peaks[1,3])){
        stop("Invalid bed file!! Check UCSC specifications in order to provide a valid bed file. Chromosome names should start with chr and start and end positions 
             should be numeric integers. No header should be provided.")
      }
      temp<-proj_default_peaks[,3]-proj_default_peaks[,2]
      if(all(temp<0)){
        stop("Invalid bed file!! Start coordinates should be smaller than end coordinates.")
      }
      colnames(proj_default_peaks)[1]<-"chr"
      colnames(proj_default_peaks)[2]<-"start"
      colnames(proj_default_peaks)[3]<-"end"
      proj_default_peaks$V4<-paste0(proj_default_peaks$chr,":",proj_default_peaks$start,"-",proj_default_peaks$end)
      colnames(proj_default_peaks)[4]<-"id"
      proj_default_peaks_filtered<-proj_default_peaks[which(proj_default_peaks$chr %in% valid_chromosomes),]
      proj_default_peaks_gr<-makeGRangesFromDataFrame(proj_default_peaks_filtered,keep.extra.columns=TRUE)
      proj_default <<- addPeakSet(proj_default, peakSet = proj_default_peaks_gr, force=TRUE)
      proj_default <<- addPeakMatrix(proj_default, 
                                     logFile =paste0(user_dir, "/Peak_Matrix_file.log"))

      print(getPeakSet(proj_default))
      saveArchRProject(proj_default)
      
      markersPeaks <- getMarkerFeatures(
        ArchRProj = proj_default,
        useMatrix = "PeakMatrix",
        groupBy = input$findMarkersPeaksGroupByATAC,#"predictedGroup_Co", #"Clusters",
        bias = c("TSSEnrichment", "log10(nFrags)"),
        testMethod = input$findMarkersPeaksTestATAC, 
        logFile =paste0(user_dir, "/Marker_Peaks_file.log")
      )
      
      markersPeaks_clusterList <- getMarkers(markersPeaks, cutOff = paste0("FDR <= ",input$findMarkersPeaksFDRATAC," & Log2FC >= ",input$findMarkersPeaksLogFCATAC))
      markersPeaks_clusters_all <- as.data.frame(unlist(markersPeaks_clusterList))
      markersPeaks_clusters_all$Clusters <- rownames(markersPeaks_clusters_all)
      markersPeaks_clusters_all$Clusters <- gsub(pattern = "[.][0-9]+", replacement = "", x = markersPeaks_clusters_all$Clusters)
      rownames(markersPeaks_clusters_all) <- NULL
      output$findMarkersPeaksTableATAC <- renderDataTable(expr = markersPeaks_clusters_all, filter = 'top', rownames = FALSE, options = list(pageLength = 10))
      export_markerPeaks_ATAC <<- markersPeaks_clusters_all
      
      markerList <- getMarkers(markersPeaks, cutOff = paste0("FDR <= ",input$findMarkersPeaksFDRATAC," & Log2FC >= ",input$findMarkersPeaksLogFCATAC), n = 10)
      top10peaks <- as.data.frame(unlist(markerList))
      top10peaks$names <- paste0(top10peaks$seqnames, ":", top10peaks$start,"-",top10peaks$end)
      
      markerListheatmapPeaks <- plotMarkerHeatmap(
        seMarker = markersPeaks, 
        cutOff = paste0("FDR <= ",input$findMarkersPeaksFDRATAC," & Log2FC >= ",input$findMarkersPeaksLogFCATAC), returnMatrix = T,
        transpose = TRUE, 
        logFile =paste0(user_dir, "/Marker_Peaks_heatmap_file.log")
      )
      to_plot <- markerListheatmapPeaks[, top10peaks$names]
      
      output$findMarkersPeaksHeatmapATAC <- renderPlotly(expr = heatmaply(to_plot, colors=paletteContinuous(set = "solarExtra", n = 100))) 
      addArchRThreads(threads = as.numeric(input$upload10xATACThreads))
      session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " FUNCTIONAL/MOTIF\nENRICHMENT ANALYSIS", " TRACKS"))
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was an error with the the detection of marker peaks.")
  }, finally = {
    end.time <- Sys.time()
    time.taken <- as.numeric (end.time - start.time, units = "mins")
    print(paste0("!!!!-----MarkerPeaks_end_Mins_taken-----", time.taken))
    
    removeModal()
    Sys.sleep(1)
    session$sendCustomMessage("handler_enableAllButtons", "findMarkersPeaksConfirmATAC")
  })	
})

observeEvent(input$findMarkersFPConfirmATAC, {
  session$sendCustomMessage("handler_disableAllButtons", "findMarkersFPConfirmATAC")
  tryCatch({
    if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else {
      showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      shinyjs::show("findMarkersFeaturePlotATAC_loader")
      
      p <- plotEmbedding(
        ArchRProj = proj_default, 
        colorBy = "GeneScoreMatrix", 
        name = input$findMarkersGeneSelectATAC, 
        embedding = input$findMarkersReductionTypeATAC,
        imputeWeights = getImputeWeights(proj_default),
        logFile =paste0(user_dir, "/Plot_embedding_file.log")
      )
    
      output$findMarkersFeaturePlotATAC <- renderPlot(expr = p)
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was an error with the the detection of marker genes.")
  }, finally = {
    removeModal()
    Sys.sleep(1)
    session$sendCustomMessage("handler_enableAllButtons", "findMarkersFPConfirmATAC")
  })
})

output$findMarkersRNAExport <- downloadHandler(
  filename = function() { 
    paste("markerGenesTableRNA-", Sys.Date(), ".txt", sep="")
  },
  content = function(file) {
    write.table(export_markerGenes_RNA, file, sep = "\t", quote = F, row.names = F)
  })

output$findMarkersGenesATACExport <- downloadHandler(
  filename = function() { 
    paste("markerGenesTableATAC-", Sys.Date(), ".txt", sep="")
  },
  content = function(file) {
    write.table(export_markerGenes_ATAC, file, sep = "\t", quote = F, row.names = F)
  })

output$findMarkersPeaksATACExport <- downloadHandler(
  filename = function() { 
    paste("markerPeaksTableATAC-", Sys.Date(), ".txt", sep="")
  },
  content = function(file) {
    write.table(export_markerPeaks_ATAC, file, sep = "\t", quote = F, row.names = F)
  })
  
#------------------Doublet detection------------------------------------------
observeEvent(input$doubletsConfirm, {
  start.time <- Sys.time()
  print(paste0("!!!!-----Doublets_start-----", start.time))
  
  session$sendCustomMessage("handler_disableAllButtons", "doubletsConfirm")
  tryCatch({
    if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else
    {
      showModal(modalDialog(div('Doublet detection in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      
      #delete previous results of doublet detection
      for_delete <- grep("DF.classifications_|pANN_", colnames(seurat_object@meta.data))
      if(length(for_delete) != 0)
      {
        seurat_object@meta.data <<- seurat_object@meta.data[, -for_delete]
      }
      
      homotypic.prop.seurat_object <- modelHomotypic(seurat_object@meta.data$seurat_clusters) 
      nExp_poi.seurat_object <- round( as.numeric(input$doubletsNExp) * nrow(seurat_object@meta.data) )  ## Assuming x% doublet formation rate
      nExp_poi.adj.seurat_object <- round(nExp_poi.seurat_object * (1-homotypic.prop.seurat_object) )
      pK_choose <- 0.09
      if(input$doubletsPKRadio == "auto")
      {
        #auto pK estimation
        sweep.res.list_seurat_object <- paramSweep_v3(seurat_object, PCs = 1:as.numeric(input$doubletsPCs), sct = FALSE)
        sweep.stats_seurat_object <- summarizeSweep(sweep.res.list_seurat_object, GT = FALSE)
        bcmvn_seurat_object <- find.pK(sweep.stats_seurat_object)
        pK=as.numeric(as.character(bcmvn_seurat_object$pK))
        BCmetric=bcmvn_seurat_object$BCmetric
        pK_choose = pK[which(BCmetric %in% max(BCmetric))]
      }
      else
      {
        pK_choose <- input$doubletsPK
      }
      
      seurat_object <<- doubletFinder_v3(seurat_object, PCs = 1:as.numeric(input$doubletsPCs), pN = as.numeric(input$doubletsPN), pK = as.numeric(pK_choose), nExp = nExp_poi.seurat_object, reuse.pANN = FALSE, sct = FALSE)
      seurat_object <<- doubletFinder_v3(seurat_object, PCs = 1:as.numeric(input$doubletsPCs), pN = as.numeric(input$doubletsPN), pK = as.numeric(pK_choose), nExp = nExp_poi.adj.seurat_object, sct = FALSE)
      
      #factor the two last columns
      seurat_object@meta.data[, ncol(seurat_object@meta.data)] <<- factor(seurat_object@meta.data[, ncol(seurat_object@meta.data)])
      seurat_object@meta.data[, ncol(seurat_object@meta.data)-2] <<- factor(seurat_object@meta.data[, ncol(seurat_object@meta.data)-2])
      
      #output 
      column_total_name <- paste0("DF.classifications_", as.numeric(input$doubletsPN), "_", pK_choose, "_", nExp_poi.seurat_object)
      column_not_homotypic_name <- paste0("DF.classifications_", as.numeric(input$doubletsPN), "_", pK_choose, "_", nExp_poi.adj.seurat_object) 
      df_total <- (table(seurat_object@meta.data[, column_total_name]))
      df_not_homotypic <- (table(seurat_object@meta.data[, column_not_homotypic_name]))
      print(head(seurat_object@meta.data))
      print(class(seurat_object@meta.data[, column_not_homotypic_name]))
      
      updateMetadata()
      updateSelInpColor()
      
      output$doubletsInfo <- renderPrint(
        {
          cat(paste0("Total number of cells: ", nrow(seurat_object@meta.data), "\nTotal number of doublets: ", df_total["Doublet"], "\nHomotypic doublets:", df_total["Doublet"]-df_not_homotypic["Doublet"],
                     "\nNon-homotypic doublets:", df_not_homotypic["Doublet"]))
        })
      
      output$doubletPCAplot <- renderPlotly(
        {
          if(!is.null(seurat_object) & input$doubletsReduction != "-")
          {
            meta <- seurat_object@meta.data
            meta$Cell_id <- rownames(meta)
            label_x <- "PC_1"
            label_y <- "PC_2"
            selected_Reduction <- input$doubletsReduction
            
            if(selected_Reduction == "umap")
            {
              label_x <- "UMAP_1"
              label_y <- "UMAP_2"
            }
            else if(selected_Reduction == "tsne")
            {
              label_x <- "tSNE_1"
              label_y <- "tSNE_2"
            }
            else if(selected_Reduction == "dfm")
            {
              label_x <- "DC_1"
              label_y <- "DC_2"
            }
            else if(selected_Reduction == "phate")
            {
              label_x <- "PHATE_1"
              label_y <- "PHATE_2"
            }
            
            plot1 <- DimPlot(seurat_object, reduction = selected_Reduction)
            plot1_data <- plot1$data
            plot1_data$Cell_id <- rownames(plot1_data)
            plot1_data <- left_join(plot1_data, meta)
            
            p <- ggplot(plot1_data, aes_string(x=label_x, y=label_y, color=column_total_name))+
              geom_point() +
              theme_bw() +
              labs(x=label_x, y=label_y)+
              theme(legend.position = "right")
            
            gp <- plotly::ggplotly(p, tooltip = c("x", "y"))
            print(gp)
          }
        }
      )
    }
  }, error = function(e) 
  {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was an error during doublets' detection.")
  }, finally = 
    {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Doublets_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "doubletsConfirm")
    }
  )
})

observeEvent(input$doubletsRemove, { 
  session$sendCustomMessage("handler_disableAllButtons", "doubletsRemove")
  tryCatch({
    if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else
    {
      ident_to_remove <- "Doublet"
      if(input$doubletsCellRemoval == "all")
      {
        Idents(seurat_object) <<- seurat_object@meta.data[, grep(x = colnames(seurat_object@meta.data), pattern = "DF.classifications")[1]]
      }
      else
      {
        Idents(seurat_object) <<- seurat_object@meta.data[, grep(x = colnames(seurat_object@meta.data), pattern = "DF.classifications")[2]]
      }
  
      seurat_object <<- subset(seurat_object, idents = ident_to_remove, invert = TRUE)
      Idents(seurat_object) <<- seurat_object$seurat_clusters
      updateInputLRclusters()
      updateMetadata()
      shinyalert("Doublets were removed successfully.", type = "success")
    }
    }, error = function(e) 
      {
        print(paste("Error :  ", e))
        shinyalert("There was an error during doublets removal. Maybe doublet detection is not completed?", type = "error")
      }, finally = 
      {
        session$sendCustomMessage("handler_enableAllButtons", "doubletsRemove")
      }
  )
})

#ATAC 
observeEvent(input$doubletsATACConfirm, {
  start.time <- Sys.time()
  print(paste0("!!!!-----Doublets_start-----", start.time))
  
  session$sendCustomMessage("handler_disableAllButtons", "doubletsATACConfirm")
  tryCatch({
    if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else
    {
      showModal(modalDialog(div('Doublet detection in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      shinyjs::show("doubletATAC_loader3")
      shinyjs::show("doubletATAC_loader4")
      
      proj_default <<- addDoubletScores(
        input = proj_default,
        k = as.numeric(input$doubletsATACk), #Refers to how many cells near a "pseudo-doublet" to count.
        knnMethod = "umap", #Refers to the embedding to use for nearest neighbor search with doublet projection.
        LSIMethod = as.numeric(input$doubletsATACLSI),
        logFile =paste0(user_dir, "/Doublet_Scores_file.log")
      )
      
      updateMetadataATAC()
      #plot1 Doublet Score
      pal <- c("grey", "#FB8861FF", "#B63679FF", "#51127CFF", "#000004FF") #grey_magma
      df<-proj_default@embeddings$umap$df
      df$score<-proj_default$DoubletScore
      df$enrichment<-proj_default$DoubletEnrichment
      
      xlim <- range(df[,1]) %>% extendrange(f = 0.05)
      ylim <- range(df[,2]) %>% extendrange(f = 0.05)
      
      pscore <- ggPoint(
        x = df[,1],
        y = df[,2],
        color = ArchR:::.quantileCut(df$score, 0, 0.95),
        xlim = xlim,
        ylim = ylim,
        discrete = FALSE,
        size = 5,
        xlab = "UMAP Dimension 1",
        ylab = "UMAP Dimension 2",
        pal = pal,
        title = "Doublet Scores -log10(P-adj.)",
        colorTitle = "Doublet Scores -log10(P-adj.)",
        rastr = TRUE,
        baseSize = 20
      ) + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), 
                axis.text.y = element_blank(), axis.ticks.y = element_blank())
      
      #plot2 Doublet enrichment
      penrich <- ggPoint(
        x = df[,1],
        y = df[,2],
        color = ArchR:::.quantileCut(df$enrichment, 0, 0.95),
        xlim = xlim,
        ylim = ylim,
        discrete = FALSE,
        size = 5,
        xlab = "UMAP Dimension 1",
        ylab = "UMAP Dimension 2",
        pal = pal,
        title = "Simulated Doublet Enrichment over Expectation",
        colorTitle = "Doublet Enrichment",
        rastr = TRUE,
        baseSize = 20
      ) + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), 
                axis.text.y = element_blank(), axis.ticks.y = element_blank())
      
      output$doubletsScoreATAC <- renderPlot( {pscore} )
      
      output$doubletEnrichmentATAC <- renderPlot( {penrich} )
    }
  }, error = function(e) 
  {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was an error with drawing the resutls.")
  }, finally = 
    {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Doublets_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "doubletsATACConfirm")
    }
  )
})

observeEvent(input$doubletsATACDelete, {
  session$sendCustomMessage("handler_disableAllButtons", "doubletsATACDelete")
  tryCatch({
    if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else
    {
      proj_default <<- filterDoublets(proj_default, filterRatio = as.numeric(input$doubletsATACfilterRatio))
      updateMetadataATAC()
      shinyalert("Doublets were removed successfully.", type = "success")
    }
  }, error = function(e) 
  {
    print(paste("Error :  ", e))
    shinyalert("There was an error during doublets removal.", type = "error")
  }, finally = 
    {
      session$sendCustomMessage("handler_enableAllButtons", "doubletsATACDelete")
    }
  )
})

#------------------Cell cycle tab---------------------------------------------
  observeEvent(input$cellCycleRun, { 
    start.time <- Sys.time()
    print(paste0("!!!!-----Cell_cycle_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "cellCycleRun")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (!"pca" %in% names(seurat_object)) session$sendCustomMessage("handler_alert", "Please, first execute PRINCIPAL COMPONENT ANALYSIS.") 
      else {
        showModal(modalDialog(div('Cell cycle phase analysis in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("cellCyclePCA_loader")
        shinyjs::show("cellCycleBarplot_loader")
        
        if(organism == "mouse")
        {
          s.genes_ver2 <- readRDS("cell_cycle_phase_helper_files/S.phase.human.mouse.conversion.rds")
          s.genes_ver2 <- s.genes_ver2$MGI.symbol
          
          g2m.genes_ver2 <- readRDS("cell_cycle_phase_helper_files/G2M.phase.human.mouse.conversion.rds")
          g2m.genes_ver2 <- g2m.genes_ver2$MGI.symbol
          s.genes <- s.genes_ver2
          g2m.genes <- g2m.genes_ver2
        }
        else
        {
          s.genes <- cc.genes$s.genes
          g2m.genes <- cc.genes$g2m.genes
        }
        
        seurat_object <<- CellCycleScoring(seurat_object, s.features = s.genes, g2m.features = g2m.genes, set.ident = F)
        seurat_object@meta.data$CC.Difference <<- seurat_object$S.Score - seurat_object$G2M.Score
        seurat_object@meta.data$Phase <<- factor(seurat_object@meta.data$Phase, levels = c("G1", "S", "G2M"))
        updateMetadata()
        
        output$cellCyclePCA <- renderPlotly(
          {
            if(!is.null(seurat_object) & input$cellCycleReduction != "-")
            {
              meta <- seurat_object@meta.data
              meta$Cell_id <- rownames(meta)
              label_x <- "PC_1"
              label_y <- "PC_2"
              selected_Reduction <- input$cellCycleReduction
              
              if(selected_Reduction == "umap")
              {
                label_x <- "UMAP_1"
                label_y <- "UMAP_2"
              }
              else if(selected_Reduction == "tsne")
              {
                label_x <- "tSNE_1"
                label_y <- "tSNE_2"
              }
              else if(selected_Reduction == "dfm")
              {
                label_x <- "DC_1"
                label_y <- "DC_2"
              }
              else if(selected_Reduction == "phate")
              {
                label_x <- "PHATE_1"
                label_y <- "PHATE_2"
              }
              
              plot1 <- DimPlot(seurat_object, reduction = selected_Reduction)
              plot1_data <- plot1$data
              plot1_data$Cell_id <- rownames(plot1_data)
              plot1_data <- left_join(plot1_data, meta)
              
              p <- ggplot(plot1_data, aes_string(x=label_x, y=label_y, color="Phase"))+
                geom_point() +
                theme_bw() +
                labs(x=label_x, y=label_y)+
                theme(legend.position = "right")
              
              gp <- plotly::ggplotly(p, tooltip = c("x", "y"))
              print(gp)
            }
          }
        )
        
        output$cellCycleBarplot <- renderPlotly(
          {
            #barplot
            fsize <- 24
            
            obj_meta <- seurat_object@meta.data
            obj_meta$cell_id <- row.names(obj_meta)
            obj_meta$Cluster_name <- Idents(seurat_object)
            obj_meta <- obj_meta[, c('cell_id', 'Cluster_name', 'Phase')]
            final_df_seurat <- obj_meta %>% dplyr::group_by(Cluster_name) %>% dplyr::count(Phase)
            colnames(final_df_seurat)[3] <- "Cells"
            
            total_cells_per_cluster <- as.data.frame(table(seurat_object$seurat_clusters)) #as.data.frame(table(Idents(seurat_object)))
            colnames(total_cells_per_cluster)[1] <- "Cluster_name"
            colnames(total_cells_per_cluster)[2] <- "Total"
            
            final_df_seurat <- left_join(final_df_seurat, total_cells_per_cluster)
            final_df_seurat$Percentage <- (final_df_seurat$Cells/final_df_seurat$Total)*100
            
            p <- ggplot(final_df_seurat, aes(x=Phase, y=Percentage, fill=Phase))+
              geom_bar(stat='identity')+
              facet_wrap(~Cluster_name) +
              theme_bw() +
              theme(axis.text.x = element_text(face = "bold", color = "black", size = fsize, angle = 90, vjust = 0.5),
                    axis.text.y = element_text(face = "bold", color = "black", size = 18, angle = 0),
                    axis.title.y = element_text(face = "bold", color = "black", size = fsize),
                    axis.title.x = element_text(face = "bold", color = "black", size = fsize),
                    legend.text = element_text(size = 18, color = "black", face = "bold.italic"),
                    legend.title = element_text(size = 18, color = "black", face = "bold.italic"),
                    strip.text.x = element_text(size = 18, color = "black", face = "bold.italic"),
                    strip.background = element_rect(color="black", fill="#FC4E07", size=1.5, linetype="solid")) +
              labs(y="% of cells", x="Phase", color="Phase")
            p
            plotly::ggplotly(p)
          }
        )
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with cell cycle phase analysis.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Cell_cycle_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "cellCycleRun")
    })
  })
  
  #------------------gProfiler tab-----------------------------------------------
  observeEvent(input$gProfilerConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----gProfiler_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "gProfilerConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@misc$markers, NULL)) session$sendCustomMessage("handler_alert", "Please, execute the gene differential analysis at the MARKERS' IDENTIFICATION tab first.")
      else {
        showModal(modalDialog(div('Functional enrichment analysis in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("gProfilerManhattan_loader")
        
        cluster_temp <- input$gProfilerList
        temp_df <- data.frame()
        # if(cluster_temp == "all_clusters")
        # {
        #   #all clusters used
        #   all_clusters <- as.character(unique(seurat_object@misc$markers$cluster))
        #   
        #   gene_lists <- list()
        #   for(i in 1:length(all_clusters))
        #   {
        #     if(input$gProfilerLFCRadio == "Up")#UP regulated
        #     {
        #       gene_lists[[i]] <- seurat_object@misc$markers[which(seurat_object@misc$markers$cluster == all_clusters[i] & 
        #                                                             seurat_object@misc$markers[, markers_logFCBase] >= as.numeric(input$gProfilerSliderLogFC) &
        #                                                             seurat_object@misc$markers[, input$gprofilerRadio] < as.numeric(input$gProfilerSliderSignificance)), 'gene']
        #     }
        #     else #down
        #     {
        #       gene_lists[[i]] <- seurat_object@misc$markers[which(seurat_object@misc$markers$cluster == all_clusters[i] & 
        #                                                             seurat_object@misc$markers[, markers_logFCBase] <= (as.numeric(input$gProfilerSliderLogFC)*(-1)) &
        #                                                             seurat_object@misc$markers[, input$gprofilerRadio] < as.numeric(input$gProfilerSliderSignificance)), 'gene']
        #     }
        #   }
        #   
        #   names(gene_lists) <- all_clusters
        #   session$sendCustomMessage("handler_startLoader", c("gprof1_loader", 60))
        #   
        #   # gostres <- gost(query = gene_lists, 
        #   #                 organism = "mmusculus", ordered_query = FALSE, 
        #   #                 multi_query = F, significant = TRUE, exclude_iea = T, 
        #   #                 measure_underrepresentation = FALSE, evcodes = TRUE, 
        #   #                 user_threshold = 0.05, correction_method = "g_SCS", 
        #   #                 domain_scope = "annotated", custom_bg = NULL, 
        #   #                 numeric_ns = "", sources = NULL, as_short_link = FALSE)
        #   
        #   gostres <- gost(query = gene_lists, 
        #                   organism = input$gProfilerOrganism, ordered_query = FALSE, 
        #                   multi_query = F, significant = TRUE, exclude_iea = F, 
        #                   measure_underrepresentation = FALSE, evcodes = TRUE, 
        #                   user_threshold = as.numeric(input$gProfilerSliderSignificanceTerms), 
        #                   correction_method = input$gprofilerRadioCorrection, 
        #                   domain_scope = "annotated", custom_bg = NULL, 
        #                   numeric_ns = "", sources = input$gProfilerDatasources, as_short_link = FALSE)
        #   
        #   temp_df <- gostres$result
        #   
        #   # output$gProfilerManhattan <- renderPlotly({ gostplot(gostres, capped = T, interactive = T) %>% 
        #   # layout(width=1500, height=length(all_clusters)*400) })
        #   if (!is.na(all_clusters)){ # TODO check if null or sth else
        #     session$sendCustomMessage("handler_fixHeight", c("gProfilerManhattan", length(all_clusters)))
        #     output$gProfilerManhattan <- renderPlotly({ print(gostplot(gostres, capped = TRUE, interactive = TRUE)) })
        #     session$sendCustomMessage("handler_startLoader", c("gprof2_loader", 80))
        #   } else session$sendCustomMessage("handler_alert", "Please, first calculate clusters") # TODO desicde sentence
        #   #output$gProfilerManhattan <- renderPlotly({ plotly::ggplotly(gostplot(gostres, capped = TRUE, interactive = TRUE)) })
        #   #output$gProfilerManhattan <- renderPlot({gostplot(gostres, capped = TRUE, interactive = FALSE)}, height = 4000)
        # }
        # else
        # {
          gene_lists <- list()
          
          if(input$gProfilerLFCRadio == "Up")#UP regulated
          {
            gene_lists[[1]] <- seurat_object@misc$markers[which(seurat_object@misc$markers$cluster == cluster_temp & 
                                                                  seurat_object@misc$markers[, markers_logFCBase] >= as.numeric(input$gProfilerSliderLogFC) &
                                                                  seurat_object@misc$markers[, input$gprofilerRadio] < as.numeric(input$gProfilerSliderSignificance)), 'gene']
          }
          else #down
          {
            gene_lists[[1]] <- seurat_object@misc$markers[which(seurat_object@misc$markers$cluster == cluster_temp & 
                                                                  seurat_object@misc$markers[, markers_logFCBase] <= (as.numeric(input$gProfilerSliderLogFC)*(-1)) &
                                                                  seurat_object@misc$markers[, input$gprofilerRadio] < as.numeric(input$gProfilerSliderSignificance)), 'gene']
          }
          
          names(gene_lists) <- cluster_temp
          
          gostres <- gost(query = gene_lists, 
                          organism = input$gProfilerOrganism, ordered_query = FALSE, 
                          multi_query = F, significant = TRUE, exclude_iea = F, 
                          measure_underrepresentation = FALSE, evcodes = TRUE, 
                          user_threshold = as.numeric(input$gProfilerSliderSignificanceTerms), 
                          correction_method = input$gprofilerRadioCorrection, 
                          domain_scope = "annotated", custom_bg = NULL, 
                          numeric_ns = "", sources = input$gProfilerDatasources, as_short_link = FALSE)
          
          temp_df <- gostres$result
          output$gProfilerManhattan <- renderPlotly({ print(gostplot(gostres, capped = TRUE, interactive = TRUE)) })
        
        
        output$gProfilerTable <- renderDataTable(temp_df[, c(1, 11, 3:6, 9, 10, 16)], options = list(pageLength = 10), rownames = F)
        export_enrichedTerms_RNA <<- temp_df[, c(1, 11, 3:6, 9, 10, 16)]
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the Enrichment Analysis.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----gProfiler_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "gProfilerConfirm")
    })
  })

observeEvent(input$sendToFlame, {
  tryCatch({
    if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else if (identical(seurat_object@misc$markers, NULL)) session$sendCustomMessage("handler_alert", "Please, execute the gene differential analysis at the MARKERS' IDENTIFICATION tab first.")
    else {
      cluster_temp <- input$gProfilerFlameSelection
      print(cluster_temp)
      class(cluster_temp)
      temp_df <- data.frame()
      gene_lists <- list()
      
      if(input$gProfilerLFCRadio == "Up")#UP regulated
      {
        temp_df <- seurat_object@misc$markers[which(seurat_object@misc$markers$cluster %in% cluster_temp & 
                                                              seurat_object@misc$markers[, markers_logFCBase] >= as.numeric(input$gProfilerSliderLogFC) &
                                                              seurat_object@misc$markers[, input$gprofilerRadio] < as.numeric(input$gProfilerSliderSignificance)), c('gene', 'cluster')]
      }
      else #down
      {
        temp_df <- seurat_object@misc$markers[which(seurat_object@misc$markers$cluster %in% cluster_temp & 
                                                              seurat_object@misc$markers[, markers_logFCBase] <= (as.numeric(input$gProfilerSliderLogFC)*(-1)) &
                                                              seurat_object@misc$markers[, input$gprofilerRadio] < as.numeric(input$gProfilerSliderSignificance)), c('gene', 'cluster')]
      }
      
      gene_lists <- split(temp_df[, 'gene'] , f = temp_df$cluster)
      
      FLAME_API_LINK <- "https://bib.fleming.gr/bib/api/flame"
      jsonBody <- gene_lists
      sendToAPI(FLAME_API_LINK, jsonBody)
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was an error during the process of sending the results to FLAME.")
  })
})

observeEvent(input$findMotifsConfirmATAC, {
  start.time <- Sys.time()
  print(paste0("!!!!-----Motif_start-----", start.time))
  
  session$sendCustomMessage("handler_disableAllButtons", "findMarkersFPConfirmATAC")
  tryCatch({
    if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
    else {
      showModal(modalDialog(div('Motif enrichment analysis in Progress. This operation may take several minutes. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
      shinyjs::show("findMotifsHeatmapATAC_loader")
      shinyjs::show("findMotifsATACTable_loader")
      
      markersPeaks <- getMarkerFeatures(
        ArchRProj = proj_default,
        useMatrix = "PeakMatrix",
        groupBy = input$findMotifsGroupByATAC, #"predictedGroup_Co", #"Clusters",
        bias = c("TSSEnrichment", "log10(nFrags)"),
        testMethod = input$findMarkersPeaksTestATAC,
        logFile =paste0(user_dir, "/Marker_Peaks_file.log")
      )

      proj_default <<- addMotifAnnotations(ArchRProj = proj_default, motifSet = input$findMotifsSetATAC, name = "Motif", force = T,
                                           logFile =paste0(user_dir, "/Motif_Annotation_file.log"))
      print("afterAddMotifAnno")
      enrichMotifs <- peakAnnoEnrichment(
        seMarker = markersPeaks,
        ArchRProj = proj_default,
        peakAnnotation = "Motif",
        cutOff = paste0("FDR <= ",input$findMotifsFDRATAC," & Log2FC >= ",input$findMotifsLogFCATAC), 
        logFile =paste0(user_dir, "/Motif_file.log")
      )
      
      print("afterPeakAno")
      heatmapEM <- plotEnrichHeatmap(enrichMotifs, n = 10, transpose = TRUE, returnMatrix = T,
                                     logFile =paste0(user_dir, "/Motif_Heatmap_file.log"))
      print(head(heatmapEM))
      full_motif_table <- t(plotEnrichHeatmap(enrichMotifs, transpose = TRUE, returnMatrix = T, 
                                              logFile =paste0(user_dir, "/Motif_Heatmap_file.log")))
      print(head(full_motif_table))
      output$findMotifsTableATAC <- renderDataTable(expr = full_motif_table, filter = 'top', options = list(pageLength = 10))
      export_motifs_ATAC <<- full_motif_table
      
      output$findMotifsHeatmapATAC <- renderPlotly(expr = heatmaply(heatmapEM, colors=paletteContinuous(set = "comet", n = 100)))
      session$sendCustomMessage("handler_enableTabs", c("sidebarMenu", " GENE REGULATORY NETWORK\nANALYSIS"))
    }
  }, error = function(e) {
    print(paste("Error :  ", e))
    session$sendCustomMessage("handler_alert", "There was an error with the the detection of eriched motifs.")
  }, finally = {
    end.time <- Sys.time()
    time.taken <- as.numeric (end.time - start.time, units = "mins")
    print(paste0("!!!!-----Motif_end_Mins_taken-----", time.taken))
    
    removeModal()
    Sys.sleep(1)
    session$sendCustomMessage("handler_enableAllButtons", "findMarkersFPConfirmATAC")
  })	
  
})  

output$gProfilerRNAExport <- downloadHandler(
  filename = function() { 
    paste("enrichedTermsTableRNA-", Sys.Date(), ".txt", sep="")
  },
  content = function(file) {
    write.table(export_enrichedTerms_RNA, file, sep = "\t", quote = F, row.names = F)
  })

output$findMotifsATACExport <- downloadHandler(
  filename = function() { 
    paste("motifsTableATAC-", Sys.Date(), ".txt", sep="")
  },
  content = function(file) {
    write.table(export_motifs_ATAC, file, sep = "\t", quote = F, row.names = T)
  })

  #------------------CIPR tab-----------------------------------------------
  observeEvent(input$annotateClustersConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----CIPR_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "annotateClustersConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@misc$markers, NULL)) session$sendCustomMessage("handler_alert", "Please, execute the gene differential analysis at the MARKERS' IDENTIFICATION tab first.")
      else {
        showModal(modalDialog(div('Cluster annotation in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("annotateClustersCIPRDotplot_loader")
        
        marker_genes <- seurat_object@misc$markers
        
        avgexp <- AverageExpression(seurat_object)
        avgexp <- as.data.frame(avgexp$RNA)
        avgexp$gene <- rownames(avgexp)
        input_CIPR <- marker_genes
        
        #Average expresssion methods
        if(input$annotateClustersMethod %in% c("all_genes_spearman", "all_genes_pearson"))
        {
          input_CIPR <- avgexp
        }
        
        CIPR(input_dat = input_CIPR,
             comp_method = input$annotateClustersMethod, 
             reference = input$annotateClustersReference, 
             keep_top_var = as.numeric(input$annotateClustersSlider),
             plot_ind = F,
             plot_top = F)
        
        CIPR_top_results$index <- as.character(CIPR_top_results$index)
        CIPR_top_results$index <- factor(CIPR_top_results$index, levels = as.character(seq(1:length(seurat_object$seurat_clusters))))
        cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(CIPR_top_results$cluster)))
        
        p <- ggplot(CIPR_top_results, aes(x=index, y=identity_score, fill = cluster, size=7)) +
          theme_bw() +
          geom_point(alpha=1, shape=21) +
          scale_size(range = c(7, 7), guide = "none") + 
          scale_x_discrete(labels=CIPR_top_results$long_name) +
          scale_fill_manual(values = cols) +
          labs(x="") + 
          theme(legend.position="top",
                axis.text.x = element_text(vjust=0.5, hjust=1, angle = 90))
        
        output$annotateClustersCIPRTable <- renderDataTable(CIPR_top_results[], options = list(pageLength = 20), rownames = F) #remove Description
        export_annotation_RNA <<- CIPR_top_results[]
        
        output$annotateClustersCIPRDotplot <- renderPlotly({ print(p)})
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with cluster annotation.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----CIPR_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "annotateClustersConfirm")
    })
  })
  
  output$annotationRNAExport <- downloadHandler(
  filename = function() { 
    paste("clusterAnnotationTableRNA-", Sys.Date(), ".txt", sep="")
  },
  content = function(file) {
    write.table(export_annotation_RNA, file, sep = "\t", quote = F, row.names = F)
  })
  
  #scATAC
  observeEvent(input$annotateClustersConfirmATAC, {
    start.time <- Sys.time()
    print(paste0("!!!!-----Integration_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "annotateClustersConfirmATAC")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        showModal(modalDialog(div('scATAC - scRNA-seq integration in progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("annotateClustersUMAP_loader")
        
        #######################################################################################################
        ############################# RNA integration - with no "group" separation ############################
        #######################################################################################################
        shinyjs::show("annotateClustersUMAP_loader")
        seRNA <<- readRDS(input$annotateClustersRDSInput$datapath)
        seRNA$seurat_clusters <<- as.character(seRNA$seurat_clusters)
        table(seRNA$seurat_clusters)
        
        proj_default <<- addGeneIntegrationMatrix(ArchRProj = proj_default, useMatrix = "GeneScoreMatrix", matrixName = "GeneIntegrationMatrix", reducedDims = "IterativeLSI",
                                                  seRNA = seRNA, addToArrow = FALSE, groupRNA = "seurat_clusters", nameCell = "predictedCell_Un", nameGroup = "predictedGroup_Un",
                                                  nameScore = "predictedScore_Un",
                                                  logFile =paste0(user_dir, "/GIM1_file.log")
                                                  )
        
        cM <- as.matrix(confusionMatrix(proj_default$Clusters, proj_default$predictedGroup_Un))
        preClust <- colnames(cM)[apply(cM, 1 , which.max)]
        cbind(preClust, rownames(cM)) #Assignments
        
        unique(unique(proj_default$predictedGroup_Un))
        
        all <- paste0(paste0(unique(unique(proj_default$predictedGroup_Un))), collapse="|")
        all
        clustall <- rownames(cM)[grep(all, preClust)]
        clustall
        
        rnaall <- colnames(seRNA)[grep(all, seRNA$seurat_clusters)]
        groupList <- SimpleList(
          all = SimpleList(
            ATAC = proj_default$cellNames[proj_default$Clusters %in% clustall],
            RNA = rnaall
          )
        )
        
        proj_default <<- addGeneIntegrationMatrix(ArchRProj = proj_default, useMatrix = "GeneScoreMatrix", matrixName = "GeneIntegrationMatrix", reducedDims = "IterativeLSI",
                                                  seRNA = seRNA, addToArrow = FALSE, groupList = groupList, groupRNA = "seurat_clusters", nameCell = "predictedCell_Co",
                                                  nameGroup = "predictedGroup_Co", nameScore = "predictedScore_Co",force = T,
                                                  logFile =paste0(user_dir, "/GIM2_file.log")
                                                  )
        
        pal <- paletteDiscrete(values = as.character(seRNA$seurat_clusters))
        
        proj_default <<- addGeneIntegrationMatrix(ArchRProj = proj_default, useMatrix = "GeneScoreMatrix", matrixName = "GeneIntegrationMatrix", reducedDims = "IterativeLSI",
                                                  seRNA = seRNA, addToArrow = TRUE, groupList = groupList, groupRNA = "seurat_clusters", nameCell = "predictedCell",
                                                  nameGroup = "predictedGroup", nameScore = "predictedScore",force = T,
                                                  logFile =paste0(user_dir, "/GIM3_file.log")
                                                  )
        saveArchRProject(ArchRProj = proj_default)
        
        #p1 <- plotEmbedding(proj_default, colorBy = "cellColData", name = "predictedGroup_Un", pal = pal, embedding = "umap")
        #p2 <- plotEmbedding(proj_default, colorBy = "cellColData", name = "predictedGroup_Co", pal = pal, embedding = "umap")
        
        meta <- as.data.frame(getCellColData(proj_default))
        meta$Cell_id <- rownames(meta)
        reduc_data <- data.frame()
        
        #prepare colors
        cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, "predictedGroup_Co"])))
        
        #for all reductions
        archr_object_reduc <- proj_default@embeddings[["umap"]]$df
        archr_object_reduc <- archr_object_reduc[, c(1, 2)]
        archr_object_reduc$Cell_id <- rownames(archr_object_reduc)
        reduc_data <- left_join(archr_object_reduc, meta)
        print(head(reduc_data))
        
        colnames(reduc_data) <- gsub("IterativeLSI#", "", colnames(reduc_data))
        colnames(reduc_data) <- gsub("Dimension_", "", colnames(reduc_data))
        
        p <- ggplot(data=reduc_data, aes_string(x="UMAP_1", y="UMAP_2", fill="predictedGroup_Co")) +
          geom_point(size= 2, shape=21, alpha= 1, stroke=0)+
          scale_fill_manual(values = cols)+
          scale_size()+
          theme_bw() +
          theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                axis.title.y = element_text(face = "bold", color = "black", size = 25),
                axis.title.x = element_text(face = "bold", color = "black", size = 25),
                legend.text = element_text(face = "bold", color = "black", size = 9),
                legend.title = element_text(face = "bold", color = "black", size = 9),
                legend.position="right",
                title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
          labs(x="UMAP 1", y="UMAP 2", color="Cell type", title = "", fill="Color")
        
        output$annotateClustersUMAPplot <- renderPlotly({plotly::ggplotly(p)})
        
        #output$annotateClustersUMAPplot <- renderPlot({ p2 })
        
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the the annotation of the dataset.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Integration_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "annotateClustersConfirmATAC")
    })
  })

  #--------------------trajectory tab----------------------------------------
  observeEvent(input$trajectoryConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----Trajectory_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "trajectoryConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first and then run UMAP, in the NON LINEAR DIMENSIONALITY REDUCTION tab.")
      else if (identical(seurat_object$umap, NULL)) session$sendCustomMessage("handler_alert", "Please, calculate UMAP first, in the NON LINEAR DIMENSIONALITY REDUCTION tab.")
      else {
        showModal(modalDialog(div('Trajectory analysis in progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("trajectoryPlot_loader")
        shinyjs::show("trajectoryPseudotimePlot_loader")
        
        #delete previous lineages columns
        for_delete <- grep("Lineage", colnames(seurat_object@meta.data))
        if(length(for_delete) != 0)
        {
          seurat_object@meta.data <<- seurat_object@meta.data[, -for_delete]
        }
        
        reduction <- Embeddings(seurat_object, input$trajectoryReduction)
        
        if(ncol(reduction) >= as.numeric(input$trajectorySliderDimensions)) #Out of bounds in the else
        {
          sds <- slingshot(reduction[, 1:as.numeric(input$trajectorySliderDimensions)], clusterLabels = seurat_object$seurat_clusters, 
                           start.clus = input$trajectoryStart, end.clus = input$trajectoryEnd, stretch = 0)
          
          metaD$all_lin <- slingLineages(sds)
          pt <- slingPseudotime(sds)
          
          for(i in 1:ncol(pt))
          {
            temp_lin <- pt[, i]
            seurat_object <<- AddMetaData(
              object = seurat_object,
              metadata = temp_lin,
              col.name = paste0("Lineage", i)
            )
          }
          
          updateInputLineageList(names(metaD$all_lin))
          updateMetadata()
          
          plot_D <- dittoDimPlot(seurat_object, "seurat_clusters",
                                 do.label = TRUE,
                                 labels.repel = TRUE, 
                                 reduction.use = "umap",
                                 add.trajectory.lineages = slingLineages(sds),
                                 trajectory.cluster.meta = "seurat_clusters",
                                 color.panel = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, 'seurat_clusters']))),
                                 data.out = F)
          
          output$trajectoryPlot <- renderPlot({ print(plot_D)})
          
          plot_P <- dittoDimPlot(seurat_object, var = "Lineage1",  
                                 do.label = F,
                                 labels.repel = F, 
                                 reduction.use = "umap",
                                 add.trajectory.lineages = list((metaD$all_lin[["Lineage1"]])),
                                 trajectory.cluster.meta = "seurat_clusters",
                                 size = 2,
                                 data.out = F) + scale_colour_gradientn(colours = plasma(100), na.value = "grey90") +
            labs(colour = "Pseudotime") 
          
          output$trajectoryPseudotimePlot <- renderPlot({ print(plot_P) })
          output$trajectoryText <- renderPrint({print(metaD$all_lin)})
        }
        else
        {
          session$sendCustomMessage("handler_alert", "Out of bounds exception. The number of dimensions selected exceed the dimensions of the chosen reduction.")
        }
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      shinyjs::hide("trajectoryPlot_loader")
      shinyjs::hide("trajectoryPseudotimePlot_loader")
      session$sendCustomMessage("handler_alert", "There was an error with Trajectory Analysis.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Trajectory_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "trajectoryConfirm")
    })
  })
  
  observeEvent(input$trajectoryConfirmLineage, {
    session$sendCustomMessage("handler_disableAllButtons", "trajectoryConfirmLineage")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        #Lineage view
        showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        plot_P <- dittoDimPlot(seurat_object, var = input$trajectoryLineageSelect,
                               do.label = F,
                               labels.repel = F,
                               reduction.use = "umap",
                               add.trajectory.lineages = list(metaD$all_lin[[input$trajectoryLineageSelect]]),
                               trajectory.cluster.meta = "seurat_clusters",
                               size = 2,
                               rename.var.groups = "Pseudotime",
                               data.out = F) + scale_colour_gradientn(colours = plasma(100), na.value = "grey90") +
          labs(colour = "Pseudotime")
        
        output$trajectoryPseudotimePlot <- renderPlot({ print(plot_P) })
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the visualization of the plot.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "trajectoryConfirmLineage")
    })
  })
  
  observeEvent(input$trajectoryConfirmATAC, {
    start.time <- Sys.time()
    print(paste0("!!!!-----Trajectory_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "trajectoryConfirmATAC")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        showModal(modalDialog(div('Trajectory analysis in Progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        cols_to_delete <- grep(pattern = "^Lineage", x = colnames(getCellColData(proj_default)))
        if(length(cols_to_delete) != 0)
        {
          metadata_copy <- getCellColData(proj_default)
          metadata_copy[cols_to_delete] <- NULL
          proj_default@cellColData <<- metadata_copy
        }
        
        #run slingshot
        rD <- getEmbedding(ArchRProj = proj_default, embedding = "umap")
        groups <- getCellColData(ArchRProj = proj_default, select = input$trajectoryGroupByATAC) #"Clusters", #predictedGroup_Co
        
        maxDims <- as.numeric(input$trajectorySliderDimensionsATAC)
        starCl <- input$trajectoryStartATAC
        endCl <- input$trajectoryEndATAC
        
        sds <- slingshot(
          data = rD[, 1:maxDims], 
          clusterLabels = groups[rownames(rD), ], 
          start.clus = starCl, 
          end.clus = endCl
        )
        
        #add lineages to the object
        for(i in 1:length(slingLineages(sds))) #sds@metadata$lineages
        {
          current_lineage <- slingLineages(sds)[[i]] #sds@metadata$lineages[[i]]
          current_name <- paste0("Lineage", i)
          proj_default <<- addTrajectory(proj_default, trajectory = current_lineage, groupBy = input$trajectoryGroupByATAC,  #"Clusters", predictedGroup_Co
                                         name = current_name, force = T, 
                                         logFile =paste0(user_dir, "/Trajectory_file.log"))
        }
        
        #for verbatim text
        output$trajectoryTextATAC <- renderPrint({ print(slingLineages(sds)) }) #sds@metadata$lineages
        
        updateSelectInput(session, "trajectoryLineageSelectATAC", choices = names(slingLineages(sds))) #sds@metadata$lineages
        updateMetadataATAC()
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error during the trajectory analysis.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----Trajectory_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "trajectoryConfirmATAC")
    })
  }) 
  
  observeEvent(input$trajectoryConfirmLineageATAC, {
    session$sendCustomMessage("handler_disableAllButtons", "trajectoryConfirmATAC")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        showModal(modalDialog(div('Rendering in progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("trajectoryPseudotimePlotATAC_loader")
        
        p <- plotTrajectory(proj_default, trajectory = input$trajectoryLineageSelectATAC, colorBy = "cellColData", name = input$trajectoryLineageSelectATAC, embedding = "UMAP", 
                            logFile = paste0(user_dir, "/Plot_trajectory_file.log"))
        output$trajectoryPseudotimePlotATAC <- renderPlot({ plot(p[[1]]) })
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the the trajectory analysis.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "trajectoryConfirmATAC")
    })
  })
  
  observeEvent(input$trajectoryGroupByATAC, 
               if(input$trajectoryGroupByATAC == "Clusters")
               {
                 cluster_names <- unique(proj_default$Clusters)
                 updateSelectInput(session, "trajectoryStartATAC", choices = cluster_names)
                 updateSelectInput(session, "trajectoryEndATAC", choices = cluster_names)
               }
               else if(input$trajectoryGroupByATAC == "predictedGroup_Co")
               {
                 cluster_names <- unique(proj_default$predictedGroup_Co)
                 updateSelectInput(session, "trajectoryStartATAC", choices = cluster_names)
                 updateSelectInput(session, "trajectoryEndATAC", choices = cluster_names)
               }
               )
  
  #--------------------Ligand Receptor tab---------------------------
  observeEvent(input$ligandReceptorConfirm, {
    start.time <- Sys.time()
    print(paste0("!!!!-----LR_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "ligandReceptorConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first.")
      else {
        showModal(modalDialog(div('Ligand - receptor analysis in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("ligandReceptorFullHeatmap_loader")
        shinyjs::show("ligandReceptorCuratedHeatmap_loader")
        #load interactions
        ligand_target_matrix = readRDS("nichenetr_helper_files/ligand_target_matrix.rds")
        lr_network = readRDS("nichenetr_helper_files/lr_network.rds")
        weighted_networks = readRDS("nichenetr_helper_files/weighted_networks.rds")
        
        weighted_networks_lr = weighted_networks$lr_sig %>% inner_join(lr_network %>% distinct(from,to), by = c("from","to"))
        
        if(organism == "mouse")
        {
          lr_network = lr_network %>% mutate(from = convert_human_to_mouse_symbols(from), to = convert_human_to_mouse_symbols(to)) %>% drop_na()
          colnames(ligand_target_matrix) = ligand_target_matrix %>% colnames() %>% convert_human_to_mouse_symbols()
          rownames(ligand_target_matrix) = ligand_target_matrix %>% rownames() %>% convert_human_to_mouse_symbols()
          ligand_target_matrix = ligand_target_matrix %>% .[!is.na(rownames(ligand_target_matrix)), !is.na(colnames(ligand_target_matrix))]
          weighted_networks_lr = weighted_networks_lr %>% mutate(from = convert_human_to_mouse_symbols(from), to = convert_human_to_mouse_symbols(to)) %>% drop_na()
        }
        
        #expressed genes
        ## receiver
        receiver = input$ligandReceptorReciever
        
        expressed_genes_receiver = get_expressed_genes(receiver, seurat_object, pct = 0.10, assay_oi = "RNA")
        ## sender
        sender = input$ligandReceptorSender
        expressed_genes_sender = get_expressed_genes(sender, seurat_object, pct = 0.10, assay_oi = "RNA")
        
        #active ligands-receptors
        ligands = lr_network %>% pull(from) %>% unique()
        receptors = lr_network %>% pull(to) %>% unique()
        
        expressed_ligands = intersect(ligands,expressed_genes_sender)
        expressed_receptors = intersect(receptors,expressed_genes_receiver)
        
        potential_ligands = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) %>% pull(from) %>% unique()
        
        ## L-R analysis
        lr_network_top = lr_network %>% filter(from %in% potential_ligands & to %in% expressed_receptors) %>% distinct(from,to)
        best_upstream_receptors = lr_network_top %>% pull(to) %>% unique()
        
        lr_network_top_df_large = weighted_networks_lr %>% filter(from %in% potential_ligands & to %in% best_upstream_receptors)
        export_ligandReceptor_full_RNA <<- lr_network_top_df_large
        lr_network_top_df = lr_network_top_df_large %>% spread("from","weight",fill = 0)
        
        lr_network_top_matrix = lr_network_top_df %>% dplyr::select(-to) %>% as.matrix() %>% magrittr::set_rownames(lr_network_top_df$to)
        
        dist_receptors = dist(lr_network_top_matrix, method = "binary")
        hclust_receptors = hclust(dist_receptors, method = "ward.D2")
        order_receptors = hclust_receptors$labels[hclust_receptors$order]
        
        dist_ligands = dist(lr_network_top_matrix %>% t(), method = "binary")
        hclust_ligands = hclust(dist_ligands, method = "ward.D2")
        order_ligands_receptor = hclust_ligands$labels[hclust_ligands$order]
        
        order_receptors = order_receptors %>% intersect(rownames(lr_network_top_matrix))
        order_ligands_receptor = order_ligands_receptor %>% intersect(colnames(lr_network_top_matrix))
        
        vis_ligand_receptor_network = lr_network_top_matrix[order_receptors, order_ligands_receptor]
        rownames(vis_ligand_receptor_network) = order_receptors %>% make.names()
        colnames(vis_ligand_receptor_network) = order_ligands_receptor %>% make.names()
        p_ligand_receptor_network = vis_ligand_receptor_network %>% t() %>% make_heatmap_ggplot("Ligands","Receptors", color = "mediumvioletred", x_axis_position = "top",legend_title = "Prior interaction potential")
        #p_ligand_receptor_network
        
        output$ligandReceptorFullHeatmap <- renderPlotly({ plotly::ggplotly(p_ligand_receptor_network) })
        
        #only curated interactions
        lr_network_strict = lr_network %>% filter(database != "ppi_prediction_go" & database != "ppi_prediction")
        ligands_bona_fide = lr_network_strict %>% pull(from) %>% unique()
        receptors_bona_fide = lr_network_strict %>% pull(to) %>% unique()
        
        lr_network_top_df_large_strict = lr_network_top_df_large %>% distinct(from,to) %>% inner_join(lr_network_strict, by = c("from","to")) %>% distinct(from,to)
        lr_network_top_df_large_strict = lr_network_top_df_large_strict %>% inner_join(lr_network_top_df_large, by = c("from","to"))
        export_ligandReceptor_short_RNA <<- lr_network_top_df_large_strict
        lr_network_top_df_strict = lr_network_top_df_large_strict %>% spread("from","weight",fill = 0)
        lr_network_top_matrix_strict = lr_network_top_df_strict %>% dplyr::select(-to) %>% as.matrix() %>% magrittr::set_rownames(lr_network_top_df_strict$to)
        
        dist_receptors = dist(lr_network_top_matrix_strict, method = "binary")
        hclust_receptors = hclust(dist_receptors, method = "ward.D2")
        order_receptors = hclust_receptors$labels[hclust_receptors$order]
        
        dist_ligands = dist(lr_network_top_matrix_strict %>% t(), method = "binary")
        hclust_ligands = hclust(dist_ligands, method = "ward.D2")
        order_ligands_receptor = hclust_ligands$labels[hclust_ligands$order]
        
        order_receptors = order_receptors %>% intersect(rownames(lr_network_top_matrix_strict))
        order_ligands_receptor = order_ligands_receptor %>% intersect(colnames(lr_network_top_matrix_strict))
        
        vis_ligand_receptor_network_strict = lr_network_top_matrix_strict[order_receptors, order_ligands_receptor]
        rownames(vis_ligand_receptor_network_strict) = order_receptors %>% make.names()
        colnames(vis_ligand_receptor_network_strict) = order_ligands_receptor %>% make.names()
        p_ligand_receptor_network_strict = vis_ligand_receptor_network_strict %>% t() %>% make_heatmap_ggplot("Ligands","Receptors", color = "brown", x_axis_position = "top",legend_title = "Prior interaction potential\n(bona fide)")
        # p_ligand_receptor_network_strict
        
        output$ligandReceptorCuratedHeatmap <- renderPlotly({ plotly::ggplotly(p_ligand_receptor_network_strict) })
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error during ligand - receptor analysis.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----LR_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "ligandReceptorConfirm")
    })
  })
  
  output$ligandReceptorFullExport <- downloadHandler(
    filename = function() { 
      paste("ligandReceptorFullTableRNA-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_ligandReceptor_full_RNA, file, sep = "\t", quote = F, row.names = F)
    })
  
  output$ligandReceptorShortExport <- downloadHandler(
    filename = function() { 
      paste("ligandReceptorStrictTableRNA-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_ligandReceptor_short_RNA, file, sep = "\t", quote = F, row.names = F)
    })
  
  #---------------------------GRN tab-------------------------------------------
  observeEvent(input$grnProduceLoom, {
    session$sendCustomMessage("handler_disableAllButtons", "grnProduceLoom")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first.")
      else {
        showModal(modalDialog(div('Analysis in Progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        genome_info <- input$grnGenomeBuild
        print(genome_info)
        seurat_object_matrix<-seurat_object@assays$RNA@counts
        
        # Filter 1: this filter counts how many UMIs per gene we have, and keeps genes with over ncol(exprMat)*.005*1 UMIs ( 30.015 in this case )
        
        nCountsPerGene.seurat_object<- rowSums(seurat_object_matrix, na.rm = T)
        genesLeft_minReads.seurat_object<- names(nCountsPerGene.seurat_object[which(nCountsPerGene.seurat_object>(ncol(seurat_object_matrix)*.01*3))])  
        length(genesLeft_minReads.seurat_object)
        
        # Filter 2: this filter counts for each gene, how many "cells" with at least one UMI  we have, and keeps genes with over ncol(exprMat)*.01 cells ( 60.03 in this case ). This filter should be very stringent in the particular project, because we will lose a very important gene - "Epcam" (with 7 cells) - that is very important in the cluster of hepatocytes (a cluster of 4 cells in 10x? ask Irina). So we can change the filter to ncol(exprMat)*.00084
        
        nCellsPerGene.seurat_object<- rowSums(seurat_object_matrix>0, na.rm = T)
        tmp <- nCellsPerGene.seurat_object[genesLeft_minReads.seurat_object]
        genesLeft_minCells.seurat_object<- names(tmp[which(tmp>(ncol(seurat_object_matrix)*.01))])
        length(genesLeft_minCells.seurat_object)
        
        # Filter 3 - Exclude genes missing from database:
        print("import rankings")
        print(user_dir)
        
        motifRankings1_mouse_human <- ""
        motifRankings2_mouse_human <- ""
        
        if(genome_info == "mm10")
        {
          print("mm10 if")
          motifRankings1_mouse_human <- importRankings("scenic_helper_files/mm10__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather") # either one, they should have the same genes
          print("finished 1st importing")
          motifRankings2_mouse_human <- importRankings("scenic_helper_files/mm10__refseq-r80__10kb_up_and_down_tss.mc9nr.feather") # either one, they should have the same genes
          print("finished importing")
        }
        else if(genome_info == "hg19")
        {
          print("hg19 if")
          motifRankings1_mouse_human <- importRankings("scenic_helper_files/hg19-500bp-upstream-10species.mc9nr.feather") # either one, they should have the same genes
          print("finished 1st importing")
          motifRankings2_mouse_human <- importRankings("scenic_helper_files/hg19-tss-centered-10kb-10species.mc9nr.feather") # either one, they should have the same genes
          print("finished importing")
        }
        else
        {
          motifRankings1_mouse_human <- importRankings("scenic_helper_files/hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather") # either one, they should have the same genes
          motifRankings2_mouse_human <- importRankings("scenic_helper_files/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather") # either one, they should have the same genes
        }
        print("get rankings")
        genesInDatabase1 <- colnames(getRanking(motifRankings1_mouse_human))
        genesInDatabase2 <- colnames(getRanking(motifRankings2_mouse_human))
        
        
        genesLeft_minCells_inDatabases1_seurat_object <- genesLeft_minCells.seurat_object[which(genesLeft_minCells.seurat_object %in% colnames(motifRankings1_mouse_human@rankings))]
        genesLeft_minCells_inDatabases2_seurat_object <- genesLeft_minCells.seurat_object[which(genesLeft_minCells.seurat_object %in% colnames(motifRankings2_mouse_human@rankings))]
        
        genesKept1_seurat_object <- genesLeft_minCells_inDatabases1_seurat_object
        genesKept2_seurat_object <- genesLeft_minCells_inDatabases2_seurat_object
        
        exprMat_filtered2_seurat_object<-seurat_object_matrix[which(rownames(seurat_object_matrix) %in% genesKept2_seurat_object),]
        
        seurat.umap <- Embeddings(seurat_object[["umap"]])[, c(1,2)]
        print(head(seurat.umap))
        print(user_dir)
        file_name <- paste0(user_dir,"/", genome_info, "_seurat.loom")
        
        loom <- build_loom(file.name= file_name,
                           dgem=exprMat_filtered2_seurat_object,
                           title="seurat_object",
                           genome= genome_info, # Just for user information, not used internally
                           default.embedding=seurat.umap,
                           default.embedding.name="UMAP")
        finalize(loom=loom)
        print("Finalized loom")
        removeModal()
        Sys.sleep(1)
        shinyalert(title = "Files for pyscenic run are ready for download!", type = "success")
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      removeModal()
      Sys.sleep(1)
      shinyalert(title = "An error has occured during the preparation of files.", type = "error")
      #session$sendCustomMessage("handler_alert", "There was an error with generating pyscenic input files.")
    }, finally = {
      session$sendCustomMessage("handler_enableAllButtons", "grnProduceLoom")
    })
  })
  
  output$grnDownloadLoom <- downloadHandler(
    filename = function() {
      paste0(user_dir,"/", input$grnGenomeBuild, "_seurat.loom")
    },
    content <- function(file) {
      file.copy(paste0(user_dir,"/", input$grnGenomeBuild, "_seurat.loom"), file)
    },
    contentType = "application/loom"
    )
  
  output$grnDownloadRDS <- downloadHandler(
    filename = function() { 
      paste("processed_seurat_object-", Sys.Date(), ".RDS", sep="")
    },
    content = function(file) {
      saveRDS(seurat_object, file)
    })
  
  output$grnMatrixScenicRNAExport <- downloadHandler(
    filename = function() { 
      paste("ScenicAUCvaluesRNA-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_scenicAUC_full_RNA, file, sep = "\t", quote = F, row.names = T)
    })
  
  #-----analysis
  observeEvent(input$grnLoomAnalysis, {
    session$sendCustomMessage("handler_disableAllButtons", "grnLoomAnalysis")
    tryCatch({
      showModal(modalDialog(div('Analysis in Progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        source('aux_rss.R')
        userId <- session$token
        user_dir_pyscenic <<- paste0("./usr_temp/pysc_", userId, gsub(pattern = "[ ]|[:]", replacement = "_", x = paste0("_", Sys.time()))) 
        dir.create(user_dir_pyscenic)
        
        numOfCores <- 1
        file.copy(from = input$grnLoomInput$datapath, to = paste0(user_dir_pyscenic, "/AUC.loom"), overwrite = TRUE)
        file.copy(from = input$grnRDSInput$datapath, to = paste0(user_dir_pyscenic, "/seurat_object.RDS"), overwrite = TRUE)
        
        seurat_object <<- readRDS(paste0(user_dir_pyscenic, "/seurat_object.RDS"))
        pyScenicLoomFile.seurat_object <- file.path(paste0(user_dir_pyscenic, "/AUC.loom"))
        loom.seurat_object <- open_loom(pyScenicLoomFile.seurat_object, mode="r+")
        
        # Read information from loom file:
        regulons_incidMat.seurat_object <-get_regulons(loom.seurat_object, column.attr.name = "Regulons", tf.as.name = TRUE, tf.sep = "_")
        regulons.seurat_object <- regulonsToGeneLists(regulons_incidMat.seurat_object)
        regulonsAUC.seurat_object <- get_regulons_AUC(loom.seurat_object, column.attr.name='RegulonsAUC')
        regulonsAucThresholds.seurat_object <- get_regulon_thresholds(loom.seurat_object)
        embeddings.seurat_object <- get_embeddings(loom.seurat_object)
        exprMat.seurat_object <- get_dgem(loom.seurat_object)
        add_col_attr(loom=loom.seurat_object, key = "seurat_clusters", value=as.character(seurat_object$seurat_clusters), as.annotation=T)
        cellInfo.seurat_object <- get_cell_annotation(loom.seurat_object)
        
        close_loom(loom.seurat_object)
        cells_rankings.seurat_object <- AUCell_buildRankings(as.matrix(exprMat.seurat_object), nCores=numOfCores)
        
        geneSets.seurat_object <- regulons.seurat_object
        names(geneSets.seurat_object)<-paste0(names(geneSets.seurat_object)," ",lengths(geneSets.seurat_object),sep="")
        
        geneSets.seurat_object<-geneSets.seurat_object[unlist(lapply(geneSets.seurat_object, function(x) length(x) > 10))]
        
        title<-c("TF","gene_set")
        write.table(t(title),paste0(user_dir_pyscenic, "/regulons.seurat_object.txt"),sep="\t",quote=FALSE,col.names=FALSE,row.names=FALSE)
        temp<-as.data.frame(geneSets.seurat_object[1])
        colnames(temp)<-names(geneSets.seurat_object[1])
        write.table(t(temp),paste0(user_dir_pyscenic, "/regulons.seurat_object.txt"),sep=",",quote=FALSE,col.names=FALSE,row.names=TRUE,append=T)
        for(i in 2:length(geneSets.seurat_object)){
          temp<-as.data.frame(geneSets.seurat_object[i])
          colnames(temp)<-names(geneSets.seurat_object[i])
          write.table(t(temp),paste0(user_dir_pyscenic, "/regulons.seurat_object.txt"),sep=",",quote=FALSE,col.names=FALSE,row.names=TRUE,append=T)
        }
        
        cells_AUC.seurat_object <- AUCell_calcAUC(geneSets.seurat_object, cells_rankings.seurat_object, aucMaxRank=nrow(cells_rankings.seurat_object)*0.1,nCores=numOfCores)
        cells_AUC.seurat_object.ordered <- orderAUC(cells_AUC.seurat_object) # added to AUCell 1.5.1
        cells_AUC.seurat_object <- cells_AUC.seurat_object[cells_AUC.seurat_object.ordered,]
        getAUC(cells_AUC.seurat_object)
        write.table(getAUC(cells_AUC.seurat_object),paste0(user_dir_pyscenic, "/AUC_per_cell.txt"),quote=F,sep="\t",row.names=T,col.names=T)
        cells_assignment.seurat_object <- AUCell_exploreThresholds(cells_AUC.seurat_object, plotHist=F, nCores=1, assign=TRUE)
        
        warningMsg.seurat_object <- sapply(cells_assignment.seurat_object, function(x) x$aucThr$comment)
        warningMsg.seurat_object[which(warningMsg.seurat_object!="")]
        
        regulonsCells.seurat_object <- getAssignments(cells_assignment.seurat_object)
        regulonsCells.seurat_object.filtered<-regulonsCells.seurat_object[unlist(lapply(regulonsCells.seurat_object, function(x) length(x) > 10))]
        cells_assignment.seurat_object.filtered<-cells_assignment.seurat_object[names(regulonsCells.seurat_object.filtered)]
        
        cellsAssigned.seurat_object <- lapply(cells_assignment.seurat_object, function(x) x$assignment)
        assignmentTable.seurat_object <- reshape2::melt(cellsAssigned.seurat_object, value.name="cell")
        colnames(assignmentTable.seurat_object)[2] <- "geneSet"
        
        assignmentMat.seurat_object <- table(assignmentTable.seurat_object[,"geneSet"], assignmentTable.seurat_object[,"cell"])
        
        library(stringr)
        maxRows <- 20
        
        aucMatrix.seurat_object <- t(assignmentMat.seurat_object)
        aucMatrix.seurat_object <- apply(aucMatrix.seurat_object, 2, as.numeric)
        system.time(rss.seurat_object <- calcRSS(t(aucMatrix.seurat_object), seurat_object$seurat_clusters))
        rss.seurat_object <- rss.seurat_object[onlyNonDuplicatedExtended(rownames(rss.seurat_object)),]
        
        
        cells_AUC.seurat_object.ordered.filtered<-cells_AUC.seurat_object
        regulonActivity_byCellType.seurat_object <- sapply(split(rownames(cellInfo.seurat_object), cellInfo.seurat_object$seurat_clusters),
                                                           function(cells) rowMeans(getAUC(cells_AUC.seurat_object.ordered.filtered)[,cells]))
        regulonActivity_byCellType.seurat_object.non.zero<-regulonActivity_byCellType.seurat_object[which(rownames(regulonActivity_byCellType.seurat_object) %in% rownames(rss.seurat_object)),]
        regulonActivity_byCellType.seurat_object.non.zero<-regulonActivity_byCellType.seurat_object.non.zero[order(rownames(regulonActivity_byCellType.seurat_object.non.zero)),]
        colnames(regulonActivity_byCellType.seurat_object.non.zero)<-paste0("AUC_",colnames(regulonActivity_byCellType.seurat_object.non.zero),sep="")
        
        geneSets.seurat_object.new<-geneSets.seurat_object[which(names(geneSets.seurat_object) %in% rownames(rss.seurat_object))]
        geneSets.seurat_object.new<-geneSets.seurat_object.new[order(names(geneSets.seurat_object.new))]
        rss.seurat_object.new<-as.data.frame(rss.seurat_object[order(rownames(rss.seurat_object)),])
        colnames(rss.seurat_object.new)<-paste0("RSS_",colnames(rss.seurat_object.new),sep="")
        rss.seurat_object.new$TF<-rownames(rss.seurat_object.new)
        
        geneSets.seurat_object.new.xls<-as.data.frame(t(c("TF","gene_set")))
        colnames(geneSets.seurat_object.new.xls)<-c("TF","gene_set")
        for(i in 1:length(geneSets.seurat_object.new)){
          temp<-as.data.frame(geneSets.seurat_object.new[i])
          colnames(temp)<-names(geneSets.seurat_object.new[i])
          temp2<-as.data.frame(t(c(colnames(temp),str_c(temp[,1],collapse=","))))
          colnames(temp2)<-c("TF","gene_set")
          geneSets.seurat_object.new.xls<-rbind(geneSets.seurat_object.new.xls,temp2)
        }
        write.table(geneSets.seurat_object.new.xls,paste0(user_dir_pyscenic, "/regulons.seurat_object.txt"),col.names=T,row.names=F,sep="\t",quote=F)
        
        cellsAssigned.seurat_object <- lapply(cells_assignment.seurat_object.filtered, function(x) x$assignment)
        assignmentTable.seurat_object <- reshape2::melt(cellsAssigned.seurat_object, value.name="cell")
        colnames(assignmentTable.seurat_object)[2] <- "geneSet"
        
        assignmentMat.seurat_object <- table(assignmentTable.seurat_object[,"geneSet"], assignmentTable.seurat_object[,"cell"])
        
        aucMatrix.seurat_object <- t(assignmentMat.seurat_object)
        aucMatrix.seurat_object <- apply(aucMatrix.seurat_object, 2, as.numeric)
        system.time(rss.seurat_object <- calcRSS(t(aucMatrix.seurat_object), seurat_object$seurat_clusters))
        rss.seurat_object <- rss.seurat_object[onlyNonDuplicatedExtended(rownames(rss.seurat_object)),]
        
        top_regulons<-unique(as.vector(apply(rss.seurat_object,2,function(x) names(sort(x,decreasing=TRUE)[1:nrow(rss.seurat_object)]))))
        cells_AUC.seurat_object.ordered.filtered<-cells_AUC.seurat_object[which(rownames(cells_AUC.seurat_object) %in% top_regulons),]
        regulonActivity_byCellType.seurat_object <- sapply(split(rownames(cellInfo.seurat_object), cellInfo.seurat_object$seurat_clusters),
                                                           function(cells) rowMeans(getAUC(cells_AUC.seurat_object.ordered.filtered)[,cells]))
        regulonActivity_byCellType_Scaled.seurat_object <- t(scale(t(regulonActivity_byCellType.seurat_object), center = T, scale=T))
        
        write.table(regulonActivity_byCellType_Scaled.seurat_object, paste0(user_dir_pyscenic, "/scaled_regulon_activity_by_cell_type_FULL_TABLE.txt"), sep = "\t", quote = F, col.names = T)
        write.table(rss.seurat_object, paste0(user_dir_pyscenic, "/rss_regulon_by_cell_type_FULL_TABLE.txt"), sep = "\t", quote = F, col.names = T)
        updateSliderInput(session, "grnTopRegulonsRNA", max = nrow(rss.seurat_object))
        export_scenicAUC_full_RNA <<- regulonActivity_byCellType_Scaled.seurat_object
        print("All done!")
        removeModal()
        Sys.sleep(1)
        shinyalert(title = "Pyscenic files were analyzed successfully! You can now visualize the results!", type = "success")
    }, error = function(e) {
      print(paste("Error :  ", e))
      removeModal()
      Sys.sleep(1)
      shinyalert(title = "An error has occured during the preparation of files.", type = "error")
      #session$sendCustomMessage("handler_alert", "There was an error with analyzing pyscenic results.")
    }, finally = {
      session$sendCustomMessage("handler_enableAllButtons", "grnLoomAnalysis")
    })
  })
    
    observeEvent(input$grnConfirmVisualizationRNA, {
      session$sendCustomMessage("handler_disableAllButtons", "grnConfirmVisualizationRNA")
      tryCatch({
        if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
        else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first.")
        else {
          showModal(modalDialog(div('Analysis in Progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
          shinyjs::show("grnHeatmapRNA_loader")
          
          rss.seurat_object <- read.delim(paste0(user_dir_pyscenic,"/rss_regulon_by_cell_type_FULL_TABLE.txt"), check.names = F)
          regulonActivity_byCellType_Scaled.seurat_object <- read.delim(paste0(user_dir_pyscenic,"/scaled_regulon_activity_by_cell_type_FULL_TABLE.txt"), check.names = F)
          
          #heatmap input
          top_regulons<-unique(as.vector(apply(rss.seurat_object,2,function(x) names(sort(x,decreasing=TRUE)[1:as.numeric(input$grnTopRegulonsRNA)]))))
          regulonActivity_byCellType_Scaled.seurat_object_top <- regulonActivity_byCellType_Scaled.seurat_object[top_regulons, ]
          #matrix input
          mat_plotted <- as.data.frame(regulonActivity_byCellType_Scaled.seurat_object)
          if(input$grnMatrixSelectionRNA == "rss")
          {
            mat_plotted <- as.data.frame(rss.seurat_object)
          }
          
          output$grnMatrixRNA <- renderDataTable(mat_plotted, options = list(pageLength = 10), rownames = T)
          
          output$grnHeatmapRNA <- renderPlotly(heatmaply(t(regulonActivity_byCellType_Scaled.seurat_object_top), colors=viridis(n = 256, option = "plasma")))
        }
      }, error = function(e) {
        print(paste("Error :  ", e))
        Sys.sleep(1)
        shinyalert(title = "An error has occured during the visualization of results.", type = "error")
        #session$sendCustomMessage("handler_alert", "There was an error with visualizing pyscenic results.")
      }, finally = {
        removeModal()
        Sys.sleep(1)
        session$sendCustomMessage("handler_enableAllButtons", "grnConfirmVisualizationRNA")
      })
    })
  
  observeEvent(input$grnRunDecoupler, {
    start.time <- Sys.time()
    print(paste0("!!!!-----GRN_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "grnRunDecoupler")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first.")
      else {
        showModal(modalDialog(div('Analysis in Progress. This operation may take several minutes, please wait...', 
                                  style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        complexes <- input$grnComplexes
        if(complexes == "keep") # split complexes = false
        {
          if (organism == "human")
            net <- read.delim("decoupleR_helper_files/collecTRI_human_complexes_false.txt")
          else
            net <- read.delim("decoupleR_helper_files/collecTRI_mouse_complexes_false.txt")
        } else
        {
          if(organism == "human")
            net <- read.delim("decoupleR_helper_files/collecTRI_human_complexes_true.txt")
          else
            net <- read.delim("decoupleR_helper_files/collecTRI_mouse_complexes_true.txt")
        }
        
        # Extract the normalized log-transformed counts
        mat <- as.matrix(seurat_object@assays$RNA@data)
        
        analysisMethod <- input$grnAnalysisMethod
        
        if(analysisMethod == "mlm")
        {
          # Run mlm
          acts <- run_mlm(mat=mat, net=net, .source='source', .target='target',
                          .mor='mor', minsize = 5)
        } else if (analysisMethod == "ulm")
        {
          # Run ulm
          acts <- run_ulm(mat=mat, net=net, .source='source', .target='target',
                          .mor='mor', minsize = 5)
        } else
        {
          # Run wsum
          acts <- run_wsum(mat=mat, net=net, .source='source', .target='target',
                           .mor='mor', minsize = 5)
          acts <- acts %>% filter(statistic == "norm_wsum")
        }
        
        # Extract ulm and store it in tfsulm in seurat object
        seurat_object[['tfactivity']] <<- acts %>%
          pivot_wider(id_cols = 'source', names_from = 'condition',
                      values_from = 'score') %>%
          column_to_rownames('source') %>%
          Seurat::CreateAssayObject(.)
        
        # Scale the data
        seurat_object <<- ScaleData(seurat_object, assay = "tfactivity")
        seurat_object@assays$tfactivity@data <<- seurat_object@assays$tfactivity@scale.data
        
        n_tfs <- 50
        # Extract activities from object as a long dataframe
        df <- t(as.matrix(seurat_object@assays$tfactivity@data)) %>%
          as.data.frame() %>%
          mutate(cluster = seurat_object$seurat_clusters) %>%
          pivot_longer(cols = -cluster, names_to = "source", values_to = "score") %>%
          group_by(cluster, source) %>%
          dplyr::summarise(mean=mean(score))
        
        # Get top tfs with more variable means across clusters
        tfs <- df %>%
          group_by(source) %>%
          dplyr::summarise(std = sd(mean)) %>%
          arrange(-abs(std)) %>%
          head(n_tfs) %>%
          pull(source)
        
        # Subset long data frame to top tfs and transform to wide matrix
        top_acts_mat <- df %>%
          filter(source %in% tfs) %>%
          pivot_wider(id_cols = 'cluster', names_from = 'source',
                      values_from = 'mean') %>%
          column_to_rownames('cluster') %>%
          as.matrix()
        
        #plot the full table
        df_toPlot <- df
        colnames(df_toPlot) <- c("cluster", "TF", "Mean activity(z-score)")
        output$grnMatrixRNA_DecoupleR <- renderDataTable(df_toPlot, options = list(pageLength = 10), rownames = T)
        export_decoupleRZscores_full_RNA <<- df_toPlot
        
        # Choose color palette
        palette_length = 100
        my_color = colorRampPalette(c("Darkblue", "white","red"))(palette_length)
        
        output$grnHeatmapRNA_DecoupleR <- renderPlotly(heatmaply(x = top_acts_mat, colors = my_color, limits=c(-3, 3)))
        updateUtilitiesAssays()
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with TF activity inference.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----GRN_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "grnRunDecoupler")
    })
  })  
  
  output$grnMatrixDecoupleRRNAExport <- downloadHandler(
    filename = function() { 
      paste("DecoupleRActivityZscoreRNA-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_decoupleRZscores_full_RNA, file, sep = "\t", quote = F, row.names = F)
    })
    
  observeEvent(input$grnConfirmATAC, {
    start.time <- Sys.time()
    print(paste0("!!!!-----GRN_start-----", start.time))
    
    session$sendCustomMessage("handler_disableAllButtons", "grnConfirmATAC")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        showModal(modalDialog(div('GRN analysis in progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("grnHeatmapATAC_loader")
        shinyjs::show("grnATACTable_loader")
        shinyjs::show("grnATACTable2_loader")
        shinyjs::show("grnATACTable3_loader")
        
        fdr_lim <- input$grnFdrATAC
        corr_lim <- input$grnCorrlationATTAC
        
        ##########################################################################
        ######################### Positive regulators ############################
        ##########################################################################
        proj_default <<- addBgdPeaks(proj_default, method = "ArchR")#---local version
        #proj_default <<- addMotifAnnotations(ArchRProj = proj_default, motifSet = input$findMotifsSetATAC, name = "Motif", force = TRUE, 
        #                                     logFile =paste0(user_dir, "/Add_motifs_file.log"))
        proj_default <<- addDeviationsMatrix(ArchRProj = proj_default, peakAnnotation = "Motif", force = TRUE, 
                                             logFile =paste0(user_dir, "/Deveations_file.log"))
        ## Deviant Motifs ##
        seGroupMotif_proj_default_condition <- getGroupSE(ArchRProj = proj_default, useMatrix = "MotifMatrix", groupBy = input$grnGroupByATAC,  #"Clusters", "predictedGroup_Co"
                                                          logFile =paste0(user_dir, "/GroupSE_file.log"))
        ## Keep z-scores or deviation scores ##
        seZ_proj_default_condition <- seGroupMotif_proj_default_condition[rowData(seGroupMotif_proj_default_condition)$seqnames=="z",]
        ## Maximum delta in z-score between all clusters ##
        rowData(seZ_proj_default_condition)$maxDelta <- lapply(seq_len(ncol(seZ_proj_default_condition)), function(x){
          rowMaxs(assay(seZ_proj_default_condition) - assay(seZ_proj_default_condition)[,x],na.rm=TRUE)
        }) %>% Reduce("cbind", .) %>% rowMaxs
        ## Correlate TF Accessibility with Genescore ##
        corGSM_MM_proj_default <<- correlateMatrices(
          ArchRProj = proj_default,
          useMatrix1 = input$grnMatrixATAC, #"GeneScoreMatrix", "GeneIntegrationMatrix"
          useMatrix2 = "MotifMatrix",
          reducedDims = "IterativeLSI",
          logFile =paste0(user_dir, "/Correlation_file.log")
        )
        ## Add Maximum Delta Deviation to the Correlation Data Frame ##
        corGSM_MM_proj_default_condition<-corGSM_MM_proj_default
        corGSM_MM_proj_default_condition$maxDelta <- rowData(seZ_proj_default_condition)[match(corGSM_MM_proj_default_condition$MotifMatrix_name, rowData(seZ_proj_default_condition)$name), "maxDelta"]
        ## Identify Positive TF Regulators ##
        corGSM_MM_proj_default_condition<-corGSM_MM_proj_default_condition[!is.na(corGSM_MM_proj_default_condition$maxDelta),]  
        corGSM_MM_proj_default_condition <- corGSM_MM_proj_default_condition[order(abs(corGSM_MM_proj_default_condition$cor), decreasing = TRUE), ]
        corGSM_MM_proj_default_condition <- corGSM_MM_proj_default_condition[which(!duplicated(gsub("\\-.*","",corGSM_MM_proj_default_condition[,"MotifMatrix_name"]))), ]
        corGSM_MM_proj_default_condition$TFRegulator <- "NO"
        #less strict---
        corGSM_MM_proj_default_condition$TFRegulator[which(corGSM_MM_proj_default_condition$cor > corr_lim & corGSM_MM_proj_default_condition$padj < fdr_lim )] <- "YES"
        #--
        #else--
        #corGSM_MM_proj_default_condition$TFRegulator[which(corGSM_MM_proj_default_condition$cor > corr_lim & corGSM_MM_proj_default_condition$padj < fdr_lim & corGSM_MM_proj_default_condition$maxDelta)] <- "YES"
        #--
        positive_regulators_GSM_MM_proj_default_condition<-sort(corGSM_MM_proj_default_condition[corGSM_MM_proj_default_condition$TFRegulator=="YES",1])
        ## Create the zscore data frames ##
        seZ_proj_default_condition_df<-assays(seZ_proj_default_condition)$MotifMatrix
        rownames(seZ_proj_default_condition_df)<-rowData(seZ_proj_default_condition)$name
        ## Keep positive regulators ##
        seZ_proj_default_condition_df_positive_regulators_GSM_MM_condition<-seZ_proj_default_condition_df[which(sub("_.*", "", rownames(seZ_proj_default_condition_df)) %in% positive_regulators_GSM_MM_proj_default_condition),]
        
        ### P2G links ###
        proj_default <<- addCoAccessibility(ArchRProj = proj_default, reducedDims = "IterativeLSI", 
                                           logFile =paste0(user_dir, "/Coaccessibility_file.log"))
        proj_default <<- addPeak2GeneLinks(ArchRProj = proj_default,reducedDims = "IterativeLSI", useMatrix = "GeneScoreMatrix", 
                                           logFile =paste0(user_dir, "/P2G_file.log")) #addPeak2GeneLinks_shiny
        p <- plotPeak2GeneHeatmap(ArchRProj = proj_default, nPlot = 5000, groupBy = input$grnGroupByATAC, #"Clusters", predictedGroup_Co
                                  returnMatrices = T, corCutOff = corr_lim, FDRCutOff = fdr_lim, 
                                  logFile =paste0(user_dir, "/P2Gheatmap_file.log"))
        p2g_table <- as.data.frame(p$Peak2GeneLinks)
        
        output$grnMatrixATAC <- renderDataTable(seZ_proj_default_condition_df_positive_regulators_GSM_MM_condition, options = list(pageLength = 10), rownames = T)
        export_positiveRegulators_ATAC <<- seZ_proj_default_condition_df_positive_regulators_GSM_MM_condition
        
        ### Peak x Motif matrix ###
        motifMatrix_proj_default<-SummarizedExperiment::assay(getMatches(proj_default))
        rownames(motifMatrix_proj_default) <- paste0(seqnames(getMatches(proj_default)), ":", start(getMatches(proj_default)), "-", end(getMatches(proj_default)))
        motifMatrix_proj_default <- as.matrix(motifMatrix_proj_default) 
        export_PeakMotifTable_ATAC <<- motifMatrix_proj_default
        print(head(motifMatrix_proj_default))
        #-----------------------------------------------------------------------
        
        pos_reg_htable <- t(seZ_proj_default_condition_df_positive_regulators_GSM_MM_condition)
        pos_reg_htable[pos_reg_htable>2] <- 2
        pos_reg_htable[pos_reg_htable< -2] <- -2
        
        output$grnHeatmapATAC <- renderPlotly(heatmaply(pos_reg_htable, colors=viridis(n = 100, option = "plasma")))
        
        p2g_table <- p2g_table[, c(1, 2, 7, 8, 3:6)]
        output$grnP2GlinksTable <- renderDataTable(p2g_table, options = list(pageLength = 10), rownames = T)
        export_peakToGenelinks_ATAC <<- p2g_table
        
        output$grnMotifTable <- renderDataTable(motifMatrix_proj_default[1:10, 1:10], options = list(pageLength = 10), rownames = T)
        
        ##########################################################################
        ##########################################################################
        ##########################################################################
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the the GRN analysis.")
    }, finally = {
      end.time <- Sys.time()
      time.taken <- as.numeric (end.time - start.time, units = "mins")
      print(paste0("!!!!-----GRN_end_Mins_taken-----", time.taken))
      
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "grnConfirmATAC")
    })
  })
  
  output$grnPositiveRegulatorsATACExport <- downloadHandler(
    filename = function() { 
      paste("positiveRegulatorsTableATAC-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_positiveRegulators_ATAC, file, sep = "\t", quote = F, row.names = T)
    })
  
  output$grnPeakToGeneLinksATACExport <- downloadHandler(
    filename = function() { 
      paste("peakToGeneLinksTableATAC-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_peakToGenelinks_ATAC, file, sep = "\t", quote = F, row.names = T)
    })
  
  output$grnPeakMotifTableATACExport <- downloadHandler(
    filename = function() { 
      paste("PeakMotifTableATAC-", Sys.Date(), ".txt", sep="")
    },
    content = function(file) {
      write.table(export_PeakMotifTable_ATAC, file, sep = "\t", quote = F, row.names = T)
    })
  
  #---------------------------Tracks tab----------------------------------------
  observeEvent(input$visualizeTracksConfirm, {
    session$sendCustomMessage("handler_disableAllButtons", "visualizeTracksConfirm")
    tryCatch({
      if (identical(proj_default, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else {
        showModal(modalDialog(div('Rendering in progress. Please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("visualizeTracksOutput_loader")
        
        p <- plotBrowserTrack(
          ArchRProj = proj_default,
          groupBy = input$visualizeTracksGroupByATAC, #"Clusters", predictedGroup_Co
          geneSymbol = input$visualizeTracksGene,
          upstream = as.numeric(input$visualizeTracksBPupstream),
          downstream = as.numeric(input$visualizeTracksBPdownstream),
          baseSize = 15, 
          facetbaseSize = 10, 
          sizes = c(10, 4, 3, 4), 
          logFile =paste0(user_dir, "/Tracks_file.log")
        )
        output$visualizeTracksOutput <- renderPlot({ plot(p[[input$visualizeTracksGene]]) })
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with the generation of tracks.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "visualizeTracksConfirm")
    })
  })
  
  
  #updates on UI elements ----------
  updateClusterTab <- function()
  {
    cluster_df <- as.data.frame(table(seurat_object$seurat_clusters)) #as.data.frame(table(Idents(seurat_object)))
    colnames(cluster_df)[1] <- "Cluster"
    colnames(cluster_df)[2] <- "Number of cells"
    cluster_df$`% of cells per cluster` <-  as.double(sprintf("%.3f", as.double( cluster_df$`Number of cells`/length(seurat_object@meta.data$orig.ident)*100 )))
    output$clusterTable <- renderDataTable(cluster_df, options = list(pageLength = 10), rownames = F)
    export_clustertable_RNA <<- cluster_df
  }
  
  updateReduction <- function()
  {
    session$sendCustomMessage("handler_disableAllButtons", "umapConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first and then re-run UMAP or tSNE or Diffusion Map above.")
      else {
        showModal(modalDialog(div('Analysis in Progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("umapPlot_loader")
        
        #get input
        dims <- as.numeric(input$umapDimensions)
        type <- input$umapType
        
        #prepare metadata
        meta <- seurat_object@meta.data
        meta$Cell_id <- rownames(meta)
        meta <- meta[, ]
        reduc_data <- data.frame()
        
        #prepare colors
        cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy])))
        
        #for all reductions
        seurat_object_reduc <- as.data.frame(seurat_object@reductions[[input$umapType]]@cell.embeddings)
        seurat_object_reduc <- seurat_object_reduc[, c(1:ncol(seurat_object_reduc))]
        seurat_object_reduc$Cell_id <- rownames(seurat_object_reduc)
        reduc_data <- left_join(seurat_object_reduc, meta)
        
        if(type == "umap" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="UMAP_1", y="UMAP_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder))+
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="UMAP 1", y="UMAP 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})
        }
        else if(type == "umap" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~UMAP_1, y=~UMAP_2, z=~UMAP_3, type="scatter3d", alpha = as.numeric(input$umapDotOpacity), mode="markers", color=as.formula(paste0('~', input$umapColorBy)),
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) )
          
          output$umapPlot <- renderPlotly({print(p)})
        }
        else if(type == "tsne" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="tSNE_1", y="tSNE_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="tSNE 1", y="tSNE 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "tsne" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~tSNE_1, y=~tSNE_2, z=~tSNE_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacity), color=as.formula(paste0('~', input$umapColorBy)), 
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) ) 
          output$umapPlot <- renderPlotly({print(p)})
        }
        else if(type == "dfm" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="DC_1", y="DC_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="DC 1", y="DC 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "dfm" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~DC_1, y=~DC_2, z=~DC_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacity), color=as.formula(paste0('~', input$umapColorBy)), 
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) ) 
          output$umapPlot <- renderPlotly({print(p)})
        }
        else if(type == "pca" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="PC_1", y="PC_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="PC 1", y="PC 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "pca" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~PC_1, y=~PC_2, z=~PC_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacity), color=as.formula(paste0('~', input$umapColorBy)), 
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) ) 
          output$umapPlot <- renderPlotly({print(p)})
        }
        else if(type == "phate" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="PHATE_1", y="PHATE_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="PHATE 1", y="PHATE 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "phate" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~PHATE_1, y=~PHATE_2, z=~PHATE_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacity), color=as.formula(paste0('~', input$umapColorBy)), 
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) ) 
          output$umapPlot <- renderPlotly({print(p)})
        }
        else
        {
          session$sendCustomMessage("handler_alert", "There was an error with the number of dimensions of the plot. Please set Dimensions to 2D and update the plot.")
        }
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with drawing the resutls.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "umapConfirm")
    })
  }
  
  updateReduction2 <- function(typeR)
  {
    session$sendCustomMessage("handler_disableAllButtons", "umapConfirm")
    tryCatch({
      if (identical(seurat_object, NULL)) session$sendCustomMessage("handler_alert", "Please, upload some data via the DATA INPUT tab first.")
      else if (identical(seurat_object@meta.data$seurat_clusters, NULL)) session$sendCustomMessage("handler_alert", "Please, execute CLUSTERING first and then re-run UMAP or tSNE or Diffusion Map above.")
      else {
        showModal(modalDialog(div('Analysis in Progress. This operation may take several minutes, please wait...', style='color:#ffffff; background-color:#222d32;'), footer = NULL, style = 'font-size:20px; text-align:center;position:absolute;')) #position:absolute;top:50%;left:50%
        shinyjs::show("umapPlot_loader")
        
        #get input
        dims <- as.numeric(input$umapDimensions)
        type <- typeR
        
        #prepare metadata
        meta <- seurat_object@meta.data
        meta$Cell_id <- rownames(meta)
        meta <- meta[, ]
        reduc_data <- data.frame()
        
        #prepare colors
        cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy])))
        
        #for all reductions
        seurat_object_reduc <- as.data.frame(seurat_object@reductions[[type]]@cell.embeddings)
        seurat_object_reduc <- seurat_object_reduc[, c(1:ncol(seurat_object_reduc))]
        seurat_object_reduc$Cell_id <- rownames(seurat_object_reduc)
        reduc_data <- left_join(seurat_object_reduc, meta)
        print(head(reduc_data)) 
        
        if(type == "umap" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="UMAP_1", y="UMAP_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder))+
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="UMAP 1", y="UMAP 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})
        }
        else if(type == "umap" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~UMAP_1, y=~UMAP_2, z=~UMAP_3, type="scatter3d", alpha = as.numeric(input$umapDotOpacity), mode="markers", color=as.formula(paste0('~', input$umapColorBy)),
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) )
          
          output$umapPlot <- renderPlotly({print(p)})
        }
        else if(type == "tsne" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="tSNE_1", y="tSNE_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="tSNE 1", y="tSNE 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "tsne" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~tSNE_1, y=~tSNE_2, z=~tSNE_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacity), color=as.formula(paste0('~', input$umapColorBy)), 
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) ) 
          output$umapPlot <- renderPlotly({print(p)})
        }
        else if(type == "dfm" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="DC_1", y="DC_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="DC 1", y="DC 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "dfm" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~DC_1, y=~DC_2, z=~DC_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacity), color=as.formula(paste0('~', input$umapColorBy)), 
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) ) 
          output$umapPlot <- renderPlotly({print(p)})
        }
        else if(type == "pca" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="PC_1", y="PC_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="PC 1", y="PC 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "pca" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~PC_1, y=~PC_2, z=~PC_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacity), color=as.formula(paste0('~', input$umapColorBy)), 
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) ) 
          output$umapPlot <- renderPlotly({print(p)})
        }
        else if(type == "phate" & dims == 2 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- ggplot(data=reduc_data, aes_string(x="PHATE_1", y="PHATE_2", fill=input$umapColorBy)) +
            geom_point(size= as.numeric(input$umapDotSize), shape=21, alpha= as.numeric(input$umapDotOpacity), stroke=as.numeric(input$umapDotBorder)) +
            scale_fill_manual(values = cols)+
            scale_size()+
            theme_bw() +
            theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                  axis.title.y = element_text(face = "bold", color = "black", size = 25),
                  axis.title.x = element_text(face = "bold", color = "black", size = 25),
                  legend.text = element_text(face = "bold", color = "black", size = 9),
                  legend.title = element_text(face = "bold", color = "black", size = 9),
                  legend.position="right",
                  title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
            labs(x="PHATE 1", y="PHATE 2", color="Cell type", title = "", fill="Color")
          output$umapPlot <- renderPlotly({plotly::ggplotly(p)})  
        }
        else if(type == "phate" & dims == 3 & dims <= ncol(seurat_object_reduc)-1)
        {
          p <- plot_ly(reduc_data, x=~PHATE_1, y=~PHATE_2, z=~PHATE_3, type="scatter3d", mode="markers", alpha = as.numeric(input$umapDotOpacity), color=as.formula(paste0('~', input$umapColorBy)), 
                       marker = list(size = as.numeric(input$umapDotSize), 
                                     line = list(color = 'black', width = as.numeric(input$umapDotBorder))
                       ),
                       colors = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(meta[, input$umapColorBy]))) ) 
          output$umapPlot <- renderPlotly({print(p)})
        }
        else
        {
          session$sendCustomMessage("handler_alert", "There was an error with the number of dimensions of the plot. Please set Dimensions to 2D and update the plot.")
        }
      }
    }, error = function(e) {
      print(paste("Error :  ", e))
      session$sendCustomMessage("handler_alert", "There was an error with drawing the resutls.")
    }, finally = {
      removeModal()
      Sys.sleep(1)
      session$sendCustomMessage("handler_enableAllButtons", "umapConfirm")
    })
  }
  
  updateFeaturePair <- function()
  {
      if (!identical(seurat_object, NULL) & input$findMarkersFeaturePairReductionType != "-" & input$findMarkersFeaturePair1 %in% rownames(seurat_object) & input$findMarkersFeaturePair2 %in% rownames(seurat_object))
      {
        geneS1 <- input$findMarkersFeaturePair1
        geneS2 <- input$findMarkersFeaturePair2
        label_x <- ""
        label_y <- ""
        show_label <- as.logical(input$findMarkersFeaturePairLabels)
        order_exp <- as.logical(input$findMarkersFeaturePairOrder)
        minq <- paste0("q", input$findMarkersFeaturePairMinCutoff)
        maxq <- paste0("q", input$findMarkersFeaturePairMaxCutoff)
        blendThr <- as.numeric(input$findMarkersBlendThreshold)
        
        if(input$findMarkersFeaturePairReductionType == "umap")
        {
          label_x <- "UMAP_1"
          label_y <- "UMAP_2"
        }
        else if(input$findMarkersFeaturePairReductionType == "tsne")
        {
          label_x <- "tSNE_1"
          label_y <- "tSNE_2"
        }
        else if(input$findMarkersFeaturePairReductionType == "dfm")
        {
          label_x <- "DC_1"
          label_y <- "DC_2"
        }
        else if(input$findMarkersFeaturePairReductionType == "pca")
        {
          label_x <- "PC_1"
          label_y <- "PC_2"
        }
        else if(input$findMarkersFeaturePairReductionType == "phate")
        {
          label_x <- "PHATE_1"
          label_y <- "PHATE_2"
        }
        
        plot_temp <- FeaturePlot(seurat_object, features = c(geneS1, geneS2), blend.threshold = blendThr, 
                            pt.size = 1.5, label = show_label, label.size = 5, cols = c("lightgrey", "red", "dodgerblue4"), 
                            order = order_exp, reduction = input$findMarkersFeaturePairReductionType, blend = TRUE, max.cutoff = maxq, min.cutoff = minq)
        
        plot <- FeaturePlot(seurat_object, features = c(geneS1, geneS2), blend.threshold = blendThr, 
                            pt.size = 1.5, label = show_label, label.size = 5, cols = c("lightgrey", "red", "dodgerblue4"), raster = F,
                            order = order_exp, reduction = input$findMarkersFeaturePairReductionType, blend = TRUE, max.cutoff = maxq, min.cutoff = minq) +
          theme_bw() +
          theme(axis.text.x = element_text(face = "bold", color = "black", size = 25, angle = 0),
                axis.text.y = element_text(face = "bold", color = "black", size = 25, angle = 0),
                axis.title.y = element_text(face = "bold", color = "black", size = 25),
                axis.title.x = element_text(face = "bold", color = "black", size = 25),
                legend.text = element_text(face = "bold", color = "black", size = 9),
                legend.title = element_text(face = "bold", color = "black", size = 9),
                legend.position="none",
                title = element_text(face = "bold", color = "black", size = 25, angle = 0)) +
        labs(x=label_x, y=label_y, color="Normalized\nexpression")
        gp1 <- plotly::ggplotly(plot[[1]] + xlim(min(plot_temp[[1]]$data[label_x]), max(plot_temp[[1]]$data[label_x])) + ylim(min(plot_temp[[1]]$data[label_y]), max(plot_temp[[1]]$data[label_y])) )
        gp2 <- plotly::ggplotly(plot[[2]] + xlim(min(plot_temp[[1]]$data[label_x]), max(plot_temp[[1]]$data[label_x])) + ylim(min(plot_temp[[1]]$data[label_y]), max(plot_temp[[1]]$data[label_y])) )
        gp3 <- plotly::ggplotly(plot[[3]] + xlim(min(plot_temp[[1]]$data[label_x]), max(plot_temp[[1]]$data[label_x])) + ylim(min(plot_temp[[1]]$data[label_y]), max(plot_temp[[1]]$data[label_y])) )
        gp4 <- plotly::ggplotly(plot[[4]]+theme(title = element_text(face = "bold", color = "black", size = 15)))
        
        output$findMarkersFPfeature1 <- renderPlotly({ gp1 })
        output$findMarkersFPfeature2 <- renderPlotly({ gp2 })
        output$findMarkersFPfeature1_2 <- renderPlotly({ gp3 })
        output$findMarkersFPcolorbox <- renderPlotly({ gp4 })  
      }
  }
  
  updateUtilitiesAssays <- function()
  {
    getAssays <- names(seurat_object@assays)
    updateSelectInput(session, "utilitiesActiveAssay", choices = getAssays, selected = DefaultAssay(seurat_object))
  }
  
  updateQC_choices <- function()
  {
    updateSliderInput(session, "minUniqueGenes", min = min(init_seurat_object$nFeature_RNA)-1, max = max(init_seurat_object$nFeature_RNA)-2, value = as.numeric(quantile(init_seurat_object$nFeature_RNA, probs = 0.02)))
    updateSliderInput(session, "maxUniqueGenes", min = min(init_seurat_object$nFeature_RNA)+1, max = max(init_seurat_object$nFeature_RNA)+1, value = as.numeric(quantile(init_seurat_object$nFeature_RNA, probs = 0.98)))
  }
  
  updateMaxHVGs <- function()
  {
    if(length(rownames(seurat_object)) < 8000)
      updateSliderInput(session = session, inputId = "nHVGs", min = 200, max = length(rownames(seurat_object)), value = min(2000, length(rownames(seurat_object))))
  }
  
  updatePCs_selection <- function(numPCs)
  {
    updateSliderInput(session = session, inputId = "snnPCs", value = numPCs) 
    updateSliderInput(session = session, inputId = "umapPCs", value = numPCs)
    updateSliderInput(session = session, inputId = "doubletsPCs", value = numPCs)
  }
  
  updateUmapTypeChoices <- function(type)
  {
    reductions_choices <<- c(reductions_choices, type)
    updateSelectInput(session, "umapType", choices = reductions_choices, selected = type)  
    updateSelectInput(session, "findMarkersReductionType", choices = reductions_choices, selected = type)
    updateSelectInput(session, "findMarkersFeaturePairReductionType", choices = reductions_choices, selected = type)
    updateSelectInput(session, "doubletsReduction", choices = reductions_choices, selected = type)
    updateSelectInput(session, "cellCycleReduction", choices = reductions_choices, selected = type)
    updateSelectInput(session, "trajectoryReduction", choices = reductions_choices, selected = type)
  }
  
  #function update selectInput
  updateSelInpColor <- function()
  {
    colnames(seurat_object@meta.data)
    tableMeta <- seurat_object@meta.data
    f <- sapply(tableMeta, is.factor)
    factors <- colnames(tableMeta[, f])
    updateSelectInput(session, "umapColorBy", choices = factors)
    updateSelectInput(session, "qcColorBy", choices = factors)
    updateSelectInput(session, "clusterGroupBy", choices = factors)
    updateSelectInput(session, "utilitiesActiveClusters", choices = factors)
  }
  
  updateRegressOut <- function()
  {
    variables <- colnames(seurat_object@meta.data)
    updateSelectInput(session, "normalizeRegressColumns", choices = variables)
  }
  
  updateGeneSearchFP <- function()
  {
    total_genes <- rownames(seurat_object[[DefaultAssay(seurat_object)]]@counts)
    
    #get signature names or numeric columns
    colnames(seurat_object@meta.data)
    tableMeta <- seurat_object@meta.data
    f <- sapply(tableMeta, is.numeric)
    sig_names <- colnames(tableMeta[, f])
      
    updateSelectizeInput(session, 'findMarkersGeneSelect', choices = c(total_genes, sig_names), server = TRUE) # server-side selectize drastically improves performance
    updateSelectizeInput(session, 'findMarkersGeneSelect2', choices = c(total_genes, sig_names), server = TRUE)
    updateSelectizeInput(session, 'findMarkersFeaturePair1', choices = total_genes, server = TRUE)
    updateSelectizeInput(session, 'findMarkersFeaturePair2', choices = total_genes, server = TRUE)
  }
  
  updateInputLineageList <- function(lin_names)
  {
    updateSelectInput(session, "trajectoryLineageSelect", choices = lin_names)
  }
  
  updateInputLRclusters <- function() #update after clustering clusternames used in --> Volcano, gProfiler, slingshot, nichenetR
  {
    all_cluster_names <- (levels(seurat_object@meta.data[, 'seurat_clusters']))
    updateSelectInput(session, "ligandReceptorSender", choices = all_cluster_names)
    updateSelectInput(session, "ligandReceptorReciever", choices = all_cluster_names)
    updateSelectInput(session, "trajectoryStart", choices = all_cluster_names)
    updateSelectInput(session, "trajectoryEnd", choices = all_cluster_names)
    updateSelectInput(session, "findMarkersClusterSelect", choices = all_cluster_names)
    updateSelectInput(session, "gProfilerList", choices = all_cluster_names)
    updateSelectizeInput(session, "gProfilerFlameSelection", choices = all_cluster_names)
    updateSelectInput(session, "utilitiesDeleteCluster", choices = all_cluster_names)
    updateSelectInput(session, "utilitiesRenameOldName", choices = all_cluster_names)
  }
  
  updateFilteredQCplots <- function()
  {
    output$filteredNFeatureViolin <- renderPlotly(
      {
        p <- VlnPlot(seurat_object, features = c("nFeature_RNA"), pt.size = 0, group.by = "orig.ident",
                     cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, 'orig.ident'])))) + 
          theme_bw() + 
          geom_hline(yintercept=c(as.numeric(input$minUniqueGenes), as.numeric(input$maxUniqueGenes)), linetype="dashed", color = "red", size=1) + 
          theme(
            plot.title = element_blank(),
            axis.title.x = element_blank(),
            legend.position = "none") +
          labs(title = "", y="Genes detected/cell")
        plotly::ggplotly(p, tooltip = c("x", "y"))
      }
    )
    
    output$filteredTotalCountsViolin <- renderPlotly(
      {
        p <- VlnPlot(seurat_object, features = c("nCount_RNA"), pt.size = 0, group.by = "orig.ident",
                     cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, 'orig.ident'])))) + 
          theme_bw() + 
          theme(
            plot.title = element_blank(),
            axis.title.x = element_blank(),
            legend.position = "none")+
          labs(title = "", y="Total counts/cell")
        plotly::ggplotly(p, tooltip = c("x", "y")) 
      }
    )
    
    output$filteredMitoViolin <- renderPlotly(
      {
        p <- VlnPlot(seurat_object, features = c("percent.mt"), pt.size = 0, group.by = "orig.ident",
                     cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, 'orig.ident'])))) + 
          theme_bw() + 
          geom_hline(yintercept= as.numeric(input$maxMtReads), linetype="dashed", color = "red", size=1) +
          theme(
            plot.title = element_blank(),
            axis.title.x = element_blank(),
            legend.position = "none")+
          labs(title = "", y="% of reads mapped to mitochonrial genome/cell")
        plotly::ggplotly(p, tooltip = c("x", "y"))
      }
    )
    
    output$filteredMtCounts <- renderPlotly(
      {
        p <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "orig.ident", raster = F,
                            cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, 'orig.ident']))))
        gp <- plotly::ggplotly(p)
        print(gp)
      }
    )
    
    output$filteredGenesCounts <- renderPlotly(
      {
        p <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "orig.ident", raster = F,
                            cols = colorRampPalette(brewer.pal(12, "Paired"))(length(unique(seurat_object@meta.data[, 'orig.ident']))))
        gp <- plotly::ggplotly(p)
        print(gp)
      }
    )
    
    output$filteredCellStats <- renderPrint(
      {
        cat(paste0("Total number of cells: ", nrow(seurat_object@meta.data)))
      }
    )
  }
  
  # Helper Functions ####
  
  sendToAPI <- function(apiLink, jsonBody) {
    response <- httr::POST(apiLink, body = jsonBody, encode = "json")
    urlFromResponse <- httr::content(response, as = "parsed")$url
    session$sendCustomMessage("handler_browseUrl", urlFromResponse)
  }
  
  
  # This function is called after a new input file has been uploaded
  # and is responsible for clearing all generated plots across all tabs
  # @param fromDataInput: If TRUE, clears all, including QC, else skips clearing QC and metadata table
  cleanAllPlots <- function(fromDataInput){
    # renderPlotly
    if (fromDataInput){
      output$nFeatureViolin <- NULL
      output$totalCountsViolin <- NULL
      output$mitoViolin <- NULL
      output$mtCounts <- NULL
      output$genesCounts <- NULL
      output$cellStats <- NULL
      output$filteredNFeatureViolin <- NULL
      output$filteredTotalCountsViolin <- NULL
      output$filteredMitoViolin <- NULL
      output$filteredMtCounts <- NULL
      output$filteredGenesCounts <- NULL
      output$filteredCellStats <- NULL
      output$metadataTable <- NULL
      output$grnMatrixRNA <- NULL
      user_dir_pyscenic <<- ""
      user_dir <<- ""
    }
    output$elbowPlotPCA <- NULL
    output$PCAscatter <- NULL
    output$PCAloadings <- NULL
    output$PCAheatmap <- NULL
    output$gProfilerManhattan <- NULL
    output$annotateClustersCIPRDotplot <- NULL
    output$ligandReceptorFullHeatmap <- NULL
    output$ligandReceptorCuratedHeatmap <- NULL
    output$hvgScatter <- NULL
    output$hvgTop10Stats <- NULL
    output$cellCyclePCA <- NULL
    output$cellCycleBarplot <- NULL
    output$clusterBarplot <- NULL
    output$umapPlot <- NULL
    output$findMarkersHeatmap <- NULL
    output$findMarkersDotplot <- NULL
    output$findMarkersFeaturePlot <- NULL
    output$findMarkersFPcolorbox <- NULL
    output$findMarkersFPfeature1_2 <- NULL
    output$findMarkersFPfeature1 <- NULL
    output$findMarkersFPfeature2 <- NULL
    output$findMarkersViolinPlot <- NULL
    output$findMarkersVolcanoPlot <- NULL
    output$grnHeatmapRNA <- NULL
    output$grnHeatmapRNA_DecoupleR <- NULL
    
    # renderPlot
    output$trajectoryPlot <- NULL
    output$trajectoryPseudotimePlot <- NULL
    
    # renderDataTable
    output$clusterTable <- NULL
    output$gProfilerTable <- NULL
    output$annotateClustersCIPRTable <- NULL
    output$findMarkersTable <- NULL
    output$grnMatrixRNA_DecoupleR <- NULL
    output$grnMatrixRNA <- NULL
    
    # renderPrint
    if (fromDataInput) output$cellStats <- NULL
    output$trajectoryText <- NULL
    
    # dittoDimPlot
    plot_D <- NULL
    plot_P <- NULL
    
    reductions_choices <<- c("-")
    updateSelectInput(session, "umapType", choices = reductions_choices)
    updateSelectInput(session, "findMarkersReductionType", choices = reductions_choices)
    updateSelectInput(session, "findMarkersFeaturePairReductionType", choices = reductions_choices)
    updateSelectInput(session, "doubletsReduction", choices = reductions_choices)
    updateSelectInput(session, "cellCycleReduction", choices = reductions_choices)
    updateSelectInput(session, "trajectoryReduction", choices = reductions_choices)
    
    #reset export table values
    export_metadata_RNA <<- ""
    export_loadingScoresTable_RNA <<- ""
    export_clustertable_RNA <<- ""
    export_markerGenes_RNA <<- ""
    export_enrichedTerms_RNA <<- ""
    export_annotation_RNA <<- ""
    export_ligandReceptor_full_RNA <<- ""
    export_ligandReceptor_short_RNA <<- ""
    export_scenicAUC_full_RNA <<- ""
    export_decoupleRZscores_full_RNA <- ""
  }
  
  cleanAllPlotsATAC <- function()
  {
    #renderPlotly
    output$TSS_plot <- NULL
    output$TSS_nFrag_plot <- NULL
    output$nFrag_plot <- NULL
    output$CellStatsATAC <- NULL
    output$lsiOutput <- NULL
    output$clusterBarplotATAC <- NULL
    output$umapPlotATAC <- NULL
    output$findMarkersGenesHeatmapATAC <- NULL
    output$findMarkersPeaksHeatmapATAC <- NULL
    output$findMotifsHeatmapATAC <- NULL
    output$trajectoryTextATAC <- NULL
    output$grnHeatmapATAC <- NULL
    
    #renderPlot
    output$findMarkersFeaturePlotATAC <- NULL
    output$trajectoryPseudotimePlotATAC <- NULL
    
    #datatable	
    output$clusterTableATAC <- NULL
    output$findMarkersGenesTableATAC <- NULL
    output$findMarkersPeaksTableATAC <- NULL
    output$findMotifsTableATAC <- NULL
    output$grnMatrixATAC <- NULL
    output$grnP2GlinksTable <- NULL
    output$metadataTableATAC <- NULL
    output$grnMotifTable <- NULL
    
    #export tables
    export_metadata_ATAC <- ""
    export_clustertable_ATAC <- ""
    export_markerGenes_ATAC <- ""
    export_markerPeaks_ATAC <- ""
    export_motifs_ATAC <- ""
    export_positiveRegulators_ATAC <- ""
    export_peakToGenelinks_ATAC <- ""
    export_PeakMotifTable_ATAC <- ""
  }
  
  #update metadata RNA and export_RNA_table
  updateMetadata <- function()
  {
    seurat_object@meta.data$Cell_id <<- rownames(seurat_object@meta.data)
    #seurat_object@meta.data$percent.mt <<- as.double(sprintf("%.3f", as.double(seurat_object@meta.data$percent.mt))) 
    
    #--3 decimal points
    for(i in 1:length(colnames(seurat_object@meta.data)))
    {
      if(is.double(seurat_object@meta.data[, i]))
      {
        seurat_object@meta.data[, i] <- as.double(sprintf("%.3f", as.double(seurat_object@meta.data[, i]))) 
      }
    }
    #--
    
    output$metadataTable <- renderDataTable(seurat_object@meta.data, options = list(pageLength = 10), rownames = F)
    export_metadata_RNA <<- seurat_object@meta.data
    updateRegressOut()
  }
  
  updateMetadataATAC <- function()
  {
    df_meta_data <- as.data.frame(getCellColData(proj_default))
    df_meta_data$cell_id <- rownames(df_meta_data)
    rownames(df_meta_data) <- NULL
    
    #--3 decimal points
    for(i in 1:length(colnames(df_meta_data)))
    {
      if(is.double(df_meta_data[, i]))
      {
        df_meta_data[, i] <- as.double(sprintf("%.3f", as.double(df_meta_data[, i]))) 
      }
    }
    #--
    
    output$metadataTableATAC <- renderDataTable(df_meta_data, options = list(pageLength = 10), rownames = F)
    
    #export table
    export_metadata_ATAC <<- df_meta_data
  }
  
  #functions for optimal number of PCs
  #Initialize Functions
  PrepDR <- function( # From Seurat
    object,
    genes.use = NULL,
    use.imputed = FALSE,
    assay.type="RNA"
  ) {
    
    if (length(VariableFeatures(object = object)) == 0 && is.null(x = genes.use)) {
      stop("Variable genes haven't been set. Run MeanVarPlot() or provide a vector
         of genes names in genes.use and retry.")
    }
    if (use.imputed) {
      data.use <- t(x = scale(x = t(x = object@imputed)))
    } else {
      data.use <- GetAssayData(object, assay.type = assay.type,slot = "scale.data")
    }
    genes.use <- if(is.null(genes.use)) VariableFeatures(object = object) else genes.use 
    genes.use <- unique(x = genes.use[genes.use %in% rownames(x = data.use)])
    genes.var <- apply(X = data.use[genes.use, ], MARGIN = 1, FUN = var)
    genes.use <- genes.use[genes.var > 0]
    genes.use <- genes.use[! is.na(x = genes.use)]
    data.use <- data.use[genes.use, ]
    return(data.use)
  }
  
  PCA_estimate_nPC<-function(data, k=10, from.nPC = 2, to.nPC=150, by.nPC=5, maxit=200, seed=617) {
    PC <-seq(from = from.nPC, to = to.nPC, by = by.nPC)
    # Init the error matrices
    error1<-matrix(0, nrow = length(c(1:k)), ncol = length(PC))
    error2<-matrix(0, nrow = length(c(1:k)), ncol = length(PC))
    print(paste0(k,"-fold paritioning..."))
    # K-Fold Partitioning
    dgem.kfold<-dismo::kfold(t(data), k=k)
    # SVD-CV based on https://stats.stackexchange.com/questions/93845/how-to-perform-leave-one-out-cross-validation-for-pca-to-determine-the-number-of
    for(i in c(1:k)) {
      print(paste0("k:",i))
      X.train<-t(data[, dgem.kfold!=i])
      X.test<-t(data[, dgem.kfold==i])
      # Find a few approximate singular values and corresponding singular vectors of a matrix.
      print("Running SVD...")
      # Seurat uses IRLBA to do PCA : https://github.com/satijalab/seurat/blob/cec7cb95c73fd6d605723e9af9a1f96eda5635de/R/dimensional_reduction.R
      pca.results<-irlba::irlba(A = X.train, nv = to.nPC, maxit = maxit) 
      gl<-pca.results$v
      for(j in 1:length(PC)) {
        print(paste0("Ndims:",PC[j]))
        P<-gl[,1:PC[j]]%*%t(gl[,1:PC[j]])
        # Naive method
        err1<-X.test %*% (diag(dim(P)[1]) - P)
        # Approximate method
        err2<-X.test %*% (diag(dim(P)[1]) - P + diag(diag(P)))
        error1[i,j]<-sum(err1^2)
        error2[i,j]<-sum(err2^2)
        rm(err1)
        rm(err2)
      }
    }
    errors1<-colSums(error1)
    errors2<-colSums(error2)
    nPC=PC[which(errors2 == min(errors2))]
    plot(PC,errors1)
    plot(PC,errors2)
    return(nPC)
  }
  
  ###disable Tabs
  
  disableTabsRNA <- function()
  {
    hideTab(inputId="uploadTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId="qcTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId="pcaTabPanel", target="scRNA-seq: PCA", session = session)
    hideTab(inputId="clusteringTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId="umapTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId="findMarkersTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId="featuresTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId ="doubletDetectionTabPanel", target = "scRNA-seq", session = session)
    hideTab(inputId="gProfilerTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId="annotateClustersTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId="trajectoryTabPanel", target="scRNA-seq", session = session)
    hideTab(inputId="grnTabPanel", target="scRNA-seq", session = session)
    
    
    ############################ show ATAC and make it selected
    
    showTab(inputId="uploadTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="qcTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="pcaTabPanel", target="scATAC-seq: LSI", session = session)
    showTab(inputId="clusteringTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="umapTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="findMarkersTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="featuresTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="doubletDetectionTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="gProfilerTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="annotateClustersTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="trajectoryTabPanel", target="scATAC-seq", session = session)
    showTab(inputId="grnTabPanel", target="scATAC-seq", session = session)
    
    updateTabsetPanel(inputId="uploadTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="qcTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="pcaTabPanel", selected="scATAC-seq: LSI", session = session)
    updateTabsetPanel(inputId="clusteringTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="umapTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="findMarkersTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="featuresTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="doubletDetectionTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId = "ATAC_markers_tabs", selected = "Marker genes (ATAC)", session = session)
    updateTabsetPanel(inputId="gProfilerTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="annotateClustersTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="trajectoryTabPanel", selected="scATAC-seq", session = session)
    updateTabsetPanel(inputId="grnTabPanel", selected="scATAC-seq", session = session)
  }
  
  disableTabsATAC <- function()
  {
    hideTab(inputId="uploadTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId="qcTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId="pcaTabPanel", target="scATAC-seq: LSI", session = session)
    hideTab(inputId="clusteringTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId="umapTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId="findMarkersTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId="featuresTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId ="doubletDetectionTabPanel", target = "scATAC-seq", session = session)
    hideTab(inputId="gProfilerTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId="annotateClustersTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId="trajectoryTabPanel", target="scATAC-seq", session = session)
    hideTab(inputId="grnTabPanel", target="scATAC-seq", session = session)
    
    ############################ show RNA and make it selected
    
    showTab(inputId="uploadTabPanel", target="scRNA-seq", session = session)
    showTab(inputId="qcTabPanel", target="scRNA-seq", session = session)
    showTab(inputId="pcaTabPanel", target="scRNA-seq: PCA", session = session)
    showTab(inputId="clusteringTabPanel", target="scRNA-seq", session = session)
    showTab(inputId="umapTabPanel", target="scRNA-seq", session = session)
    showTab(inputId="findMarkersTabPanel", target="scRNA-seq", session = session)
    showTab(inputId="featuresTabPanel", target="scRNA-seq", session = session)
    showTab(inputId = "doubletDetectionTabPanel", target = "scRNA-seq", session = session)
    showTab(inputId="gProfilerTabPanel", target="scRNA-seq", session = session)
    showTab(inputId="annotateClustersTabPanel", target="scRNA-seq", session = session)
    showTab(inputId="trajectoryTabPanel", target="scRNA-seq", session = session)
    showTab(inputId="grnTabPanel", target="scRNA-seq", session = session)
    
    updateTabsetPanel(inputId="uploadTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="qcTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(session, inputId = "qc_tabs_rna", selected = "Pre-filtering plots")
    updateTabsetPanel(inputId="pcaTabPanel", selected="scRNA-seq: PCA", session = session)
    updateTabsetPanel(inputId="clusteringTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="umapTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="findMarkersTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="featuresTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="doubletDetectionTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="gProfilerTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="annotateClustersTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="trajectoryTabPanel", selected="scRNA-seq", session = session)
    updateTabsetPanel(inputId="grnTabPanel", selected="scRNA-seq", session = session)
  }
  
  
  #ATAC
  #function update selectInput
  updateSelInpColorATAC <- function()
  {
    tableMeta <- as.data.frame(getCellColData(proj_default))
    
    fac <- sapply(tableMeta, is.factor)
    
    chars <- sapply(tableMeta, is.character)
    
    selected <- (fac | chars)
    
    factors_and_chars <- colnames(tableMeta)[which(selected==T)]
    print(factors_and_chars)
    updateSelectInput(session, "umapColorByATAC", choices = factors_and_chars)
  }
  
  #update trajectory clusters (to be deleted)
  #updateInpuTrajectoryClustersATAC <- function()
  #{
  #  cluster_names <- unique(proj_default$Clusters)
  #  updateSelectInput(session, "trajectoryStartATAC", choices = cluster_names)
  #  updateSelectInput(session, "trajectoryEndATAC", choices = cluster_names)
  #}
  
  #alert
  observeEvent(input$removeProject, {
    #remove RNAseq project
    init_seurat_object <<- NULL
    seurat_object <<- NULL
    cleanAllPlots(TRUE)
    
    #remove ATAC-seq project
    proj_default <<- NULL
    cleanAllPlotsATAC()
    
    hideAllLoaders()
    Sys.sleep(1)
    shinyalert(title = "Project removed succesfully! You can now upload a new project.", type = "success")
    removeModal()
  })
  
  observeEvent(input$cancelProjectRemoval, {
    removeModal()
  })
  
  modal_confirm <- modalDialog(
    "This operation cannot be completed, because another project is already loaded. Do you want to discard it?",
    title = "New project",
    footer = tagList(
      actionButton(inputId="removeProject", label="Yes"),
      actionButton("cancelProjectRemoval", "No")
    )
  )
  
  
  #alert2 rename cluster
  observeEvent(input$confirmUtilitiesRename, {
    old_idents_for_change <- seurat_object$seurat_clusters
    seurat_object$seurat_clusters <<- gsub(pattern = paste0("^",input$utilitiesRenameOldName,"$"), replacement = input$utilitiesRenameNewName, x = old_idents_for_change)
    seurat_object$seurat_clusters <<- as.factor(seurat_object$seurat_clusters)
    Idents(seurat_object) <<- seurat_object$seurat_clusters
    updateInputLRclusters()
    updateMetadata()
    #showNotification("Successful operation.", type = "message", closeButton = T)
    shinyalert(title = "Cluster rename operation completed succesfully!", type = "success")
    removeModal()
  })
  
  observeEvent(input$cancelUtilitiesRename, {
    removeModal()
  })
  
  modal_confirm_utilities_rename <- modalDialog(
    "By renaming a cluster, you will need to update existing plots that have already been produced with the previous clustering scheme. 
    By changing the name of a cluster to an already existing cluster name will lead to merging of the two clusters. 
    Do you want to proceed?",
    title = "Cluster rename",
    footer = tagList(
      actionButton(inputId="confirmUtilitiesRename", label="Yes"),
      actionButton("cancelUtilitiesRename", "No")
    )
  )
  
  #alert3 delete cluster
  observeEvent(input$confirmUtilitiesDelete, {
    ident_to_remove <- input$utilitiesDeleteCluster
    seurat_object <<- subset(seurat_object, idents = ident_to_remove, invert = TRUE)
    seurat_object$seurat_clusters <<- Idents(seurat_object)
    updateInputLRclusters()
    updateMetadata()
    #showNotification("Successful operation.", type = "message", closeButton = T)
    shinyalert(title = "Cluster deletion operation completed succesfully!", type = "success")
    removeModal()
  })
  
  observeEvent(input$cancelUtilitiesDelete, {
    removeModal()
  })
  
  modal_confirm_utilities_delete <- modalDialog(
    "By deleting a cluster, the user is advised to repeat various steps of the analysis and/or update existing plots. 
    In more detail, the scaled expression values calculated in the tab \"DATA NORMALIZATION & SCALING\" should be recalculated as the
    final number of cell has changed.
    Do you want to proceed?",
    title = "Cluster deletion",
    footer = tagList(
      actionButton(inputId="confirmUtilitiesDelete", label="Yes"),
      actionButton("cancelUtilitiesDelete", "No")
    )
  )
}