inst/Phenotypic_Correlation/app.R

In TBA: Collection of 'shiny' Apps for Tree Breeding Analysis

library(shiny)
library(shinybusy)
library(ggplot2)
library(reshape2)
library(readxl)
# Function to calculate phenotypic correlation matrix and generate plot

Phenotypic_Correlation <- function(data, low_color = "blue", high_color = "green", correlation_values_color = "black") {

  old_options <- options(scipen = 999)  # Save current options
  on.exit(options(old_options))         # Restore options when function exits

  # Convert the first two columns to factor type
  data[, 1:2] <- lapply(data[, 1:2], as.factor)

  # Convert the remaining columns to numeric
  data[, -c(1, 2)] <- lapply(data[, -c(1, 2)], as.numeric)

  # Extract trait names (excluding the first two columns)
  traits <- names(data)[-c(1, 2)][sapply(data[-c(1, 2)], is.numeric)]

  # Prepare a matrix to store correlations
  correlation_matrix1 <- matrix("", nrow = length(traits), ncol = length(traits),
                                dimnames = list(traits, traits))

  # Calculate correlations for each pair of traits
  for (i in 1:length(traits)) {
    for (j in 1:length(traits)) {
      trait1 <- traits[i]
      trait2 <- traits[j]

      if (i == j) {
        correlation_matrix1[i, j] <- "1"  # Set correlation to 1 if it's the same trait
      } else {
        # Perform linear regression for trait1
        formula1 <- as.formula(paste0("`", trait1, "` ~ `", names(data)[1], "` + `", names(data)[2], "`"))
        model1 <- lm(formula1, data = data)
        anova_result1 <- anova(model1)

        # Perform linear regression for trait2
        formula2 <- as.formula(paste0("`", trait2, "` ~ `", names(data)[1], "` + `", names(data)[2], "`"))
        model2 <- lm(formula2, data = data)
        anova_result2 <- anova(model2)

        # Calculate phenotypic variance for trait1 and trait2
        replication_levels <- nlevels(data[[1]])
        genotypic_variance1 <- round((anova_result1$`Mean Sq`[2] - anova_result1$`Mean Sq`[3]) / replication_levels,6)

        phenotypic_variance1<-genotypic_variance1+anova_result1$`Mean Sq`[3]
        genotypic_variance2 <- round((anova_result2$`Mean Sq`[2] - anova_result2$`Mean Sq`[3]) / replication_levels,6)
        phenotypic_variance2<-genotypic_variance2+anova_result2$`Mean Sq`[3]

        # Calculate covariance sums
        total_of_genotypes_trait1 <- tapply(data[[trait1]], data[[2]], sum)
        total_of_genotypes_trait2 <- tapply(data[[trait2]], data[[2]], sum)
        total_of_replication_trait1 <- tapply(data[[trait1]], data[[1]], sum)
        total_of_replication_trait2 <- tapply(data[[trait2]], data[[1]], sum)
        number_of_replication <- nlevels(data[[1]])
        number_of_genotype <- nlevels(data[[2]])
        Grand_total_trait1 <- sum(data[[trait1]])
        Grand_total_trait2 <- sum(data[[trait2]])
        CF <- ( Grand_total_trait1 * Grand_total_trait2 ) / ( number_of_replication * number_of_genotype )
        Total_SP <- round(sum(data[[trait1]] * data[[trait2]]) - CF,6)
        Genotypic_SP <- round((sum(total_of_genotypes_trait1 * total_of_genotypes_trait2) / number_of_replication) - CF,6)

        Replication_SP <- round((sum(total_of_replication_trait1 * total_of_replication_trait2) / number_of_genotype) - CF,6)
        Error_SP <- Total_SP - Genotypic_SP - Replication_SP
        DF_Replication <- number_of_replication - 1
        DF_Genotypes <- number_of_genotype - 1
        DF_Error <- DF_Replication * DF_Genotypes
        Replication_MP <- round(Replication_SP / DF_Replication,6)
        Genotypic_MP <- round(Genotypic_SP / DF_Genotypes,6)
        Error_MP <- round(Error_SP / DF_Error,6)
        Genotypic_Covariance <- round((Genotypic_MP - Error_MP) / number_of_replication,6)
        Phenotypic_Covariance <- Error_MP + Genotypic_Covariance

        # Calculate correlation
        correlation1 <- round(Phenotypic_Covariance / sqrt(phenotypic_variance1 * phenotypic_variance2), 4)

        # Calculate t-statistic and p-value
        n <- number_of_replication * number_of_genotype  # Number of observations
        df <- n - 2  # Degrees of freedom for correlation

        # Check for valid correlation
        if (!is.nan(correlation1)&& !is.na(correlation1)) {
          t_stat <- (correlation1) * (sqrt(df / (1 - (correlation1)^2)))  # Calculate t-statistic
          p_value <- 2 * pt(abs(t_stat), df = df, lower.tail = FALSE)  # Calculate two-tailed p-value
        } else {
          t_stat <- NA
          p_value <- NA
        }

        # Determine significance level symbol
        if (!is.nan(t_stat) && !is.na(t_stat)) {
          if (p_value < 0.05) {
            significance_symbol <- "*"   # Significant at 5%
          } else {
            significance_symbol <- "NS"  # Non-significant
          }
        } else {
          significance_symbol <- ""  # No significance symbol if t_stat is NA and NaN
        }

        # Format correlation value without scientific notation
        formatted_correlation1 <- format(correlation1, scientific = FALSE)

        # Fill in the correlation matrix with correlation value and significance symbol
        correlation_matrix1[i, j] <- paste(formatted_correlation1, significance_symbol, sep = "")
      }
    }
  }

  # Prepare a matrix to store formatted correlation values without significance symbols
  # Extract numeric correlations for plotting
  numeric_correlation_matrix1 <- matrix(as.numeric(gsub("[^0-9.-]", "", correlation_matrix1)),
                                        nrow = length(traits), ncol = length(traits),
                                        dimnames = list(traits, traits))

  # Melt the numeric correlation matrix for ggplot
  melted_corr_mat1 <- melt(numeric_correlation_matrix1)
  colnames(melted_corr_mat1) <- c("Var1", "Var2", "value")

  # Replace NaNs with NA for ggplot handling
  melted_corr_mat1$value[is.nan(melted_corr_mat1$value)] <- NA

  # Determine bar height based on the number of traits
  number_of_traits1 <- length(traits)
  bar_height1 <- ifelse(number_of_traits1 < 8, 9, 14)  # Adjust height for fewer than 8 traits

  # Plotting the correlation heatmap using ggplot
  plot1 <- ggplot(data = melted_corr_mat1, aes(x = Var1, y = Var2, fill = as.numeric(value))) +
    geom_tile(color = "white", linewidth = 0.5) +
    geom_text(aes(Var2, Var1, label = ifelse(
      value == 1, "1",
      ifelse(value == -1, "-1",ifelse(value == 0, "0", sub("\\.?0+$", "", format(value, scientific = FALSE))))), fontface = "bold"), color = correlation_values_color, size = 6.5) + # Reduced text size for correlation values
    scale_fill_gradient(low = low_color, high = high_color, limits = c(-1, 1), name = "Correlation Range") +
    guides(fill = guide_colorbar(barwidth = 3, barheight = bar_height1, title.position = "top", title.hjust = 0.5,title.vjust = 2,
                                 title.theme = element_text(size = 18,face = "bold"), label.theme = element_text(size = 16, face = "bold"))) +  # Increase size of color bar
    theme(panel.border = element_rect(color = "black", fill = NA, linewidth = 2)) +
    theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
    theme(axis.text.x = element_text(size = 15, angle = 45, hjust = 1, color = "black", face = "bold")) +
    theme(axis.text.y = element_text(size = 15, color = "black", face = "bold")) +
    theme(axis.ticks.length = unit(0.3, "cm")) +
    theme(axis.ticks = element_line(linewidth = 1)) +
    scale_x_discrete(expand = c(0, 0)) +  # Remove gap on x-axis
    scale_y_discrete(expand = c(0, 0))    # Remove gap on y-axis

  # Return the correlation matrix and plot
  return(list(correlation_matrix1 = correlation_matrix1, plot1 = plot1))
}



