inst/Genotypic_Path/app.R

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

library(shiny)
library(readxl)
library(shinybusy)

Genotypic_Path<- function(data) {
  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_matrix <- matrix(NA, nrow = length(traits), ncol = length(traits))
  formatted_correlation_matrix <- matrix(NA, nrow = length(traits), ncol = length(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_matrix[i, j] <- 1 # Set correlation to 1 if it's the same trait
        formatted_correlation_matrix[i, j] <- 1  # Set formatted correlation value to 1
      } 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,4)


        genotypic_variance2 <- round((anova_result2$`Mean Sq`[2] - anova_result2$`Mean Sq`[3]) / replication_levels,4)


        # 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,4)
        Genotypic_SP <- round((sum(total_of_genotypes_trait1 * total_of_genotypes_trait2) / number_of_replication) - CF,4)
        Replication_SP <- round((sum(total_of_replication_trait1 * total_of_replication_trait2) / number_of_genotype) - CF,4)
        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,4)
        Genotypic_MP <- round(Genotypic_SP / DF_Genotypes,4)
        Error_MP <- round(Error_SP / DF_Error,4)

        Genotypic_Covariance <- round((Genotypic_MP - Error_MP) / number_of_replication,4)



        # Calculate correlation
        correlation <- round(Genotypic_Covariance / sqrt(genotypic_variance1 * genotypic_variance2), 4)

        # Perform significance test
        n <- nlevels(data[[2]]) # Number of observations means genotypes
        df <- n - 2  # Degrees of freedom for Pearson correlation
        if (!is.nan(correlation)&& !is.na(correlation)) {
          t_stat <- (correlation) * (sqrt(df / (1 - (correlation)^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 1%
          }else {
            significance_symbol <- "NS"  # Non-significant
          }
        } else {
          significance_symbol <- ""  # No significance symbol if t_stat is NA
        }

        # Store correlation value in the matrices

        formatted_correlation_matrix[i, j] <- correlation
        correlation_matrix[i, j] <- paste0(format(correlation, scientific = FALSE), significance_symbol)
      }
    }
  }
  genotypic_correlation_matrix <- noquote(correlation_matrix)
  correlation_only <- noquote(formatted_correlation_matrix)

  # Path Analysis
  dependent_variable <- correlation_only[1:(length(traits) - 1), length(traits)]
  dependent_variable_matrix <- matrix(dependent_variable, ncol = 1)
  independent_variable <- correlation_only[1:(length(traits) - 1), 1:(length(traits) - 1)]
  direct_effect <- solve(independent_variable,dependent_variable_matrix)
  Direct_and_indirect_effect <- matrix(nrow = (length(traits) - 1), ncol = (length(traits) - 1))

  for (i in 1:(length(traits) - 1)) {
    for (j in 1:(length(traits) - 1)) {
      Direct_and_indirect_effect[i, j] <- round(direct_effect[j] * independent_variable[i, j], 4)
    }
  }

  Path_effects <- cbind(Direct_and_indirect_effect, genotypic_correlation_matrix[1:(length(traits) - 1), length(traits)])
  rownames(Path_effects) <- traits[1:(length(traits) - 1)]
  colnames(Path_effects) <- traits[1:(length(traits))]

  residual <- 1 - t(direct_effect) %*% dependent_variable_matrix
  Residual_effect <- round(sqrt(residual), 4)

  rownames(Direct_and_indirect_effect) <- traits[1:(length(traits) - 1)]
  colnames(Direct_and_indirect_effect) <- traits[1:(length(traits) - 1)]
  rownames(Residual_effect)<-"Residual Effect"

  # Convert Path_effects to data frame
  Path_effects <- as.data.frame(Path_effects)

  # Return the data frame with row names and residual effect
  return(list(Path_effects = Path_effects, Residual_effect = Residual_effect))
}


ui<-fluidPage(

  sidebarLayout(
    sidebarPanel(
      h3("Genotypic Path Analysis", style = "color: blue; font-weight: bold;font-size: 30px;"),
      h3("Upload the data file", style = "font-weight: bold"),
      fileInput("file_path_analysis_geno", "Choose Excel File (.xlsx , .xls)", accept = c(".xlsx", ".xls")),
      actionButton("analyze_path_analysis_geno", "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(" The last column must be the dependent trait for path analysis. For example, if OC (Oil Content) is the dependent trait, it should appear in the last column.", 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_gp_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 (5 independent traits and 1 dependent trait: OC)", style = "color: red;font-weight: bold;font-size: 14px;"),
      h3("Download Results", style = "font-weight: bold"),
      downloadButton("downloadPathEffectsCSVGenotype", "Genotypic Path Analysis Results (CSV)", style = "color: #00008B; font-weight: bold; width: 100%;white-space: normal;margin-bottom: 15px;"),
      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("pathEffectsTitleGenotype"),
      # Wrap tableOutput in a div with CSS for vertical scrolling
      div(style = "overflow-y: auto;overflow-x: auto; height: 400px;",  # Adjust height as needed
          tableOutput("pathEffectsTableGenotype")
      ),
      uiOutput("annotations_path_effects_genotype"),
      uiOutput("residualEffectGenotype")

    )
  )
)

server<-function(input, output, session) {
  ###### Genotypic Path Analysis logic  #######
  output$download_gp_example <- downloadHandler(
    filename = function() {
      "Genotypic_Path_Data.xlsx"  # File name when user downloads
    },
    content = function(file) {
      # Locate the file in the package's inst/Genotypic_Correlation folder
      example_path <- system.file("Genotypic_Path", "example_GP_data.xlsx", package = "TBA")

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

  analyzePathAnalysisGeno <- function(file) {
    req(file)
    data <- readxl::read_excel(file$datapath)
    res <- Genotypic_Path(data)
    path_effects_geno(res$Path_effects)
    residual_effect_geno(res$Residual_effect)
    return(res)
  }

  # Reset previous outputs when a new file is uploaded
  observeEvent(input$file_path_analysis_geno, {
    path_effects_geno(NULL)  # Clear the analysis results
    residual_effect_geno(NULL) # Clear the analysis results
    output$pathEffectsTableGenotype <- renderUI(NULL)  # Clear table output
    output$residualEffectGenotype <- renderUI(NULL)     # Clear output
    output$pathEffectsTitleGenotype<-renderUI(NULL)  # Clear Title
    output$annotations_path_effects_genotype <- renderUI(NULL)   # Clear annotations
  })

  observeEvent(input$analyze_path_analysis_geno, {

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

    res <- analyzePathAnalysisGeno(input$file_path_analysis_geno)
    remove_modal_spinner()


    output$pathEffectsTableGenotype <- renderTable({
      if (!is.null(path_effects_geno())) {
        path_effects_geno()
      }
    }, rownames = TRUE)

    output$residualEffectGenotype <- renderUI({
      if (!is.null(residual_effect_geno())) {
        tagList(
          h4("Residual Effect", style = "color: purple; font-weight: bold;"),
          p(residual_effect_geno())
        )
      }
    })

    output$annotations_path_effects_genotype <- renderUI({
      HTML("<br><b>Note:</b> Diagonal values show direct effects")

    })

    output$pathEffectsTitleGenotype <- renderUI({
      tagList(
        h3("Genotypic Path Analysis Results", style = "color: purple; font-weight: bold;")
      )
    })
  })

  output$downloadPathEffectsCSVGenotype <- downloadHandler(
    filename = function() {
      paste("genotypic_path_effects", Sys.Date(), ".csv", sep = "")
    },
    content = function(file) {
      write.csv(path_effects_geno(), file, row.names = TRUE)

    }
  )

}

shinyApp(ui, server)