## ----global-setup-appendix-b, include = FALSE---------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 10,
    fig.height = 6.5,
    message = FALSE,
    warning = FALSE,
    dev = "png",
    dpi = 300
)
suppressPackageStartupMessages(library(kableExtra))

## ----load-packages, message=FALSE---------------------------------------------
suppressPackageStartupMessages({
    library(ggplot2)
    library(SummarizedExperiment)
    library(dplyr)
    library(gridExtra)
})

set.seed(42)

if (!requireNamespace("devtools", quietly = TRUE)) {
    stop("Please install devtools to run this vignette against the local package source.")
}
if ("package:TSENAT" %in% search()) {
    detach("package:TSENAT", unload = TRUE, character.only = TRUE)
}

pkg_root <- normalizePath("..", mustWork = FALSE)
devtools::load_all(pkg_root)

# Load preprocessed dataset
data(readcounts)
readcounts <- as.matrix(readcounts)
mode(readcounts) <- "numeric"

# Load metadata
metadata_df <- read.table(
    system.file("extdata", "metadata.tsv", package = "TSENAT"),
    header = TRUE, sep = "\t"
)

gff3_dataset <- system.file("extdata", "annotation.gff3.gz", package = "TSENAT")

# Configure analysis parameters first (best practice: fail-fast validation)
config <- TSENAT_config(
  sample_col = "sample",
  condition_col = "condition",
  subject_col = "paired_samples",
  q = seq(0, 2, by = 0.05),
  nthreads = 3,
  paired = TRUE,
  control = "normal"
)

# Build analysis object: creates SummarizedExperiment and initializes TSENATAnalysis
analysis <- build_analysis(
    readcounts = readcounts,
    tx2gene = gff3_dataset,
    metadata = metadata_df,
    tpm = tpm,
    effective_length = effective_length,
    config = config
)

# Apply filtering for quality control
analysis <- filter_analysis(
    analysis,
    stringency = "medium"
)

## ----compute-tsallis-entropy, message=FALSE, results='hide'-------------------
# Compute diversity using S4 wrapper with bootstrap confidence intervals
analysis <- calculate_diversity(
  analysis, 
  norm = TRUE,
  pseudocount = "auto"
)

## ----validate-rank-transform-assumptions, message=FALSE, results='hide'-------
# Validate statistical assumptions (including GAMM diagnostics)
analysis <- calculate_assumptions(
    analysis,
    checks = "all"  # Include core checks + GAMM diagnostics
)

# Get assumptions
assumptions_text <- results(analysis, type = "assumptions")
print(assumptions_text)

## ----display-assumption-results, message=FALSE, echo=FALSE--------------------
# Extract assumption check results as structured list
check_results <- results(analysis, type = "assumptions", format = "list")

# Display assumptions table as formatted kable
if (!is.null(check_results) && !is.null(check_results$assumptions_table)) {
  assumptions_df <- check_results$assumptions_table
  
  knitr::kable(
    assumptions_df,
    caption = "**Supplementary Table 14 | Assumption Checks for Aligned Rank Transform (ART).** Evaluates homogeneity of variance and normality for rank-based interaction test validity."
  )
}

## ----rank-transform-test-analysis, message=FALSE, results='hide'--------------
# Run Aligned Rank Transform (ART, default method)
# This preserves the GAMM results in 'analysis' for later concordance comparison

# Run ART for q × condition interaction
analysis_rank <- calculate_rank_transform(
    analysis,
    multicorr = "hochberg"
)

# Extract ALL results (no filtering) for summary statistics
rank_transform_results_all <- results(analysis_rank, type = "rank_test", rankBy = "pvalue")

# Extract top significant genes for display
rank_transform_results <- results(analysis_rank, type = "rank_test", rankBy = "pvalue", n = 20, filterFDR = 0.05)
print(head(rank_transform_results, n = 10))

## ----rank-transform-results-summary, message=FALSE, echo=FALSE----------------
# Create summary statistics table with additional metrics using ALL tested genes
# Handle NA and NaN values properly
mean_effect <- mean(rank_transform_results_all$effect_size_eta2, na.rm = TRUE)
median_effect <- median(rank_transform_results_all$effect_size_eta2, na.rm = TRUE)

