## ----global-setup, include = FALSE--------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 10,
    fig.height = 6.5,
    echo = TRUE,
    eval = TRUE,
    results = 'hide',
    message = TRUE,
    warning = TRUE
)
suppressPackageStartupMessages(library(kableExtra))

## ----install-bioconda, eval = FALSE-------------------------------------------
# # Install from Bioconda
# conda install -c bioconda r-tsenat

## ----install-github, eval = FALSE---------------------------------------------
# remotes::install_github("gallardoalba/TSENAT")

## ----quick-start-example, eval = FALSE----------------------------------------
# library(TSENAT)
# library(SummarizedExperiment)
# 
# # Load example data
# data("readcounts", package = "TSENAT")
# 
# # Load readcounts object, which contains Salmon-quantified
# # fragment counts with EM-based ambiguity resolution applied
# readcounts <- as.matrix(readcounts)
# 
# # Load sample metadata and annotation
# metadata_df <- read.table(
#   system.file("extdata", "metadata.tsv", package = "TSENAT"),
#   header = TRUE, sep = "\t"
# )
# gff3_file <- system.file("extdata", "annotation.gff3.gz", package = "TSENAT")
# 
# # Create analysis object with transcript counts, annotation, and metadata
# config <- TSENAT_config(
#   sample_col = "sample",
#   condition_col = "condition",
#   paired = TRUE,
#   subject_col = "paired_samples",
#   control = "normal",
#   q = seq(0, 2, by = 0.05)
# )
# 
# ## Build TSENATAnalysis object
# analysis <- build_analysis(
#   readcounts = readcounts,
#   tx2gene = gff3_file,
#   metadata = metadata_df,
#   config = config,
#   tpm = tpm,
#   effective_length = effective_length)
# 
# # Filter low-abundance transcripts
# analysis <- filter_analysis(analysis)
# 
# # Compute Tsallis entropy using S4 wrapper (using single q value for quick start)
# analysis <- calculate_diversity(analysis)
# 
# # Plot overall q-curve
# p_qcurve <- plot_diversity_spectrum(
#     analysis,
#     dev_width = 12,
#     dev_height = 8)
# print(p_qcurve)

## ----workflow-setup-----------------------------------------------------------
# Load packages
suppressPackageStartupMessages({
    library(TSENAT)
    library(ggplot2)
    library(SummarizedExperiment)
})

# Setup random seed for reproducibility
set.seed(42)


## ----prepare-example-data-----------------------------------------------------
# Load example dataset (lazy-loaded as 'readcounts' by default)
# This includes SALMON preprocessing outputs:
# - readcounts: Salmon NumReads (raw fragment counts with EM-based ambiguity resolution)
# - tpm: Salmon TPM matrix (length and library-normalized expression)
# - effective_length: Salmon EffectiveLength vector (accounts for fragment length distribution)
data(readcounts)

readcounts <- as.matrix(readcounts)

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

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

## ----initialize-analysis-object-----------------------------------------------
## Configure analysis parameters first (best practice: fail-fast principle)
## This validates all parameters before object creation
config <- TSENAT_config(
  sample_col = "sample",
  condition_col = "condition",
  subject_col = "paired_samples",
  q = seq(0, 2, by = 0.05),
  nthreads = 2,
  paired = TRUE,
  control = "normal"
)

## Build a complete `TSENATAnalysis` object 
analysis <- build_analysis(
  config = config,
  readcounts = readcounts,
  metadata = metadata_df,
  tx2gene = gff3_file, 
  tpm = tpm,
  effective_length = effective_length
)


## ----filter-validate-main-----------------------------------------------------
# Filter lowly-expressed transcripts
analysis <- filter_analysis(analysis, stringency = "medium")

## ----tsallis-entropy-computation, message = FALSE, warning = FALSE------------
# Set random seed for reproducible bootstrap CI calculations
set.seed(12345)

# Compute diversity
analysis <- calculate_diversity(
    analysis,
    norm = TRUE
)

