## ----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(TSENAT) library(ggplot2) library(SummarizedExperiment) library(dplyr) library(gridExtra) }) set.seed(42) # 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-srh-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 Scheirer-Ray-Hare Test.** Evaluates homogeneity of variance and normality for rank-based interaction test validity." ) } ## ----srh-rank-test-analysis, message=FALSE, results='hide'-------------------- # IMPORTANT: Create a copy of analysis before running SRH # This preserves the GAMM results in 'analysis' for later concordance comparison # Run SRH test for q * condition INTERACTION on the copy analysis <- calculate_srh( analysis, multicorr = "hochberg" ) # Extract ALL results (no filtering) for summary statistics srh_results_all <- results(analysis, type = "rank_test", rankBy = "pvalue") # Extract top significant genes for display srh_results <- results(analysis, type = "rank_test", rankBy = "pvalue", n = 20, filterFDR = 0.05) print(head(srh_results, n = 10)) ## ----srh-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(srh_results_all$effect_size_eta2, na.rm = TRUE) median_effect <- median(srh_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(srh_results_all), sum(srh_results_all$p_value < 0.05, na.rm = TRUE), sum(srh_results_all$adj_p_value < 0.05, na.rm = TRUE), sum(is.na(srh_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(srh_results_all$p_value < 0.05 & srh_results_all$effect_size_eta2 > 0.10, na.rm = TRUE) ) ) knitr::kable(rank_summary_stats, caption = "**Supplementary Table 15 | Scheirer-Ray-Hare Test Results Summary.** Rank-based nonparametric test of *q* × condition interactions. Compare with GAMM results (main vignette Table 3) to assess parametric vs. nonparametric method concordance.", label = "tab:srh-test-summary", align = c("l", "r")) # Show top genes results rank_sig <- srh_results_all[!is.na(srh_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 == "srh_paired") "Scheirer-Ray-Hare (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 Scheirer-Ray-Hare rank-based test.** 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)]) } } ## ----extended-fig-1-srh-q-curves, fig.width=12, fig.height=8, message=FALSE, out.width="98%", fig.pos="H", fig.cap="**Supplementary Figure 1 | SRH rank test results for scale-dependent isoform switching.** Q-curves for top genes identified by Scheirer-Ray-Hare rank-based test; compare visual patterns with GAMM results (main vignette Figure 2) to assess method agreement."---- # Plot top genes from SRH test (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 = srh_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 = "tab: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, df cols_to_show <- intersect(c("gene", "p_interaction", "adj_p_interaction", "effect_size", "test_statistic", "df_residual"), 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) }) } # Format df as integer if ("df_residual" %in% colnames(gam_display)) { gam_display$df_residual <- sapply(gam_display$df_residual, function(x) { if (is.na(x)) "NA" else as.integer(x) }) } # Rename columns for better readability colnames(gam_display) <- c("Gene", "p (interaction)", "Adjusted p", "Effect size", "Test statistic", "df") 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 SRH results (Supplementary Table 16) for method agreement assessment.", label = "tab:gam-top-genes", align = c("l", "r", "r", "r", "r", "r")) } } } else { cat("WARNING: p_interaction column not found in results\n") cat("Available columns:", paste(colnames(sait_results), collapse=", "), "\n") } ## ----compare-srh-gam-methods, message=FALSE----------------------------------- # Compute concordance analysis using new two-object API analysis <- calculate_concordance( analysis_sait = analysis_sait, analysis_rank = analysis, verbose = TRUE ) # Extract concordance results with list format for component display concordance_results <- results(analysis, 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) { knitr::kable(concordance_results$summary_table, caption = "**Supplementary Table 19 | Global Concordance Metrics: GAMM vs Scheirer-Ray-Hare Methods.** Quantifies agreement between parametric (GAMM) and nonparametric (SRH) approaches for detecting q × condition interactions.", align = c("l", "r")) } # Display agreement distribution if (!is.null(concordance_results$agreement_dist) && nrow(concordance_results$agreement_dist) > 0) { cat("\n") knitr::kable(concordance_results$agreement_dist, caption = "**Supplementary Table 20 | Agreement Distribution: Concordance by Gene Category.** Distribution of genes across categories: both methods significant, GAM only, SRH only, and neither method.", align = c("l", "r", "r")) } # 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) 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 GAM and SRH provide strongest evidence for q × condition effects."), align = c("l", "r", "r", "r", "r")) } } 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 SRH rank-based and GAMM approaches for detecting q × condition interactions.** Venn diagram or scatter plot showing overlap in genes identified as significant by each method; reveals parametric-nonparametric agreement patterns."---- # Create comparison plot using S4 wrapper # Concordance analysis results are stored in analysis metadata after calculate_concordance() p <- plot_concordance(analysis, verbose = TRUE) print(p) ## ----session-info-appendix-b-------------------------------------------------- sessionInfo()