ui<-fluidPage(

  sidebarLayout(
    sidebarPanel(
      h3("Phenotypic Correlation", style = "color: blue; font-weight: bold;font-size: 30px;"),
      h3("Upload the data file", style = "font-weight: bold"),
      fileInput("file3", "Choose Excel File (.xlsx , .xls)", accept = c(".xlsx", ".xls")),
      h3("Results Visualization Options", style = "font-weight: bold"),
      selectInput("low_color_pheno", "Low Value Color:", c("blue", "green", "yellow", "orange", "red", "brown", "white", "pink", "purple", "lightgreen", "gray", "maroon", "lightblue")),
      selectInput("high_color_pheno", "High Value Color:", c("green", "red", "yellow", "orange", "blue", "brown", "white", "pink", "purple", "lightgreen", "gray", "maroon", "lightblue", "black")),
      selectInput("correlation_values_color", "Correlation Value Color:", c("black", "green", "red","yellow", "orange", "blue", "brown", "white", "pink", "purple", "lightgreen", "gray", "maroon", "lightblue")),
      actionButton("analyze_phenotypic", "Analyze", style = "color: #FFFFFF; background-color: #007BFF; border-color: #007BFF;margin-bottom: 10px;"),
      p("Instructions for data format:", style = "color: orange; font-weight: bold;font-size: 16px;"),
      p("Excel file name should not contain spaces (e.g., use 'Sample_Data.xlsx' instead of 'Sample Data.xlsx')", style = "color: red;font-weight: bold;font-size: 14px;"),
      p("First column: Replication", style = "color: red;font-weight: bold;font-size: 14px;"),
      p("Second column: Genotypes", style = "color: red;font-weight: bold;font-size: 14px;"),
      p("Subsequent columns: Trait values (e.g., DBH, PH, FW, SW, KW, OC)", style = "color: red;font-weight: bold;font-size: 14px;"),
      p("Trait names should be short (e.g., 'DBH' for Diameter at Breast Height)", style = "color: red;font-weight: bold;font-size: 14px;"),
      p("Note: The analysis is based on the Randomized Block Design (RBD)", style = "color: purple; font-weight: bold;font-size: 16px;"),
      downloadButton("download_pc_example", "Download Example Data",
                     style = "color: #FFFFFF; background-color: #28A745; border-color: #28A745; margin-bottom: 10px;"),
      p("The example dataset includes:170 genotypes, 3 replications for each genotype and 6 traits", style = "color: red;font-weight: bold;font-size: 14px;"),
      h3("Download Results", style = "font-weight: bold"),
      downloadButton("downloadMatrix_pheno", "Phenotypic Correlation Matrix (CSV)", style = "color: #00008B; font-weight: bold; width: 100%;white-space: normal;"),
      downloadButton("downloadPlotJPEG_pheno", "Phenotypic Correlation Heatmap Plot (JPEG)", style = "color: #00008B; font-weight: bold;width: 100%;white-space: normal;"),
      downloadButton("downloadPlotPNG_pheno", "Phenotypic Correlation Heatmap Plot (PNG)", style = "color: #00008B; font-weight: bold;width: 100%;white-space: normal;margin-bottom: 10px; "),
      p("For feedback, queries or suggestions,   email:  tbacafri@gmail.com",style = "color: darkgreen; font-weight: bold; font-size: 14px; width: 100%; white-space: normal;")

    ),
    mainPanel(

      uiOutput("titles_phenotypic"),


      # Wrap tableOutput in a div with CSS for scrolling
      div(style = "overflow-y: auto;auto;overflow-x: auto; height: 400px;",   # Set height and enable both scrolls
          tableOutput("correlation_output_pheno")
      ),

      uiOutput("filtering_ui_pheno"),  # Filtering UI

      uiOutput("annotations_output_pheno"),
      uiOutput("plot_loading_message1"),

      div(
        style = "overflow-x: auto; overflow-y: auto; width: 100%; ",
        uiOutput("dynamic_plot_ui_pheno")  # Dynamically generated plot UI
      )

    )
  )
)
server <- function(input, output, session) {
  ######## Phenotypic correlation analysis logic    ############

  output$download_pc_example <- downloadHandler(
    filename = function() {
      "Phenotypic_Correlation_Data.xlsx"  # File name when user downloads
    },
    content = function(file) {
      # Locate the file in the package's inst/Phenotypic_Correlation folder
      example_path <- system.file("Phenotypic_Correlation", "example_PC_data.xlsx", package = "TBA")

      # Copy that file to the temp download location
      file.copy(example_path, file)
    }
  )
  correlation_analysis_pheno <- reactiveVal()

  analyzePhenotypic <- function(file, low_color, high_color, correlation_values_color) {
    req(file)
    data <- readxl::read_excel(file$datapath)
    res <- Phenotypic_Correlation(data, low_color, high_color, correlation_values_color)
    correlation_analysis_pheno(res)
    return(res)
  }

  # Reset previous outputs when a new file is uploaded
  observeEvent(input$file3, {
    correlation_analysis_pheno(NULL)  # Clear the analysis results
    output$correlation_output_pheno <- renderUI(NULL)  # Clear correlation table output
    output$dynamic_plot_ui_pheno <- renderUI(NULL)     # Clear dynamic plot UI
    output$annotations_output_pheno <- renderUI(NULL)   # Clear annotations
    output$titles_phenotypic <- renderUI(NULL)     # Clear titles
  })

  observeEvent(input$analyze_phenotypic, {

    show_modal_spinner(
      spin = "circle",
      color = "#007BFF",
      text = "Analyzing, please wait..."   # (Optional text under spinner)
    )

    res <- analyzePhenotypic(input$file3, input$low_color_pheno, input$high_color_pheno, input$correlation_values_color)

    remove_modal_spinner()

    # Render the correlation table
    output$correlation_output_pheno <- renderTable({
      if (!is.null(res)) {
        res$correlation_matrix1  # Make sure this matches the correct matrix for pheno
      }
    }, rownames = TRUE)


    # Dynamically calculate the plot width and height based on the number of traits

    output$dynamic_plot_ui_pheno <- renderUI({
      req(correlation_analysis_pheno())  # Ensure correlation analysis result is available

      # Get the number of traits (columns) in the correlation matrix
      num_traits <- length(colnames(correlation_analysis_pheno()$correlation_matrix1))

      # Dynamically calculate plot dimensions based on the number of traits
      plot_width1 <- paste0(150 + num_traits * 100, "px")  # Adjust width for more traits
      plot_height1 <- paste0(150 + num_traits * 50, "px")  # Adjust height for more traits

      # Render plot with dynamic dimensions
      plotOutput("plot_output_pheno", height = plot_height1, width = plot_width1)
    })
    # Show message before rendering plot
    output$plot_loading_message1 <- renderUI({
      tags$p("Generating heatmap plot, please wait...", style = "color: orange; font-weight: bold; font-size: 16px;")
    })


    # Render the plot itself

    output$plot_output_pheno <- renderPlot({
      if (!is.null(correlation_analysis_pheno())) {
        # Clear the message once the plot is ready
        output$plot_loading_message1 <- renderUI(NULL)
        print(correlation_analysis_pheno()$plot1)
      }
    })



    output$annotations_output_pheno <- renderUI({
      if (!is.null(correlation_analysis_pheno())) {
        tagList(
          p(tags$b(tags$span("NS: Non-significant", style = "font-size: 18px;"))),
          p(tags$b(tags$span("*: Significant at 5% (p<0.05)", style = "font-size: 18px;"))),
          p(tags$b(tags$span("Note: The significance of phenotypic correlation tested using t test (two-tail)."))),
          h3("Phenotypic Correlation Heatmap Plot", style = "color: purple; font-weight: bold;")
        )
      }
    })
    output$titles_phenotypic <- renderUI({
      tagList(
        h3("Phenotypic Correlation Matrix", style = "color: purple; font-weight: bold;")
      )
    })

    #Select a Variable to see its correlation with other variables

    output$filtering_ui_pheno <- renderUI({
      req(correlation_analysis_pheno())

      tagList(
        selectInput(
          inputId = "selected_variable",
          label = "Select a Trait to see its correlation with other Traits:",
          choices = c("", rownames(correlation_analysis_pheno()$correlation_matrix1)),  # Add an empty option
          selected = ""  # Set the default to an empty value
        ),
        tableOutput("filtered_correlation_table_pheno")
      )
    })

    output$filtered_correlation_table_pheno <- renderTable({
      req(input$selected_variable)  # Ensure a variable is selected
      req(input$selected_variable != "")  # Ensure the selected variable is not the placeholder
      req(correlation_analysis_pheno())   # Ensure the correlation analysis results are available
      selected_variable <- input$selected_variable
      corr_matrix <- correlation_analysis_pheno()$correlation_matrix1

      # Construct the data frame manually
      variables <- colnames(corr_matrix)
      correlations <- corr_matrix[selected_variable, ]


      # Remove the row corresponding to the selected variable itself
      filter_result <- data.frame(
        Traits = variables,
        Correlation = correlations,
        stringsAsFactors = FALSE
      )
      filter_result <- filter_result[filter_result$Traits!= selected_variable, ]

      return(filter_result)
    }, rownames = FALSE)

  })

  output$downloadMatrix_pheno <- downloadHandler(
    filename = function() {
      paste("Phenotypic_Correlation_Matrix", Sys.Date(), ".csv", sep = "")
    },
    content = function(file) {
      write.csv(correlation_analysis_pheno()$correlation_matrix1, file, row.names = TRUE)

    }
  )


  output$downloadPlotPNG_pheno <- downloadHandler(
    filename = function() {
      paste("Phenotypic_Correlation_Heatmap_plot", Sys.Date(), ".png", sep = "")
    },
    content = function(file) {
      req(correlation_analysis_pheno())
      # Get the number of traits (columns) in the correlation matrix
      num_traits <- length(colnames(correlation_analysis_pheno()$correlation_matrix1))

      # Dynamically calculate the plot dimensions in inches for ggsave
      plot_width_inch1 <- 1.5 + num_traits * 1.5 # Adjust width (inches) based on traits
      plot_height_inch1 <- 1.5 + num_traits * 0.5  # Adjust height (inches) based on traits

      # Save the plot with dynamic dimensions
      ggsave(file, correlation_analysis_pheno()$plot1, device = "png",
             width = plot_width_inch1, height = plot_height_inch1, dpi = 600)

    }
  )

  output$downloadPlotJPEG_pheno <- downloadHandler(
    filename = function() {
      paste("Phenotypic_Correlation_Heatmap_plot", Sys.Date(), ".jpeg", sep = "")
    },
    content = function(file) {
      req(correlation_analysis_pheno())
      # Get the number of traits (columns) in the correlation matrix
      num_traits <- length(colnames(correlation_analysis_pheno()$correlation_matrix1))

      # Dynamically calculate the plot dimensions in inches for ggsave
      plot_width_inch1 <- 1.5 + num_traits * 1.5  # Adjust width (inches) based on traits
      plot_height_inch1 <- 1.5 + num_traits * 0.5  # Adjust height (inches) based on traits

      # Save the plot with dynamic dimensions
      ggsave(file, correlation_analysis_pheno()$plot1, device = "jpeg",
             width = plot_width_inch1, height = plot_height_inch1, dpi = 600)

    }
  )


}

shinyApp(ui, server)