rank_summary_stats <- data.frame(
  Metric = c(
    "Genes tested", 
    "Significant (p < 0.05)", 
    "Significant (adj_p < 0.05, FWER-controlled)", 
    "NAs",
    "Mean effect size ($\\eta^2$)",
    "Median effect size ($\\eta^2$)",
    "Strong effect genes ($\\eta^2$ > 10%)"
  ),
  Value = c(
    nrow(rank_transform_results_all),
    sum(rank_transform_results_all$p_value < 0.05, na.rm = TRUE),
    sum(rank_transform_results_all$adj_p_value < 0.05, na.rm = TRUE),
    sum(is.na(rank_transform_results_all$p_value)),
    if (is.na(mean_effect) || !is.finite(mean_effect)) "Not available" else sprintf("%.1f%%", mean_effect * 100),
    if (is.na(median_effect) || !is.finite(median_effect)) "Not available" else sprintf("%.1f%%", median_effect * 100),
    sum(rank_transform_results_all$p_value < 0.05 & rank_transform_results_all$effect_size_eta2 > 0.10, na.rm = TRUE)
  )
)

knitr::kable(rank_summary_stats, caption = "**Supplementary Table 15 | Aligned Rank Transform (ART) Results Summary.** Non-parametric test of *q* × condition interactions via ART. Compare with GAMM results (main vignette Table 3) to assess parametric vs. non-parametric method concordance.", label = "rank-transform-test-summary", align = c("l", "r"))

# Show top genes results
rank_sig <- rank_transform_results_all[!is.na(rank_transform_results_all$p_value), ]
if (nrow(rank_sig) > 0) {
  rank_sig <- rank_sig[order(rank_sig$adj_p_value), ]
  rank_sig <- rank_sig[1:min(10, nrow(rank_sig)), ]
  
  # Select columns: gene, p-values, f-statistic, effect size, interaction class
  cols_to_show <- intersect(c("gene", "p_value", "adj_p_value", "f_statistic", "effect_size_eta2", 
                              "interaction_class"), 
                            colnames(rank_sig))
  if (length(cols_to_show) > 0) {
    rank_display <- rank_sig[, cols_to_show, drop = FALSE]
    
    # Format p-values with scientific notation for very small values
    if ("p_value" %in% colnames(rank_display)) {
      rank_display$p_value <- sapply(rank_display$p_value, function(x) {
        if (is.na(x)) "NA" else if (x >= 0.9999) "1" else if (x < 0.0001) format(signif(x, 2), scientific=TRUE) else sprintf("%.6f", x)
      })
    }
    
    if ("adj_p_value" %in% colnames(rank_display)) {
      rank_display$adj_p_value <- sapply(rank_display$adj_p_value, function(x) {
        if (is.na(x)) "NA" else if (x >= 0.9999) "1" else if (x < 0.0001) format(signif(x, 2), scientific=TRUE) else sprintf("%.6f", x)
      })
    }
    
    # Replace test method name for clarity (removed - column no longer included)
    # if ("test_method" %in% colnames(rank_display)) {
    #   rank_display$test_method <- sapply(rank_display$test_method, function(x) {
    #     if (is.na(x)) "NA" else if (x == "rt_paired") "Conover-Iman Rank Transform (paired)" else x
    #   })
    # }
    
    # Format other numeric columns
    if ("f_statistic" %in% colnames(rank_display)) {
      rank_display$f_statistic <- sapply(rank_display$f_statistic, function(x) {
        if (is.na(x)) "NA" else sprintf("%.4f", x)
      })
    }
    
    if ("effect_size_eta2" %in% colnames(rank_display)) {
      rank_display$effect_size_eta2 <- sapply(rank_display$effect_size_eta2, function(x) {
        if (is.na(x)) "NA" else sprintf("%.4f", x)
      })
    }
    
    knitr::kable(rank_display, 
                 caption = "**Supplementary Table 16 | Top genes identified by Aligned Rank Transform (ART).** Ranked by adjusted p-value; comparison with GAMM (main vignette Table 3) reveals method-specific sensitivities.",
                 col.names = c("Gene", "P-value", "Adj. P-value", "F-Statistic", 
                               "Effect Size (η²)", "Interaction Class")[1:ncol(rank_display)],
                 row.names = FALSE)
  }
}

## ----extended-fig-1-rank-transform-q-curves, fig.width=12, fig.height=8, message=FALSE, out.width="98%", fig.pos="H", fig.cap="**Supplementary Figure 1 | Aligned Rank Transform (ART) results for scale-dependent isoform switching.** Q-curves for top genes identified by ART; compare visual patterns with GAMM results (main vignette Figure 2) to assess method agreement."----
# Plot top genes from Aligned Rank Transform (ART) (using all genes, not just significant ones)
# This ensures the plot displays n_top genes regardless of significance threshold
top_genes_plot <- plot_diversity_spectrum(analysis, sait_res = rank_transform_results_all, n_top = 4)
print(top_genes_plot)