# Extract and display diversity results
diversity <- results(analysis,
    type = "diversity",
    n_genes = 4,
    sample = "SRR14800481")

print(diversity)

## ----kable-tsallis-calc-sequence, echo = FALSE, results = "asis"--------------
# Extract diversity results for three specific q values using the results() accessor
# Show entropy values for ONE sample across three q values (0.05, 1, 2)
# Note: q=0 (richness) is computed deterministically outside bootstrap to avoid discrete data artifacts

q_values_to_show <- c(0, 0.5, 1, 1.5, 2)
sample_data <- NULL
first_sample_name <- NULL

for (q_val in q_values_to_show) {
    # Use tryCatch to gracefully handle q-values that don't exist
    div_q <- tryCatch({
        results(analysis, type = "diversity", q = q_val, format = "se")
    }, error = function(e) {
        NULL
    })
    
    if (!is.null(div_q)) {
        # Extract matrix from SummarizedExperiment
        if (is(div_q, "SummarizedExperiment")) {
            mat <- assay(div_q)
        } else if (is.matrix(div_q)) {
            mat <- div_q
        } else {
            next
        }
        
        # Get first 4 genes and first sample only
        first_sample_idx <- 1
        first_sample_name <- colnames(mat)[first_sample_idx]
        mat_subset <- mat[1:min(4, nrow(mat)), first_sample_idx, drop = FALSE]
        
        # Add values for this q with descriptive column names
        if (is.null(sample_data)) {
            sample_data <- data.frame(Gene = rownames(mat_subset))
        }
        
        col_name <- paste0("q=", sprintf("%.1f", q_val))
        sample_data[[col_name]] <- as.numeric(mat_subset)
    }
}

if (!is.null(sample_data) && nrow(sample_data) > 0) {
    knitr::kable(
        sample_data,
        caption = paste0("**Table 1:** Tsallis entropy for first 4 genes across three diversity scales (sample: ", 
                         first_sample_name, ")"),
        digits = 5,
        row.names = FALSE
    )
} else {
  cat("Note: Diversity results not available for the requested q-values. ",
      "Ensure calculate_diversity() was run with appropriate q-value configuration.")
}

## ----fig-1-isoform-diversity-profiles, fig.width=10, fig.height=5.2, out.width="98%", fig.pos="H", fig.cap = "**Figure 1:** Isoform diversity profiles across entropic indices. Lines show normalized Tsallis entropy (0-1) for each sample, blue control and red treatment."----
# Plot overall q-curve
p_qcurve <- plot_diversity_spectrum(analysis)
print(p_qcurve)

## ----qc-entropy-outlier-detection, echo=TRUE----------------------------------
# Perform multi-q sample influence analysis via S4 wrapper
analysis <- calculate_m_estimator(
    analysis,
    loss_type = "huber",
    q_combine_method = "mean",
    influence_threshold = 0.75
)

# Extract results from metadata using accessor function
sample_qc <- metadata(analysis, "m_estimate_results")
print(sample_qc)

## ----display-qc-outlier-results, echo=FALSE, results="asis"-------------------
if (!is.null(sample_qc) && nrow(sample_qc) > 0) {
  # Round numeric columns for better readability
  sample_qc$Proportion_Affected <- round(sample_qc$Proportion_Affected, 3)
  sample_qc$Distance_from_Centroid <- round(sample_qc$Distance_from_Centroid, 4)
  
  # Select columns: exclude Robustness_Weight and Pair_ID
  cols_to_show <- c("Sample", "Condition", "Proportion_Affected", 
                    "Distance_from_Centroid", "Status")
  cols_to_show <- cols_to_show[cols_to_show %in% colnames(sample_qc)]
  
  knitr::kable(
      sample_qc[, cols_to_show],
      caption = "**Table 2 | Sample Influence Assessment via M-estimation.** Ranked by magnitude of resampling-based influence values. Columns: Sample identifier; condition; proportion affected; distance from centroid; outlier status. M-estimation identifies influential samples for sensitivity analysis.",
      digits = 4,
      row.names = FALSE
  )
}

