--- title: "TSENAT: Tsallis Entropy Analysis Toolbox" author: "Cristobal Gallardo " date: "`r Sys.Date()`" output: BiocStyle::html_document: toc: true toc_float: true toc_depth: 2 bibliography: TSENAT.bib suppress-bibliography: true vignette: | %\VignetteIndexEntry{1. TSENAT: Tsallis Entropy Analysis Toolbox} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r 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)) ``` # Introduction TSENAT (Tsallis Entropy Analysis Toolbox) allows the modelization and quantification of isoform-usage complexity in RNA-seq data using Tsallis entropy, a scale-dependent diversity measure distinct from differential abundance tools (DESeq2, edgeR) and differential usage tools (DEXSeq, DRIMSeq). By tuning the entropy index parameter (q), TSENAT enables to examine transcriptome heterogeneity at different scales: rare variants (low q) or dominant isoforms (high q). This package enables computing Tsallis entropy and Tsallis divergence from transcript-level abundance estimates, comparing measures between groups, and visualizing scale-dependent differences via q-curves. ## Motivation Common RNA-seq tools focus on either total gene abundance changes or individual transcript usage shifts, yet genes often remodel their isoform landscapes without changing overall expression-a phenomenon missed by these approaches. TSENAT makes possible to quantify this isoform-usage diversty directly via Tsallis entropy, allowing to capture biological signals invisible to abundance- or proportion-based summaries. In this guide, we demonstrate the complete workflow: preprocessing transcript counts, computing entropy across entropic indices, testing for between-group differences, and visualizing scale-dependent complexity. ## High-level workflow 1. **Load and preprocess** transcript counts and filter low-abundance transcripts. 2. **Compute diversity** (Tsallis entropy) and **divergence** (Tsallis divergence). 3. **Test for differences** in entropy patterns between groups across entropic indices. 4. **Visualize results**. Design assumptions: This guide assumes paired or longitudinal designs with >=6-8 samples per group for adequate power to detect entropy shifts while controlling false discovery rate. Smaller sample sizes may be underpowered to detect subtle isoform complexity changes. ## Installation The easiest way to install TSENAT is through Bioconda: ```{r install-bioconda, eval = FALSE} # Install from Bioconda # conda install -c bioconda r-tsenat ``` To install the latest development version from GitHub: ```{r install-github, eval = FALSE} remotes::install_github("gallardoalba/TSENAT") ``` ## Quick Start For users who want to see results quickly, here is a minimal workflow: ```{r 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) ``` # What is Entropy and Tsallis Entropy? ## Historical Foundation and Intellectual Progression The concept of entropy originated in Claude Shannon's landmark 1948 paper on information theory [@I051], which established that information content could be quantified mathematically. Shannon entropy became the foundation for understanding complexity across mathematics, physics, and biology—if you draw a transcript from a distribution, how predictable is the outcome? However, Shannon entropy treats all elements equally, regardless of their frequency: it weights rare and common elements identically. This limitation led mathematicians and physicists to explore generalized entropy families, most notably Rényi's parametric family and Tsallis entropy. Both introduced tunable parameters that allow sensitivity to shift between rare and abundant elements. Remarkably, despite appearing different mathematically, Tsallis and Rényi entropies can be unified within a coherent framework through generalized logarithmic and exponential functions [@I063]—both answer the fundamental question: What organizational scales matter? More recently, Hill numbers [@I021] provided a modern ecological framework that reinterprets all generalized entropy measures as "true diversity" of different orders. This formulation clarified a crucial insight: diversity questions have scale-dependent answers. The Hill numbers framework unified richness (q=0), Shannon entropy (q=1), Simpson/Gini (q=2), and higher-order generalizations under a single mathematical umbrella—formalizing the idea that rare species and dominant species reveal different ecological (or in our case, transcriptomic) truths. ## Why This Matters for RNA-seq: The Isoform Complexity Problem Standard RNA-seq analysis measures whether transcript abundance changes between conditions. But a critical complementary question remains underexplored: how does the diversity of isoforms change? A gene may show little change in total abundance while dramatically reshuffling its isoform repertoire—a phenomenon that current methods largely miss. Entropy quantifies precisely this: the complexity, richness, and balance of isoform heterogeneity. By measuring entropy across different biological scales (using the parameter q discussed below), researchers can detect whether changes are driven by shifts in rare isoforms or reorganization of dominant variants. The beauty of the Tsallis framework is that you obtain a complete picture of isoform heterogeneity by computing across a range of entropic indices—a "q-curve"—revealing which aspects of isoform organization change between conditions. ## Mathematical Foundation and Interpretation ### Tsallis Entropy: Definition and Intuition For a discrete probability vector p=(p₁,…,pₙ) representing isoform proportions within a gene, Tsallis entropy is defined as: $$S_q(p)=\frac{1-\sum_{i=1}^{n}p_i^q}{q-1}$$ This parametric family unites diverse entropy concepts under a single framework: - **Generalization**: Extends beyond Shannon entropy to capture scale-dependent phenomena - **Mathematical elegance**: Reduces to well-known diversity indices at specific entropic indices - **Practical flexibility**: Enables data-driven exploration across the full diversity spectrum ### The q Parameter: A Sensitivity Dial for Distribution Scales From an information theory perspective [@I051; @I002], entropy measures the uncertainty when drawing a single transcript from an isoform distribution. Higher entropy means the draw is less predictable (many similarly abundant isoforms), while lower entropy means one or a few isoforms dominate. The q parameter acts as a sensitivity dial that controls which aspects of the distribution become visible: - **q < 1** (e.g., 0.5): Emphasizes rare, low-abundance isoforms; useful for discovering cryptic or condition-specific variants. - **q = 1**: Recovers Shannon entropy; provides balanced sensitivity across all abundance scales. - **q > 1**: Emphasizes dominant, abundant isoforms; captures core expression architecture. This principle formalizes what cannot be implemented in classical Shannon analysis: the ability "not to set rare and common events on the same footing, as in standard statistics, but to enhance or depress them according to the parameter chosen" [@I005; @I007; @I009]. ### Biological Interpretation: Richness and Evenness The concept of "true diversity" [@I021] emphasizes that diversity decomposes into two independent components: species richness (how many distinct isoforms exist) and evenness (how evenly distributed they are across the population). Two genes can have identical Shannon entropy yet differ dramatically in isoform structure: one might be dominated by a single abundant isoform with many rare variants (revealing complexity at low q sensitivity), while another distributes transcripts equally across many isoforms (maintaining heterogeneity across all q values). This distinction is crucial: no single entropic index captures the complete complexity landscape. ## Why Multi-Scale Entropy Matters: Biological Contexts The multi-scale nature of Tsallis entropy makes it suited for exploring isoform complexity across diverse biological contexts. Different biological processes prioritize different organizational scales: **Isoform complexity as a biological signal**: Isoform switching—reorganization of the isoform landscape without necessarily changing total gene abundance—reflects strategic shifts in protein function driven by splicing regulation. Evidence from single-cell transcriptomics demonstrates that transcript-level complexity varies systematically across cell types and developmental states [@ISO010], validating that isoform heterogeneity is a genuine biological phenomenon rather than noise. By measuring information content at different entropic indices, researchers can detect patterns invisible to traditional transcript abundance measures alone: - Changes in rare isoform usage (revealed through low q sensitivity) might reflect exploratory or error-correction mechanisms - Shifts in dominant isoform selection (revealed through high q sensitivity) might reflect functional specialization or robustness demands - The full q-curve reveals whether cellular transitions involve wholesale reorganization or targeted adjustments TSENAT enables detection of these changes through entropy-based approaches, which capture whether complexity is increasing (diversity spreading across isoforms) or decreasing (consolidation onto dominant isoforms). The recent emphasis on information-theoretic approaches in computational biology [@I029; @I034] reflects broader recognition that complex biological systems encode information across multiple organizational scales. ## Beyond Classical Abundance Measures Traditional RNA-seq analysis focuses on fold-changes and differential abundance of specific transcripts. Tools like edgeR [@S195; @S199], DESeq2 [@S069], and DEXSeq [@S070] measure whether individual transcripts increase/decrease in expression or shift in relative proportion.Entropy-based approaches complement these methods by detecting complexity shifts —whether genes consolidate onto fewer dominant isoforms or maintain balanced distributions. TSENAT detects this reorganization answering a fundamentally different biological question: not "which transcripts change?" but "how is the isoform complexity changing?" ## Mechanistic Evidence: Disease and Evolution Peer-reviewed literature provides strong empirical support for entropy-based analysis in biological contexts: **Entropy and cancer heterogeneity**: Cancer cells accumulate genetic and epigenetic perturbations that systematically increase disorder in gene regulatory networks. Tarabichi and colleagues demonstrate that "Increased entropy of signaling (or gene interaction networks) has been well studied as a cancer characteristic: Network entropy increases along with cancer progresses" [@Y005]. Nijman's complementary analysis reveals the mechanism: "cancer-associated perturbations collectively disrupt normal gene regulatory networks by increasing their entropy. Importantly, in this model both somatic driver and passenger alterations contribute to 'perturbation-driven entropy', thereby increasing phenotypic heterogeneity and evolvability" [@Y006]. This framework elegantly explains observed cancer heterogeneity without requiring that every genetic change confers an advantage—some mutations contribute entropy directly through network disruption. Increased entropy in gene regulatory networks thus drives phenotypic heterogeneity and cellular plasticity, suggesting that transcript-level entropy captures similar organizational principles [@Y006]. # TSENAT Main Workflow We now demonstrate a complete workflow for detecting **isoform switching**-when cells reorganize their isoform landscape without necessarily changing total gene abundance. This often reflects strategic shifts in protein function driven by splicing regulation. In this workflow, you will: 1. Load transcript counts and sample metadata 2. Identify genes with significant scale-dependent isoform complexity changes (using Scale-Adaptive Interaction Tests) 3. Assess which individual transcripts drive those changes (using jackknife resampling) 4. Interpret results in the context of paired experimental designs We'll work with a paired experimental design where treated and control samples are linked, enabling detection of robust biological signals while controlling for subject-level variability. By the end, you'll know not just *which genes* undergo isoform switching, but *which transcripts* are responsible and *at which diversity scales* the switching occurs. ## Load data and metadata An example dataset is included for demonstration. ```{r workflow-setup} # Load packages suppressPackageStartupMessages({ library(TSENAT) library(ggplot2) library(SummarizedExperiment) }) # Setup random seed for reproducibility set.seed(42) ``` Now we will load the example dataset and associated metadata: ```{r 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") ``` ## Data preprocessing and filtering Next, we will create a `TSENATAnalysis` object containing a `SummarizedExperiment` data container (a standard Bioconductor class for storing high-throughput assay data along with associated metadata). You can read more about it in the [Bioconductor documentation](https://bioconductor.org/packages/release/bioc/html/SummarizedExperiment.html). We will use the GFF3.gz annotation file for extracting the transcript-to-gene mapping and building the analysis object. We'll also include sample metadata for improved analysis. Following Bioconductor best practices, we create the configuration FIRST, then pass it to the builder function (fail-fast principle): ```{r 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 ) ``` The `TSENATAnalysis` object contains a `SummarizedExperiment` with transcript-level counts and gene annotations extracted from the GFF3 file. To reduce noise and improve statistical power, we should filter out lowly-expressed transcripts. Expression estimates of transcript isoforms with zero or low expression might be highly variable [@S202; @S203]. For more details on the effect of transcript isoform prefiltering on differential transcript usage, see [this paper](https://doi.org/10.1186/s13059-015-0862-3). ```{r filter-validate-main} # Filter lowly-expressed transcripts analysis <- filter_analysis(analysis, stringency = "medium") ``` The `stringency = "medium"` parameter keeps transcripts present in >=50% of samples with TPM values above the median of mean gene TPM. This balances noise reduction with preservation of isoform diversity needed for meaningful entropy calculations. ## Compute Tsallis Entropy We now compute Tsallis entropy across your configured q-spectrum. Diversity measures provide a comprehensive framework for assessing isoform heterogeneity at multiple scales [@I007; @I010]. We normalize entropy to [0,1] range (`norm = TRUE`) for comparable cross-q assessment. ```{r 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) ``` ```{r 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.") } ``` With the q-spectrum we can produce q-curves at two levels: 1. **Overall diversity profile (Figure 1)**: Shows how entropy changes across q-values for all samples aggregated, colored by condition (control vs. treatment). 2. **Gene-specific q-curves (Figure 2)**: For genes with significant q×condition interactions, shows how each gene's diversity pattern differs by condition across the entropic spectrum. ```{r 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) ``` Diverging curves between groups indicate scale-dependent diversity differences: separation at low `q` implies differences in rare isoforms, while separation at high `q` signals differences in dominant isoforms. ## Quality Control: Sample Influence Assessment Before proceeding with group-level comparisons, we should assess whether individual samples exert disproportionate influence on our entropy estimates [@S115]. To address this, we employ a leave-one-out influence assessment combined with robust M-estimation [@S019] (iteratively re-weighted least squares with Huber loss). This approach quantifies how much each sample's removal affects the estimated location differences across the q-spectrum, providing a sample-level quality control metric independent of group assignment. Samples with high influence scores exert disproportionate leverage on entropy estimates and may warrant deeper investigation for technical quality issues, subject-specific outliers, or genuine biological signals that require careful interpretation. We configure three key parameters: `loss_type = "huber"` employs Huber M-estimation (robust to outliers with bounded influence), `q_combine_method = "mean"` averages influence metrics across the entire q-spectrum for a single score per sample, and `influence_threshold = 0.75` designates the percentile cutoff above which samples are flagged (such that the top 25% most influential samples are reviewed). ```{r 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) ``` ```{r 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 ) } ``` **Interpretation:** Samples with distance from centroid >1.5 or proportion of affected genes >0.85 are flagged for quality control review. Flagged samples may reflect genuine biological heterogeneity or technical artifacts requiring careful examination before proceeding with downstream analysis. ## Scale-Adaptive Interaction Tests Each gene produces a **q-curve**: a trajectory showing how entropy changes across entropic scales from rare (low `q`) to abundant (high `q`) isoforms. Our goal is to detect whether this curve differs between groups—in other words, whether the effect of `q` on entropy **depends on** which group a sample belongs to. This is captured statistically as a **q x condition interaction**: if significant, it reveals genes with condition-specific diversity patterns that change shape across the diversity spectrum. Scale-Adaptive Interaction Tests (SAIT) naturally formalize this multi-scale comparison through adaptive statistical frameworks, making them ideal for identifying such scale-dependent biological signals. We can test for these interactions using one of four methods, each with specific strengths: - GAM (generalized additive model): flexibly captures nonlinear q-response patterns via adaptive smoothing splines. Default method. - LMM (linear mixed models): models subject-level random intercepts to account for within-subject correlation in repeated q-ordered measurements. - GEE (generalized estimating equations): particularly useful for paired and longitudinal designs with repeated q-measures. - FPCA (functional principal component analysis): treats q-values as ordered functional data, implicitly capturing q-ordering structure. All methods automatically account for the AR(1) correlation structure inherent in Tsallis entropy measurements. ```{r 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) ``` ```{r 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). Benjamini-Hochberg correction applied; 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) } ``` **Interpretation:** TRUE in the Heteroscedasticity column indicates that the model satisfies homogeneity-of-variance assumptions across the q-spectrum; FALSE values suggest variance heterogeneity and warrant caution in result interpretation. Now we will plot the q-curve profile for the top genes identified by the scale-adaptive interaction test. ```{r 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 (Benjamini-Hochberg q < 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) ``` ## Transcript Switching Across Diversity Scales The SAIT interaction tests whether condition effects depend on which diversity scale (q-value) you examine. This section identifies **which individual transcripts** drive these scale-dependent patterns, revealing whether the same transcripts switch across all scales or whether different regulatory mechanisms dominate at rare versus abundant isoform scales. ### Two-Stage Analysis Approach This workflow combines two complementary statistical methods: 1. **Multi-*q* SAIT interaction test**: Tests whether the condition effect on diversity depends on *q* (overall pattern across diversity scales) 2. **Single-*q* jackknife resampling**: Tests which individual transcripts drive those patterns at each *q*-value, with robustness assessment **Delta influence** (from jackknife resampling [@S115]) quantifies *how each transcript's relative importance changes between normal and tumor conditions*, weighted by stability across bootstrap iterations (positive = more influential in normal condition; negative = more influential in tumor condition). ### Key Interpretation Questions - **Consistent switching**: Do the same transcripts show significant delta influence across all *q*-values? Suggests robust, scale-independent splicing shift. - **Mixed/scale-dependent switching**: Do different transcripts matter at different *q*-values? Suggests regulatory complexity where mechanisms differ by scale. - **Classification**: Does isoform reorganization involve the same transcripts across scales, or a layered process where different regulatory inputs dominate at rare vs. abundant scales? ```{r 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 ) ``` We configure three key parameters for the switching analysis: `nboot = 100` specifies 100 bootstrap resamples to estimate robust confidence intervals for delta influence values [@S115]; `lm_p_threshold = 0.05` filters genes to those with significant q$\times$condition interaction effects (p < 0.05) before assessing transcript switching; and `threshold = 90` designates transcripts with support $\geq$ 90% across bootstrap samples as robust switches. These parameters balance sensitivity (detecting switching) with specificity (avoiding false positives from noise). The tables below show the transcript switching patterns for the top genes with strongest q×condition interactions, automatically computed via bootstrap resampling. The `Direction Consistency` column classifies each transcript: "Consistent" indicates stable switching across entropic indices, while "Mixed" indicates scale-dependent switching. This reveals whether isoform shifts involve the same transcripts across scales or whether different entropic indices emphasize different transcripts. ```{r extract-switching-results, echo=TRUE} # Retrieve gene switching tables tables_result <- results(analysis, type = "switching_tables") ``` ```{r 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") } ``` ### Delta Influence Across Diversity Scales Visualize switching patterns for the top genes identified by the SAIT interaction test. This shows which transcripts are switching in genes with significant q × condition interaction effects. The heatmaps below display **jackknife delta influence** [@S232]-a resampling-based measure of how robustly each transcript's relative contribution changes between conditions-across different entropic indices and transcripts for each gene. Delta influence values are computed from bootstrap resampling iterations, providing robust, non-parametric estimates of transcript switching significance weighted by consistency across replicates. ```{r 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) ``` **Interpretation:** Rows represent q-values across the diversity spectrum (0 to 2.0, with low q emphasizing rare isoforms and high q emphasizing abundant isoforms). - **Red cells** indicate transcripts with strong positive delta influence in the **normal condition** (first condition)—meaning these transcripts became *more important* in normal samples compared to tumor samples. - **Blue cells** indicate negative delta influence in the normal condition (or equivalently, positive influence in the tumor condition)—meaning these transcripts became *less important* in normal samples but *more important* in tumor samples. - **Color intensity** reflects robustness across bootstrap resamples: bright red/blue = consistent signal, pale/white = noisy or uncertain signal. This robustness-weighted visualization reveals not just *which* transcripts switch, but *how reliably* they switch, allowing distinction between robust biological signals and artifacts of sampling variation. Finally we can visualize the isoform expression profiles. ```{r 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") ``` This heatmap reveals a critical blind spot in standard transcript quantification approaches: raw abundance measures reduce each transcript to a single number, completely obscuring the relational structure that defines isoform switching. TSENAT's entropy-based approach captures precisely what this heatmap reveals: not just which transcripts change in abundance, but how the diversity and balance of the isoform repertoire itself changes. ## Effect Size Analysis To understand which diversity scales show the most pronounced biological differences between groups, we will compute Tsallis divergence $D_q$ across the q-spectrum [@K001; @I002; @I003; @I004; @I059; @I061]]. This metrics allows to reveal whether group differences are driven by rare isoforms (low q), dominant isoforms (high q), or uniformly across scales. Values $D > 0.1$ indicate meaningful information-theoretic separation between conditions, revealing that isoform complexity patterns fundamentally differ between treatment groups across all q-dependent scales [@I060]. For two probability distributions $P$ and $Q$ representing isoform proportions in control and treatment conditions, Tsallis divergence is: $$D_q(P||Q) = \frac{1 - \sum_i p_i^q \cdot q_i^{1-q}}{q-1}$$ where $p_i$ and $q_i$ are the probability values at position $i$. ```{r effect-size-multiq} # Computes pairwise information-theoretic distance (Tsallis divergence) analysis <- calculate_divergence(analysis) ``` ```{r 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 ) } } } ``` ```{r 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) ``` ```{r 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)] # 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 ) ``` ### Interpretation of Pattern Types: - Rare driven: D(q=0.5) >> D(q=2.0). Low-abundance isoforms are condition-specific; treatment group preferentially expresses rare transcripts not seen in control. - Abundant driver: D(q=0.5) << D(q=2.0). High-abundance isoforms shift between conditions; treatment remodels the dominant transcript landscape without rare variants changing. - Balanced: D(q=0.5) ~= D(q=2.0). All isoforms shift proportionally; no preferential weighting to rare or abundant forms; suggests systematic rebalancing. The effect sizes reveal which genes show the most information-theoretic separation across the q-spectrum **between paired treatment groups**. Genes with Tsallis divergence D > 0.1 show substantial divergence, indicating that entropy distributions differ fundamentally between conditions across the full q-spectrum. **Biological Interpretation:** Following Tsallis entropy theory and information-theoretic principles [@I002; @I003; @I004; @S204], effect size quantification [@I029] genes with large D values are those where isoform complexity distributions differ qualitatively between conditions when examined at all sensitivity scales (q = rare -> abundant isoforms). For example, a gene might show high entropy at low q (emphasizing rare isoforms) in one condition but low entropy at high q (emphasizing abundant isoforms) in another, revealing scale-dependent regulatory mechanisms. These genes are candidates for investigation of condition-specific splicing architecture and dynamic isoform switching. ```{r 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) ``` ```{r 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) ``` Interpreting the Q-Spectrum Curve: - **Shape matters more than single value**: The **q-spectrum** curve is the complete biological story [@I021]. A flat curve (balanced) vs. declining curve (rare driven) vs. rising curve (abundant driven) encode fundamentally different mechanisms [@I059]. As described in the pattern table above, genes falling into each category reveal different regulatory strategies. - **q=1 (KL divergence)**: The middle point corresponds to ordinary Kullback-Leibler divergence [@K001; @I060], where the Tsallis divergence reduces to standard KL divergence in the limit q→1. This weights all isoforms equally. We can also plot the global divergence spectrum across all genes. This curve is the full biological signature of **isoform switching**-the pattern encodes which abundance scales (rare vs. abundant) are affected. ```{r 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) ``` Statistical Validation During Interpretation: 1. Bootstrap CI validity: Per-q CIs should narrow as q increases (higher q = aggregation effect, less variance). If CIs grow with q, check for data quality issues [@S232; @S115]. 2. Monotonicity check: By Tsallis entropy theory , entropy is monotone decreasing in q. Divergence should NOT show erratic increases with q. Small fluctuations are normal, but large spikes indicate numerical instability [@I002]. 3. Comparison with genome-wide patterns: Compute q-spectra for housekeeping genes (GAPDH, ACTB, etc.). These should show BALANCED patterns. If not, revisit normalization parameters [@S023]. # Appendices This main vignette is complemented by two comprehensive appendices: **Appendix A: Equivalence Validation - TSENAT vs SplicingFactory** Objetive: Validates that TSENAT's Shannon and Simpson entropy implementations are mathematically equivalent to SplicingFactory. **[View Appendix A](TSENAT_appendix_A.html)** **Appendix B: Non-Parametric Validation of Statistical Test Results via GAM and Rank-Based Methods** Objective: Validate q-value * group interaction detection results from statistical tests using complementary non-parametric modeling approaches. **[View Appendix B](TSENAT_appendix_B.html)** ## References 1. Adami, C. (2004). "Information theory in molecular biology." *Physics of Life Reviews*, 1(1), 3–22. https://doi.org/10.1016/j.plrev.2004.01.002 2. Alomani, G., & Kayid, M. (2023). "Further Properties of Tsallis Entropy and Its Application." *Entropy*, 25(2), 199. https://doi.org/10.3390/e25020199 3. Anastasiadis, A. (2012). "Entropy Properties and Multiple Tsallis Distributions." *Entropy*, 14, 174–176. https://doi.org/10.3390/e14020174 4. Bajic, D. (2024). "Information Theory, Living Systems, and Communication Engineering." *Entropy*, 26(5), 430. https://doi.org/10.3390/e26050430 5. Bartal, A., & Jagodnik, K. M. (2022). "Progress in and Opportunities for Applying Information Theory to Computational Biology and Bioinformatics." *Entropy*, 24(7), 925. https://doi.org/10.3390/e24070925 6. Benjamini, Y., & Hochberg, Y. (1995). "Controlling the false discovery rate: a practical and powerful approach to multiple testing." *Journal of the Royal Statistical Society B*, 57(1), 289–300 7. Cao, J., Packer, J. S., Ramani, V., Cusanovich, D. A., Huynh, C., Daza, R., ... & Shendure, J. (2017). "Comprehensive single-cell transcriptional profiling of a multicellular organism." *Science*, 357, 661–667. https://doi.org/10.1126/science.aam8940 8. Chakraborty, T. (2019). *Introductory Time Series Analysis*. Indian Statistical Institute, Kolkata 9. Chanda, P., Costa, E., Hu, J., Sukumar, S., Van Hemert, J., & Walia, R. (2020). "Information Theory in Computational Biology: Where We Stand Today." *Entropy*, 22(6), 627. https://doi.org/10.3390/e22060627 10. Chao, A., Chiu, C.-H., & Jost, L. (2010). "Phylogenetic Diversity Measures Based on Hill Numbers." *Philosophical Transactions of the Royal Society B*, 365(1558), 3599–3609. https://doi.org/10.1098/rstb.2010.0272 11. Cover, T. M., & Thomas, J. A. (2006). *Elements of Information Theory* (2nd ed.). Wiley-Interscience. 12. Efron, B., & Tibshirani, R. J. (1993). *An Introduction to the Bootstrap* (2nd ed.). Chapman and Hall 13. Erhard, F., Hense, B., Jafari, M., Siebourg-Polster, J., Dölken, L., & Zimmer, R. (2018). "Improved Ribo-seq puromycin target reliability using Bayesian nonparametrics." *Bioinformatics*, 34(12), 2096–2102 14. Ernst, M. D. (2004). "Permutation Methods: A Basis for Exact Inference." *Statistical Science*, 19(4), 676–685 15. Furuichi, S. (2006). "Information theoretical properties of Tsallis entropies." *Journal of Mathematical Physics*, 47, 023302. https://doi.org/10.1063/1.2165744 16. Gandrillon, O., Gaillard, M., Espinasse, T., Garnier, N. B., Dussiau, C., Kosmider, O., & Sujobert, P. (2021). "Entropy as a Measure of Variability and Stemness in Single-Cell Transcriptomics." *Entropy*, 24(1), e93532. https://doi.org/10.3390/e24010018 17. Gao, X., Tsai, S.-B., Liu, F., Pan, L., & Deng, Y. (2019). "Uncertainty Measure Based on Tsallis Entropy in Evidence Theory." *International Journal of Intelligent Systems*, 34(6), 1626–1647. https://doi.org/10.1002/int.22185 18. Golomb, R., Yoles, M., Fishilevich, S., Cohen, O., Savariego Peled, E., Dahary, D., ... & Pilpel, Y. (2026). "An Information Content Principle Explains Regulatory Patterns." *bioRxiv*, ahead of print. https://doi.org/10.1101/2026.02.19.706555 19. Hyndman, R. J., & Athanasopoulos, G. (2018). *Forecasting: Principles and Practice* (2nd ed.). OTexts. https://otexts.com/fpp2/ 20. Jost, L. (2006). "Entropy and Diversity." *Oikos*, 113(2), 363–375. https://doi.org/10.1111/j.2006.0030-1299.14714.x 21. Jose, J., & Lal, P. S. (2013). "Application of ARIMA(1,1,0) Model for Predicting Time Delay of Search Engine Crawlers." *Informatica Economică*, 17(4), 26–39. https://doi.org/10.12948/issn14531305/17.4.2013.03 22. Kullback, S., & Leibler, R. A. (1951). "On information and sufficiency." *Annals of Mathematical Statistics*, 22(1), 79–86. https://doi.org/10.1214/aoms/1177729694 23. Nijman, S. M. B. (2020). "Perturbation-driven entropy as a source of cancer cell heterogeneity." *Trends in Cancer*, 6(6), 454–462. https://doi.org/10.1016/j.trecan.2020.02.016 24. Phipson, B., & Smyth, G. K. (2010). "Permutation P-values should never be zero." *Statistical Applications in Genetics and Molecular Biology*, 9(1), Article 39 25. Ramírez-Reyes, A., Hernández-Montoya, A. R., Herrera-Corral, G., & Domínguez-Jiménez, I. (2016). "Determining the Entropic Index q of Tsallis Entropy in Images through Redundancy." *Entropy*, 18(8), 299. https://doi.org/10.3390/e18080299 26. Ré, M. A., & Azad, R. K. (2014). "Generalization of Entropy Based Divergence Measures for Symbolic Sequence Analysis." *PLoS ONE*, 9(4), e93532. https://doi.org/10.1371/journal.pone.0093532 27. Sason, I. (2022). "Divergence Measures: Mathematical Foundations and Applications in Information-Theoretic and Statistical Problems." *Entropy*, 24(5), 712. https://doi.org/10.3390/e24050712 28. Seweryn, M. T., Pietrzak, M., & Ma, Q. (2020). "Application of Information Theoretical Approaches to Assess Diversity and Similarity in Single-Cell Transcriptomics." *Computational and Structural Biotechnology Journal*, 18, 1830–1837. https://doi.org/10.1016/j.csbj.2020.05.006 29. Shannon, C. E. (1948). "A Mathematical Theory of Communication." *The Bell System Technical Journal*, 27(3–4), 379–423 30. Shiner, J. S., Emelyanova, N. A., & Gafarov, F. M. (2002). *Entropy and Entropy Generation: Fundamentals and Applications*. Kluwer Academic Publishers 31. Simpson, E. H. (1949). "Measurement of diversity." *Nature*, 163, 688. https://doi.org/10.1038/163688a0 32. Tarabichi, M., Antoniou, A., Saiselet, M., Pita, J. M., Andry, G., Dumont, J. E., ... & Maenhaut, C. (2013). "Systems biology of cancer: entropy, disorder, and selection-driven evolution to independence, invasion and 'swarm intelligence'." *Cancer Metastasis Reviews*, 32, 403–421. https://doi.org/10.1007/s10555-013-9431-y 33. Tsallis, C. (2017). "On the foundations of statistical mechanics." *European Physical Journal Special Topics*, 226, 1433–1443. https://doi.org/10.1140/epjst/e2016-60252-2 34. Van Erven, T., & Harremoes, P. (2014). "Rényi Divergence and Kullback-Leibler Divergence." *IEEE Transactions on Information Theory*, 60(7), 3797–3820. https://doi.org/10.1109/TIT.2014.2320500 35. Yulmetyev, R. M., Emelyanova, N. A., & Gafarov, F. M. (2004). "Dynamical Shannon Entropy and Information." *Physica A*, 341, 649–676. https://doi.org/10.1016/j.physa.2004.03.094 ## Session Information ```{r session-info, results = 'markup'} sessionInfo() ```