## ----generalized-additive-model-analysis, message=FALSE, results='hide'-------
# Load precomputed GAMM analysis object from RDS file
analysis_sait <- readRDS(
    system.file("extdata", "analysis_sait.rds", package = "TSENAT")
)

# Extract GAMM results for inspection using accessor function
sait_results <- results(analysis_sait, type = "sait", rankBy = "pvalue")
print(sait_results)


## ----gam-results-summary-table, message=FALSE, echo=FALSE---------------------
# Create summary statistics table with additional metrics
# Handle NA and NaN values properly
mean_effect <- mean(sait_results$effect_size, na.rm = TRUE)
median_effect <- median(sait_results$effect_size, na.rm = TRUE)

summary_stats <- data.frame(
  Metric = c(
    "Genes tested", 
    "Significant (p < 0.05)", 
    "Concordant (p < 0.05 AND adj_p < 0.05)",
    "NAs",
    "Mean effect size",
    "Median effect size",
    "Strong effect genes (effect_size > 30%)",
    "Model convergence rate"
  ),
  Value = c(
    nrow(sait_results),
    sum(sait_results$p_interaction < 0.05, na.rm = TRUE),
    sum(sait_results$p_interaction < 0.05 & sait_results$adj_p_interaction < 0.05, na.rm = TRUE),
    sum(is.na(sait_results$p_interaction)),
    if (is.na(mean_effect) || !is.finite(mean_effect)) "Not available" else sprintf("%.1f%%", mean_effect * 100),
    if (is.na(median_effect) || !is.finite(median_effect)) "Not available" else sprintf("%.1f%%", median_effect * 100),
    sum(sait_results$p_interaction < 0.05 & sait_results$effect_size > 0.30, na.rm = TRUE),
    sprintf("%.1f%%", mean(sait_results$model_converged, na.rm = TRUE) * 100)
  )
)

knitr::kable(summary_stats, caption = "**Supplementary Table 17 | GAMM Results Summary.** Overview of generalized additive mixed model statistics for q × condition interaction tests across full q-spectrum with paired sample structure.", label = "gam-summary-stats", align = c("l", "r"))

# Show top genes with additional information
if ("p_interaction" %in% colnames(sait_results)) {
  gam_sig <- sait_results[!is.na(sait_results$p_interaction), ]
  
  if (nrow(gam_sig) > 0) {
    gam_sig <- gam_sig[order(gam_sig$p_interaction), ]
    gam_sig <- gam_sig[1:min(10, nrow(gam_sig)), ]
    
    # Select columns for display: gene, p-values, effect size, test statistic
    cols_to_show <- intersect(c("gene", "p_interaction", "adj_p_interaction", "effect_size", 
                                "test_statistic"), 
                              colnames(gam_sig))
    if (length(cols_to_show) > 0) {
      gam_display <- gam_sig[, cols_to_show, drop = FALSE]
      
      # Format numeric columns (but NOT p-value columns - preserve full precision for later formatting)
      numeric_cols <- sapply(gam_display, is.numeric)
      # Only round non-p-value columns to 3 decimal places
      pval_cols <- c("p_interaction", "adj_p_interaction")
      cols_to_round <- setdiff(names(numeric_cols[numeric_cols]), pval_cols)
      if (length(cols_to_round) > 0) {
        gam_display[, cols_to_round] <- round(gam_display[, cols_to_round], 3)
      }
      
      # Format p-value columns with scientific notation for very small values
      if ("p_interaction" %in% colnames(gam_display)) {
        gam_display$p_interaction <- sapply(gam_display$p_interaction, function(x) {
          if (is.na(x)) "NA" else if (x < 0.001) sprintf("%.2e", x) else sprintf("%.6f", x)
        })
      }
      
      if ("adj_p_interaction" %in% colnames(gam_display)) {
        gam_display$adj_p_interaction <- sapply(gam_display$adj_p_interaction, function(x) {
          if (is.na(x)) "NA" else if (x < 0.001) sprintf("%.2e", x) else sprintf("%.6f", x)
        })
      }
      
      # Format effect size as percentage
      if ("effect_size" %in% colnames(gam_display)) {
        gam_display$effect_size <- sapply(gam_display$effect_size, function(x) {
          if (is.na(x)) "NA" else sprintf("%.1f%%", x * 100)
        })
      }
      
      # Format test statistic
      if ("test_statistic" %in% colnames(gam_display)) {
        gam_display$test_statistic <- sapply(gam_display$test_statistic, function(x) {
          if (is.na(x)) "NA" else sprintf("%.2f", x)
        })
      }
      
      # Rename columns for better readability — dynamic to match actual columns present
      col_name_map <- c(
        "gene" = "Gene",
        "p_interaction" = "p (interaction)",
        "adj_p_interaction" = "Adjusted p",
        "effect_size" = "Effect size",
        "test_statistic" = "Test statistic"
      )
      present_cols <- intersect(names(col_name_map), colnames(gam_display))
      gam_display <- gam_display[, present_cols, drop = FALSE]
      colnames(gam_display) <- unname(col_name_map[present_cols])
      n_cols <- ncol(gam_display)
      
      knitr::kable(gam_display, 
                   caption = "**Supplementary Table 18 | Top 10 Genes by GAMM p-value (with Effect Size and Test Statistics).** Ranked by adjusted p-value; compare with ART results (Supplementary Table 16) for method agreement assessment.",
                   label = "gam-top-genes",
                   align = c("l", rep("r", max(0, n_cols - 1))),
                   row.names = FALSE)
    }
  }
} else {
  cat("WARNING: p_interaction column not found in results\n")
  cat("Available columns:", paste(colnames(sait_results), collapse=", "), "\n")
}