## ----test-q-condition-interaction---------------------------------------------
# Scale-adaptive interaction test across q values using S4 wrapper
analysis <- calculate_sait(
    analysis,
    method = "gam",
    multicorr = "hochberg"
)
       
# Extract SAIT interaction results
sait_results <- results(analysis, type = "sait", rankBy = "pvalue", n = 10)

print(sait_results)

## ----display-interaction-results, echo = FALSE, results="asis"----------------
sait_res <- results(analysis, type = "sait")
if (!is.null(sait_res) && is.data.frame(sait_res) && nrow(sait_res) > 0) {
    display <- head(sait_res, 6)
    
    # Select only most informative columns
    key_cols <- c("gene", "p_interaction", "adj_p_interaction", 
                  "effect_size", "test_statistic", "model_converged", "heteroscedasticity_detected")
    key_cols <- key_cols[key_cols %in% colnames(display)]
    display <- display[, key_cols]
    
    # Round numeric columns (but NOT p-value columns to preserve significance)
    numeric_cols <- sapply(display, is.numeric)
    pval_cols <- grepl("^p_|^adj_p_", colnames(display))
    display[, numeric_cols & !pval_cols] <- round(display[, numeric_cols & !pval_cols], 4)
    
    # Format p-values to show more precision (scientific notation for very small values)
    for (col in colnames(display)) {
        if (is.numeric(display[[col]]) && grepl("^p_|^adj_p_", col)) {
            display[[col]] <- format(display[[col]], scientific = TRUE, digits = 3)
        }
    }
    
    knitr::kable(display, 
                caption = "**Table 3 | Leading genes with scale-dependent condition effects from generalized additive models.** We display the 6 genes with lowest adjusted p-values (*q* < 0.05). Hochberg step-up correction applied (FWER control); columns show gene identifier, effect size, test statistic, convergence status, and heteroscedasticity detection. Methods: GAMs with *q* and condition as smooth predictors [@S190].",
                col.names = c("Gene", "P-value", "Adj. P-value", "Effect Size", "Test Statistic", "Model Converged", "Heteroscedasticity"),
                digits = 4, row.names = FALSE)
}

## ----fig-2-scale-dependent-genes, fig.width=14, fig.height=12, out.width="98%", fig.pos="H", fig.cap = "**Figure 2:** Scale-dependent interaction analysis. GAM-identified genes showing significant q$\times$condition effects (Hochberg-adjusted p < 0.05)."----
# Plot q-curve profiles for the top 4 genes using the S4 wrapper
combined_plot <- plot_sait(
    analysis,
    n_top = 4
)

print(combined_plot)

## ----identify-isoform-switching-transcripts, echo=TRUE------------------------
# Multi-Q analysis: Test switching at different q values
analysis <- calculate_jis(
    analysis,
    q = c(0, 0.5, 1, 1.5, 2),
    nboot = 100,
    lm_p_threshold = 0.05,
    threshold = 90
)

## ----extract-switching-results, echo=TRUE-------------------------------------
# Retrieve gene switching tables
tables_result <- results(analysis, type = "switching_tables")

