diff --git a/Rfunctions/CompPlot.R b/Rfunctions/CompPlot.R index f049fd2..6bb5214 100644 --- a/Rfunctions/CompPlot.R +++ b/Rfunctions/CompPlot.R @@ -194,7 +194,7 @@ Plot_Comp_Logit <- function(input, BaseContrast, resDiff, SelectTaxoPlotCompDebo axis.text = element_text(size = input$axisFontSizeLogit), legend.title = element_text(size = input$legendFontSizeLogit), legend.text = element_text(size = input$legendFontSizeLogit)) - plot <- plot + geom_text(aes(label = labels, color = color_var), size = input$labelsSizeLogit %/% 3, hjust = 1, vjust = 1, show.legend = FALSE) + plot <- plot + geom_text_repel(aes(label = labels, color = color_var), size = input$labelsSize %/% 3, hjust = 1, vjust = 1, show.legend = FALSE) plot <- plot + scale_color_manual(name = "Legend title", values = c("Not significant" = input$colour01, "Significant for contrast 1" = input$colour02, "Significant for contrast 2" = input$colour03, "Significant for both contrasts" = input$colour04)) if (input$fixed11) { diff --git a/Rfunctions/VisuPlot.R b/Rfunctions/VisuPlot.R index e816239..12d0f8e 100644 --- a/Rfunctions/VisuPlot.R +++ b/Rfunctions/VisuPlot.R @@ -72,7 +72,6 @@ Plot_Visu_Barplot <- ## All the modalities for all the var of interest val = c(val, eval(expr)) } - ## Create the others data frame for the plot function all_dataBarPlot_mat = c() all_tmp_mat = matrix(0, ncol = 3, nrow = length(allnamesTax)) @@ -1504,120 +1503,6 @@ Plot_Visu_Tree <- function(input, resDiff, CT_Norm_OTU, taxo_table) ## NETWORK ## ## -#Function that computes a permutation test and return an adjacency matrix -#input: mat, max iteration number -#output: adjacency matrix - -# compute_pcor <- function(input, mat, n_iter = 100) { -# -# n_rows <- nrow(mat) -# correlation_storage <- vector("list", length = n_rows) -# for (i in 1:n_rows) { -# correlation_storage[[i]] <- numeric(n_iter * (n_rows - 1)) -# } -# adjacency_matrix <- matrix(0, nrow = n_rows, ncol = n_rows) -# correlation_matrix <- matrix(NA, nrow = n_rows, ncol = n_rows) -# -# # Loop through each pair of rows -# for (i in 1:(n_rows - 1)) { -# for (j in (i + 1):n_rows) { -# # Initialize vectors to store correlations between rows i and j for each shuffle -# correlations_i <- numeric(n_iter) -# correlations_j <- numeric(n_iter) -# -# # Shuffle and compute correlations n_iter times -# for (shuffle in 1:n_iter) { -# # Shuffle only the i-th and j-th row -# mat_shuffled_i <- mat -# mat_shuffled_j <- mat -# mat_shuffled_i[i, ] <- sample(mat[i, ]) -# mat_shuffled_j[j, ] <- sample(mat[j, ]) -# -# # Compute correlation of the shuffled i-th row with row j and vice versa -# correlations_i[shuffle] <- cor(mat_shuffled_i[i, ], mat[j, ]) -# correlations_j[shuffle] <- cor(mat_shuffled_j[j, ], mat[i, ]) -# } -# -# # Store the shuffled correlations -# correlation_storage[[i]][((j - 1) * n_iter + 1):(j * n_iter)] <- correlations_i -# correlation_storage[[j]][((i - 1) * n_iter + 1):(i * n_iter)] <- correlations_j -# -# # Compute the observed correlation once for each unique pair of rows -# observed_correlation <- cor(mat[i, ], mat[j, ]) -# correlation_matrix[i, j] <- observed_correlation -# correlation_matrix[j, i] <- observed_correlation # symmetry -# -# # Calculate p-values and fill the adjacency matrix using the observed_correlation -# p_value_i <- min(sum(correlations_i < observed_correlation), sum(correlations_i > observed_correlation)) / n_iter -# p_value_j <- min(sum(correlations_j < observed_correlation), sum(correlations_j > observed_correlation)) / n_iter -# -# # Determine significance based on the p-value and the sign of the observed correlation -# if ((p_value_i <= 0.025 || p_value_i >= 0.975) && observed_correlation > 0) { -# adjacency_matrix[i, j] <- 1 -# adjacency_matrix[j, i] <- 1 -# } else if ((p_value_i <= 0.025 || p_value_i >= 0.975) && observed_correlation < 0) { -# adjacency_matrix[i, j] <- -1 -# adjacency_matrix[j, i] <- -1 -# } -# } -# } -# -# # Return the adjacency matrix and the correlation matrix -# return(list(adjacency_matrix = adjacency_matrix, correlation_matrix = correlation_matrix)) -# } - - -# Register the parallel backend -# numCores <- detectCores() - 1 -# doParallel::registerDoParallel(cores=numCores) -# # Function to compute partial correlations -# compute_partial_correlations <- function(mat, n_iter, method, threshold) { -# n_rows <- nrow(mat) -# corr_storage <- vector("list", length = n_rows) -# adjacency <- matrix(0, nrow = n_rows, ncol = n_rows) -# corr_matrix <- matrix(NA, nrow = n_rows, ncol = n_rows) -# -# # Perform parallel computation -# results <- foreach(i = 1:(n_rows - 1), .combine = rbind, .multicombine = TRUE, .packages = c("stats")) %dopar% { -# local_corrs <- numeric(n_iter) -# for (shuffle in 1:n_iter) { -# mat_shuffled <- mat -# mat_shuffled[i, ] <- sample(mat[i, ]) -# local_corrs[shuffle] <- cor(mat_shuffled[i, ], mat[j, ], method = method) -# } -# -# # Compute the observed correlation for the i-th row -# observed_corr <- sapply((i + 1):n_rows, function(j) cor(mat[i, ], mat[j, ], method = method)) -# p_values <- sapply((i + 1):n_rows, function(j) { -# local_pval <- min(sum(local_corrs < observed_corr[j-i]), sum(local_corrs > observed_corr[j-i])) / n_iter -# local_pval -# }) -# -# cbind(i, (i + 1):n_rows, observed_corr, p_values) -# } -# -# # Process results -# for (result in results) { -# i <- result[, 1] -# j <- result[, 2] -# observed_corr <- result[, 3] -# p_value <- result[, 4] -# -# # Fill the adjacency matrix -# adjacency[i, j] <- ifelse(p_value <= threshold/2 || p_value >= 1 - threshold/2, sign(observed_corr), 0) -# corr_matrix[i, j] <- observed_corr -# } -# -# # Ensure symmetry -# adjacency <- adjacency + t(adjacency) -# corr_matrix <- corr_matrix + t(corr_matrix) -# -# list(adjacency = adjacency, corr_matrix = corr_matrix) -# } -# -# # Stop the parallel cluster when you are done -# doParallel::stopImplicitCluster() - Plot_network <- function(input, resDiff, @@ -1671,56 +1556,10 @@ Plot_network <- if (!is.null(counts_tmp_combined)) { ################################# - ###### partial correlation ###### + ###### correlation ###### ################################# countsMatrix <- as.matrix(counts_tmp_combined) - # ppcor <- - # ppcor::pcor(countsMatrix, method = input$pcorrMethod) - # pcor <- ppcor$estimate - # rownames(pcor) <- colnames(countsMatrix) - # colnames(pcor) <- colnames(countsMatrix) - # - # # Step 1: Calculate the variance-covariance matrix - # covMatrix <- cov(countsMatrix) - # - # # Step 2: Compute the inverse - # invCovMatrix <- MASS::ginv(covMatrix) - # - # - # # Step 3: Calculate partial correlations - # pcorMatrix <- -invCovMatrix / sqrt(outer(diag(invCovMatrix), diag(invCovMatrix))) - # - # # Set diagonal to zero (self-correlation is not meaningful) - # diag(pcorMatrix) <- 0 - # - # # Step 4: Calculate p-values - # t_values <- pcorMatrix * sqrt((n - 2) / (1 - pcorMatrix^2)) - # p_values <- 2 * pt(-abs(t_values), df = n - 2) - # - - ################################ - ####### adjacency matrix ####### - ################################ - - - #if (input$colorCorr == "pcorr") { - # adjacency <- matrix(apply(p_values, c(1, 2), function(p_val, pcor_val) { - # if (!is.na(p_val) && p_val <= 0.05) { - # return(1) - # } else { - # return(0) - # } - # }, pcor_val = pcor), nrow = nrow(p_values), ncol = ncol(p_values)) - - # adjacency <- matrix(apply(pcor, c(1, 2), function(x) { - # # Check for NA values in x and ensure input$pcorrThreshold is not NA - # if (!is.na(x) && !is.na(input$pcorrThreshold) && abs(x) > as.numeric(input$pcorrThreshold)) { - # return(x) - # } else { - # return(0) - # } - # }), nrow = nrow(pcor), ncol = ncol(pcor)) permutation <- compute_pcor adjacency <- permutation$adjacency @@ -1766,7 +1605,6 @@ Plot_network <- req(igraphGraph) dataVN <- toVisNetworkData(igraphGraph) dataVN$nodes$title <- paste0("", dataVN$nodes$id, "") - dataVN$nodes$label <- dataVN$nodes$id # dataVN$nodes$label <- sapply(dataVN$nodes$id, function(x)if(is.element(x, list_to_label)){x}else{""}) dataVN$nodes$cor <- 0 dataVN$edges$pcor <- 0 @@ -1920,7 +1758,7 @@ Plot_network <- dataVN$edges[which(!dataVN$edges$from %in% dataVN$nodes),] dataVN$edges[which(!dataVN$edges$to %in% dataVN$nodes),] - + dataVN$nodes$label <- paste0(dataVN$nodes$id, " - C", dataVN$nodes$community) plot <- visNetwork(nodes = dataVN$nodes, edges = dataVN$edges) plot <- diff --git a/server.R b/server.R index 0772944..f80de64 100644 --- a/server.R +++ b/server.R @@ -7604,8 +7604,17 @@ CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT" # resDiff = ResDiffAnal() # if(!is.null(resDiff$dds)) Plot_Visu_Diversity(input,resDiff,type="box") # }) + +observe({ + if(input$PlotVisuSelect == "Network") + updateRadioButtons(session, "SelectSpecifTaxo", choiceNames = c("Most abundant","All", "Differential features", "Non differential features", "By Network communities"), choiceValues = c("Most", "All", "Diff", "NoDiff", "ByNetwork")) + if(input$NormOrRaw == "vst") + updateRadioButtons(session, "SelectSpecifTaxo", choiceNames = c("Most abundant","All", "Differential features", "Non differential features"), choiceValues = c("Most", "All", "Diff", "NoDiff")) + +}) get_others <- reactive({ + data = dataInput()$data taxo = input$TaxoSelect resDiff = ResDiffAnal() @@ -7622,6 +7631,37 @@ CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT" Available_taxo = rownames(counts)[ord] selTaxo = Available_taxo[1:min(12, length(Available_taxo))] + if (input$SelectSpecifTaxo == "ByNetwork") + { + taxoCommunities = Plot_network( + input, + ResDiffAnal(), + Available_taxo, + compute_pcor = compute_pcor(), + SelectTaxoPlotNetworkDebounce(), + qualiVariable, + colors = colorsVisuPlot(), + dataInput = dataInput() + )$taxo + + req(input$CommunityVisuPlot) + selTaxo <- taxoCommunities %>% + dplyr::filter(Community == as.integer(input$CommunityVisuPlot)) %>% + dplyr::select(input$TaxoSelect) %>% + dplyr::distinct() + + selTaxo <- selTaxo[, 1] + res = selectizeInput( + "selectTaxoPlot", + h6(strong( + paste("Select the", input$TaxoSelect, "to plot") + )), + Available_taxo, + selected = selTaxo, + multiple = TRUE + ) + + } if (input$SelectSpecifTaxo == "Most") res = selectizeInput( "selectTaxoPlot", @@ -7761,7 +7801,7 @@ CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT" selected = Available_taxo, multiple = TRUE ) - others = setdiff(Available_taxo, selTaxo) + others = setdiff(Available_taxo, selTaxo) } return(list(res = res, others = others)) }) @@ -7859,17 +7899,37 @@ CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT" }) output$SelectValueQualiVar <- renderUI({ - target = isolate(values$TargetWorking) + VarInt = input$VisuVarInt + target = values$TargetWorking + val = NULL + if (length(VarInt) > 0) + { + + #Get the modalities + for (i in 1:length(VarInt)) + { + ## Replace "-" by "." + target[, VarInt[i]] = gsub("-", ".", target[, VarInt[i]]) + + Tinput = paste("input$", "ModVisu", VarInt[i], sep = "") + expr = parse(text = Tinput) + ## All the modalities for all the var of interest + val = c(val, eval(expr)) + } + } + req(val) if (!is.null(target)) { conditionalPanel( condition = "input.colorCorr == 'corr'", selectInput( "values1", 'Values to consider as "1"', - choices = as.character(unique(as.factor(target[, input$sec_variable]))), - selected = min(as.character(unique( - as.factor(target[, input$sec_variable]) - ))), + choices = val, + selected = ifelse(!is.null(val), val, NULL), + # choices = as.character(unique(as.factor(target[, input$sec_variable]))), + # selected = min(as.character(unique( + # as.factor(target[, input$sec_variable]) + # ))), multiple = TRUE ), h6('Other values will be considered as "0"', align = "center") @@ -8615,8 +8675,7 @@ CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT" ) updateNumericInput(session, "pcorrThreshold", - value = input$ValAlpha, - max = input$valAlpha + value = input$ValAlpha ) } else @@ -8746,18 +8805,40 @@ CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT" #It means that we got p-values if(!is.null(ppcor$statistic)){ - updateRadioButtons(session, "colorCorr", choiceNames = c("Color nodes according to correlation with a variable", paste0("Color edges according to partial correlation")), choiceValues = c("corr", "pcorr"), selected = "pcorr") + #updateRadioButtons(session, "colorCorr", choiceNames = c("Color nodes according to correlation with a variable", paste0("Color edges according to partial correlation")), choiceValues = c("corr", "pcorr")) p_val <- ppcor$p.value + + if(isTRUE(input$correctPval)){ + p_val_vector <- as.vector(p_val) + + # Apply the Bonferroni correction + p_val_adjusted_vector <- p.adjust(p_val_vector, method = "bonferroni") + # Reshape the vector back to the original matrix dimensions + p_val_adjusted_matrix <- matrix(p_val_adjusted_vector, nrow = nrow(p_val), ncol = ncol(p_val)) + + p_val <- p_val_adjusted_matrix + } adjacency <- matrix(0, nrow = nrow(p_val), ncol = ncol(p_val)) adjacency[p_val > (1 - as.numeric(min(input$AlphaVal, input$pcorrThreshold)))] <- 1 } else{ - updateRadioButtons(session, "colorCorr", choiceNames = c("Color nodes according to correlation with a variable", paste0("Color edges according to correlation")), choiceValues = c("corr", "pcorr"), selected = "pcorr") + #updateRadioButtons(session, "colorCorr", choiceNames = c("Color nodes according to correlation with a variable", paste0("Color edges according to correlation")), choiceValues = c("corr", "pcorr")) resCorrTest <- corr.test(countsMatrix, ci = FALSE, method = input$pcorrMethod) corr_matrix <- resCorrTest$r pval <- resCorrTest$p - pval_bool <- + if(isTRUE(input$correctPval)){ + p_val_vector <- as.vector(pval) + + # Apply the Bonferroni correction + p_val_adjusted_vector <- p.adjust(p_val_vector, method = "bonferroni") + # Reshape the vector back to the original matrix dimensions + p_val_adjusted_matrix <- matrix(p_val_adjusted_vector, nrow = nrow(pval), ncol = ncol(pval)) + + pval <- p_val_adjusted_matrix + } + + pval_bool <- t(apply(pval, 1, function(v) { sapply(v, function(x) { x < min(input$AlphaVal, input$pcorrThreshold) @@ -8777,11 +8858,20 @@ CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT" } } else{ - updateRadioButtons(session, "colorCorr", choiceNames = c("Color nodes according to correlation with a variable", paste0("Color edges according to correlation")), choiceValues = c("corr", "pcorr"), selected = "pcorr") + #updateRadioButtons(session, "colorCorr", choiceNames = c("Color nodes according to correlation with a variable", paste0("Color edges according to correlation")), choiceValues = c("corr", "pcorr")) resCorrTest <- corr.test(countsMatrix, ci = FALSE, method = input$pcorrMethod) corr_matrix <- resCorrTest$r pval <- resCorrTest$p - + if(isTRUE(input$correctPval)){ + p_val_vector <- as.vector(pval) + + # Apply the Bonferroni correction + p_val_adjusted_vector <- p.adjust(p_val_vector, method = "bonferroni") + # Reshape the vector back to the original matrix dimensions + p_val_adjusted_matrix <- matrix(p_val_adjusted_vector, nrow = nrow(pval), ncol = ncol(pval)) + + pval <- p_val_adjusted_matrix + } pval_bool <- t(apply(pval, 1, function(v) { sapply(v, function(x) { @@ -8896,6 +8986,32 @@ CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT" # input$modifwidthVisu, input$widthVisu, "auto" # ))) + output$ByNetworkCommunities <- renderUI({ + counts = dataMergeCounts()$counts + sumTot = rowSums(counts) + ord = order(sumTot, decreasing = TRUE) + Available_taxo = rownames(counts)[ord] + res = unique(Plot_network( + input, + ResDiffAnal(), + Available_taxo, + compute_pcor = compute_pcor(), + SelectTaxoPlotNetworkDebounce(), + qualiVariable, + colors = colorsVisuPlot(), + dataInput = dataInput() + )$taxo$Community) + # res_vector <- unlist(res) + # + # res_table <- table(res_vector) + # + # res <- names(res_table[res_table > 1]) + + selectizeInput("CommunityVisuPlot", + "Select the community", + choices = res) + }) + output$Community <- renderUI({ counts = dataMergeCounts()$counts sumTot = rowSums(counts) diff --git a/ui.R b/ui.R index 4f525c2..4c7f598 100644 --- a/ui.R +++ b/ui.R @@ -1202,6 +1202,9 @@ function(request) { conditionalPanel(condition="input.PlotVisuSelect!='Network' && input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Diversity' && input.PlotVisuSelect!='Scatterplot' && input.PlotVisuSelect!='Krona'", radioButtons("SelectSpecifTaxo","Select the features",c("Most abundant"="Most","All"="All", "Differential features" = "Diff", "Non differential features" = "NoDiff")) ), + conditionalPanel(condition = "input.SelectSpecifTaxo == 'ByNetwork'", + uiOutput("ByNetworkCommunities") + ), conditionalPanel(condition="input.PlotVisuSelect!='Network' && input.PlotVisuSelect!='Rarefaction' && input.PlotVisuSelect!='Diversity' && input.PlotVisuSelect!='Scatterplot' && (input.SelectSpecifTaxo=='Diff' || input.SelectSpecifTaxo=='NoDiff') && input.PlotVisuSelect!='Krona' ", selectizeInput("ContrastList_table_Visu","",choices = "", multiple = TRUE), radioButtons("UnionInterContrasts","Union or intersection ?",c("Union"="Union","Intersection"="Inter"),inline = TRUE) @@ -1223,7 +1226,7 @@ function(request) { ## _______network #### ### ### conditionalPanel(condition="input.PlotVisuSelect=='Network'", - radioButtons("colorCorr", label = "", choices = c("Color nodes according to correlation with a variable" = "corr", "Color edges according to partial correlation" = "pcorr"), selected = "pcorr"), + radioButtons("colorCorr", label = "", choices = c("Color nodes according to correlation with a variable" = "corr", "Color edges according to correlation" = "pcorr"), selected = "pcorr"), conditionalPanel(condition = "input.colorCorr == 'corr'", uiOutput("SelectSecVariable"), div(id = "ValueOfQualitativeVaribale", uiOutput("SelectValueQualiVar"))) ), @@ -1234,8 +1237,9 @@ function(request) { radioButtons("clusterWeightsParam", label = "Add weights as a clustering parameter", choices = c("Yes" = "Yes", "No" = "No"), selected = "No") ), conditionalPanel(condition="input.colorCorr=='pcorr' && input.PlotVisuSelect=='Network'", - radioButtons("pcorrMethod", label = "Select the partial correlation coefficient", choices = c("Pearson" = "pearson", "Kendall" = "kendall", "Spearman" = "spearman"), selected = "pearson"), - numericInput("pcorrThreshold", label = "Select the p-value threshold", min = 0.01, max= 0.05, value = 0.05, step = 0.001) + radioButtons("pcorrMethod", label = "Select the partial correlation coefficient", choices = c("Pearson" = "pearson", "Kendall" = "kendall", "Spearman" = "spearman"), selected = "spearman"), + numericInput("pcorrThreshold", label = "Select the p-value threshold", value = 0.05), + checkboxInput("correctPval", "use Bonferroni correction", value = FALSE) ), # conditionalPanel(condition="input.colorCorr=='pcorr' && input.PlotVisuSelect=='Network'", # sliderInput("permThreshold", label = "Select the number of permutation test iteration", min = 1, max= 100, value = 10, step = 1)