## ----compare-rank-transform-gam-methods, message=FALSE------------------------
# Compute concordance analysis using new two-object API
rank_results <- results(analysis_rank, type = "rank_test")
if (is.null(rank_results) || length(rank_results) == 0) {
    stop("Cannot run calculate_concordance(): rank_test results are empty. Ensure calculate_rank_transform() completed successfully.")
}

analysis_rank <- calculate_concordance(
    analysis_sait = analysis_sait,
    analysis_rank = analysis_rank,
    verbose = TRUE
)

# Extract concordance results with list format for component display
concordance_results <- results(analysis_rank, type = "concordance", format = "list")

## ----display-concordance-metrics, message=FALSE, echo=FALSE, results='asis'----
# Display concordance summary - results() with format="list" returns structured components
if (!is.null(concordance_results) && is.list(concordance_results)) {
  # Display summary table
  if (!is.null(concordance_results$summary_table) && nrow(concordance_results$summary_table) > 0) {
    print(knitr::kable(concordance_results$summary_table, 
                 caption = "**Supplementary Table 19 | Global Concordance Metrics: GAMM vs Aligned Rank Transform (ART).** Quantifies agreement between parametric (GAMM) and non-parametric (ART) approaches for detecting q × condition interactions.",
                 align = c("l", "r"),
                 row.names = FALSE))
  }
  
  # Display agreement distribution
  if (!is.null(concordance_results$agreement_dist) && nrow(concordance_results$agreement_dist) > 0) {
    cat("\n")
    print(knitr::kable(concordance_results$agreement_dist, 
                 caption = "**Supplementary Table 20 | Agreement Distribution: Concordance by Gene Category.** Distribution of genes across categories: both methods significant, GAMM only, ART only, and neither method.",
                 align = c("l", "r", "r"),
                 row.names = FALSE))
  }
  
  # Display high-confidence genes (concordant in both methods)
  if (!is.null(concordance_results$high_conf_table) && nrow(concordance_results$high_conf_table) > 0) {
    cat("\n")
    n_hc <- nrow(concordance_results$high_conf_table)
    print(knitr::kable(concordance_results$high_conf_table, 
                 caption = paste0("**Supplementary Table 21 | Robust Entropic Order Index Interactions: High-Confidence Genes Detected by Both Methods (n=", n_hc, ", ranked by statistical significance).** Cross-method validation: genes significant in both GAMM and ART provide strongest evidence for q × condition effects."),
                 align = c("l", "r", "r", "r", "r"),
                 row.names = FALSE))
  }
} else {
  cat("Concordance results not available\n")
}

## ----visualize-method-concordance, fig.width=10, fig.height=8, out.width="98%", fig.pos="H", fig.cap="**Supplementary Figure 2 | Method concordance visualization comparing Aligned Rank Transform (ART) and GAMM approaches for detecting q × condition interactions.** Scatter plot showing overlap in genes identified as significant by each method; reveals parametric-non-parametric agreement patterns."----
# Create comparison plot using S4 wrapper
# Concordance analysis results are stored in analysis metadata after calculate_concordance()
p <- plot_concordance(analysis_rank, verbose = TRUE)
grid::grid.draw(p)

## ----session-info-appendix-b--------------------------------------------------
sessionInfo()