## ----display-switching-table, echo=FALSE, results='asis'----------------------
# Only display tables if tables_result was successfully created
if (exists("tables_result") && is.list(tables_result) && 
    !is.null(tables_result$gene_headers) && length(tables_result$gene_headers) > 0) {
  # Iterate over the prepared comparison tables and display top 2 genes
  for (gene_idx in 1:min(2, length(tables_result$gene_headers))) {
    # Print gene header with proper spacing and formatting
    if (knitr::is_latex_output()) {
      cat("\n\n\\subsection*{Gene: ", tables_result$gene_headers[gene_idx], "}\n\n", sep="")
    } else {
      cat("\n\n### Gene: ", tables_result$gene_headers[gene_idx], "\n\n", sep="")
    }
    
    comparison_data <- tables_result$comparison_tables[[gene_idx]]
    q_metadata <- tables_result$q_metadata[[gene_idx]]
    
    if (!is.null(comparison_data) && !is.null(q_metadata)) {
      # Create display labels for column names (for knitr::kable)
      q_values_available <- q_metadata$q_values_available
      # Filter to show only q=0.50, q=1.00, q=1.50
      q_indices <- which(sapply(q_values_available, function(qk) {
        q_num <- q_metadata$q_key_to_value[[qk]]
        q_num %in% c(0.50, 1.00, 1.50)
      }))
      q_values_to_show <- q_values_available[q_indices]
      
      col_display <- c(
        "Transcript",
        sapply(q_values_to_show, function(qk) {
          q_num <- q_metadata$q_key_to_value[[qk]]
          paste0("q=", sprintf("%.2f", q_num))
        }),
        "Direction Consistency"
      )
      
      # Select only these columns from comparison_data
      cols_to_keep <- c(1, q_indices + 1, ncol(comparison_data))
      comparison_data_filtered <- comparison_data[, cols_to_keep, drop = FALSE]
      
      # Render table with kableExtra
      table_output <- knitr::kable(comparison_data_filtered, 
                         digits = 3, 
                         col.names = col_display,
                         align = c("l", rep("c", ncol(comparison_data_filtered) - 2), "r"))
      print(table_output)
    }
    
    cat("\n")
  }
} else {
  cat("Multi-Q switching tables not available; skipping this section.\n")
}

## ----fig-3-transcript-jackknife-influence, echo=TRUE, results='hide', fig.width=12, fig.height=10, eval=TRUE, out.width="98%", fig.pos="H", fig.cap = "**Figure 3:** Transcript-level switching patterns via jackknife resampling. Rows q-values (0-2.0), columns transcripts. Red/blue indicates condition influence."----
# Generate multi-Q heatmaps for top 4 genes using S4 wrapper
plot_jis_delta(analysis, n_genes = 4)

## ----fig-4-transcript-abundance-heatmap, echo=TRUE, results='hide', fig.width=12, fig.height=8, eval=TRUE, out.width="98%", fig.pos="H", fig.cap = "**Figure 4:** Transcript abundance heatmap with hierarchical clustering. Rows are transcripts, columns are conditions. Blue low, red high expression."----
# Generate transcript abundance heatmap
plot_expression(analysis, top_n = 4, metric = "median")

## ----effect-size-multiq-------------------------------------------------------
# Computes pairwise information-theoretic distance (Tsallis divergence)
analysis <- calculate_divergence(analysis)

## ----divergence-assay-display, echo = FALSE, results="asis"-------------------
# Display the raw divergence computation results
if (exists("analysis")) {
    div_results <- results(analysis, type = "divergence")
    
    # Extract divergence data - format depends on how it's stored in divergence_results slot
    if (is.matrix(div_results) && nrow(div_results) > 0) {
        div_table <- div_results[1:min(4, nrow(div_results)), 1:min(3, ncol(div_results))]
        
        knitr::kable(
            div_table,
            caption = "**Table 5 | Pairwise Tsallis divergence estimates.** Divergence (distance) between conditions from Tsallis entropy framework. Columns: gene identifier; pairwise comparison; divergence value; confidence interval (95%). Ranked by magnitude.",
            digits = 5,
            row.names = TRUE
        )
    } else if (is.list(div_results) && length(div_results) > 0) {
        # If stored as list per q-value, extract first q-value results
        first_q_result <- div_results[[1]]
        if (is.matrix(first_q_result)) {
            div_table <- first_q_result[1:min(4, nrow(first_q_result)), 1:min(3, ncol(first_q_result))]
            
            knitr::kable(
                div_table,
                caption = "**Table 6 | Pairwise Tsallis divergence from list-based results.** Alternative representation of divergence metrics for leading genes. Columns: gene identifier; pairwise comparison; divergence value; uncertainty. Asymmetry ($D(A \\to B) \\ne D(B \\to A)$) reflects information-theoretic properties. Ordered by divergence magnitude.",
                digits = 5,
                row.names = TRUE
            )
        }
    }
}

## ----effect-size-merge, echo=TRUE---------------------------------------------
# Compute effect sizes
analysis <- calculate_effect_sizes(
    analysis,
    significance_threshold = 0.05,
    enrich_per_q_pattern = TRUE
)

# Extract top 6 genes by statistical significance
top_genes_result <- results(
    analysis, 
    type = "effect_sizes_divergence",
    top_n = 6,
    sort_by = "p_value_interaction"
)

print(top_genes_result)

## ----divergence-effect-sizes-table, echo=FALSE, results="asis"----------------
# Select key columns for display (effect sizes at key q-values, pattern, and divergence metrics)
display_cols <- c(
    "gene", "p_value_interaction", "slope_diff",
    "effect_size_D_q0_5", "effect_size_D_q1", "effect_size_D_q2",
    "per_q_pattern"
)

# Only keep columns that exist in the data
display_cols <- display_cols[display_cols %in% colnames(top_genes_result)]
top_genes_display <- top_genes_result[, display_cols, drop = FALSE]

# Rename columns for better display
col_names <- c(
    "Gene", "P-value", "Slope Diff",
    "D(q=0.5)", "D(q=1.0)", "D(q=2.0)",
    "Pattern"
)
colnames(top_genes_display) <- col_names[1:ncol(top_genes_display)]

# Format p-values with scientific notation to avoid displaying 0
if ("P-value" %in% colnames(top_genes_display)) {
    top_genes_display[["P-value"]] <- format(
        top_genes_display[["P-value"]],
        scientific = TRUE, digits = 3
    )
}

# Display the result as a formatted table
knitr::kable(
    top_genes_display, 
    caption = "**Table 7 | Top 6 Genes by Effect Size (Tsallis Divergence).** Genes ranked by interaction p-value; columns show divergence at key q-values (0.5, 1.0, 2.0) and inferred pattern type (rare-driven, abundant-driven, or balanced).",
    digits = 4,
    row.names = FALSE,
    format.args = list(scientific = TRUE)
)

## ----effect-size-plot-multiq, fig.width=10, fig.height=6, out.width="98%", fig.align='center', fig.pos="H", fig.cap='**Figure 5:** Distribution of Tsallis divergence effect sizes across genes. X-axis shows divergence values; y-axis shows gene density. Threshold line indicates substantial divergence (D > 0.1) distinguishing signal genes from background.'----
# Visualize the distribution of Tsallis divergence effect sizes across genes
plot_obj <- plot_divergence_distribution(
    analysis, 
    threshold = 0.1
)

print(plot_obj)

## ----compare-q-spectra-plot, fig.width=10, fig.height=6, out.width='100%', fig.align='center', fig.cap='**Figure 6:** Q-spectrum curves for top genes by significance. Each line represents divergence as a function of q-value; shape reveals biological pattern (rare-driven, balanced, or abundant-driven).'----
# Visualize q-spectrum curves for top 4 genes by significance
p_multi <- plot_divergence_spectrum(
    analysis, 
    n_genes = 4, 
    use_pvalue_ranking = TRUE
)

print(p_multi)

## ----visualize-q-spectrum, fig.width=10, fig.height=5, out.width='100%', fig.align='center', fig.cap='**Figure 7:** Global divergence spectrum across all genes. Aggregated q-spectrum pattern shows which abundance scales (rare vs. abundant isoforms) are affected genome-wide.'----
# Visualize global divergence q-spectrum across all genes (aggregated)
# Plot is rendered directly by knitr (no file dependency)
p <- plot_divergence_spectrum(analysis)

print(p)

## ----session-info, results = 'markup'-----------------------------------------
sessionInfo()

