Test for Differential Splicing Between Groups

Performs differential splicing analysis by comparing orthologous splice junction usage between two groups of species while accounting for phylogenetic relationships. The function fits a statistical model, and identifies differentially used orthologous splice junctions.

Usage

runDS(
  data,
  se,
  minSpecies = "max",
  minCounts = 1,
  methodName = "phyloDEOU",
  alpha = 0.1,
  batch = NULL,
  nbThreads = 1,
  nbSimuEVE = 5000,
  seed = 42,
  ...
)

Arguments

data

List. Experiment object from newExperiment containing:

sampleInfo: data.frame with sample metadata including sampleID, speciesID, and group columns.
speciesTree: phylo object representing the species phylogeny (required for phylogenetic methods).

se

SummarizedExperiment. Orthologous splice junction count data from findOSJ with size factors computed by estSizeFactors. Must contain sizeFactor column in colData.

minSpecies

Named numeric vector or "max". Minimum number of species per group required to retain an orthologous splice junction. Format: c("1" = n1, "2" = n2) where n1 and n2 are the thresholds for groups 1 and 2 respectively. Default: "max", resolved to the number of species available in each group (i.e. an orthologous splice junction must be present in every species of both groups, the most stringent setting). Provide an explicit named vector to override, or set to NULL to skip species filtering entirely.

minCounts

Numeric. Minimum average normalized counts required in both splicedIn and splicedOut states to retain an orthologous splice junction. Set to NULL to skip count filtering. Default: 1.

methodName

Character. Statistical method for differential testing. One of:

"phyloDEOU" (default): Regularized phylogenetic (Ornstein-Uhlenbeck, OU) linear model. Recommended for most analyses.
"phyloDEBM": Regularized phylogenetic (Brownian Motion, BM) linear model.
"phylolmOU": Phylogenetic (OU) linear model.
"phylolmBM": Phylogenetic (BM) linear model.
"limma": Standard limma.
"limmaCorr": limma with duplicateCorrelation to account for intra-species correlation.
"DESeq2": Standard DESeq2.
"EVE": Phylogenetic ANOVA model.

alpha

Numeric. Significance threshold for Benjamini-Hochberg adjusted p-values. Orthologous splice junctions with padj <= alpha are flagged as significant. Default: 0.1.

batch

Character or NULL. Column name in colData(se) containing batch information for batch effect correction. When specified:

For limma-based methods: batch is added to the design matrix formula.
For DESeq2: batch is added to the design formula.
For EVE: batch effects are removed using limma::removeBatchEffect before analysis.

Default: NULL (no batch correction).

nbThreads

Integer. Number of parallel threads for model fitting. Default: 1.

nbSimuEVE

Number of orthologous splice junctions simulated to compute empirical p-values (EVE only).

seed

Integer. Random seed used for the EVE null simulations (EVE only), for reproducibility. Default: 42.

...

Additional arguments passed to the underlying method functions. For phylogenetic methods (phyloDEOU, phyloDEBM, phylolmOU, phylolmBM), these are passed to phyloDE::phylolmFit. See phylolmFit for available options.

Value

A SummarizedExperiment (filtered subset of input) with additional mcols columns:

significant: Logical. TRUE if orthologous splice junction shows significant differential splicing (padj <= alpha).
padj: Numeric. Benjamini-Hochberg adjusted p-values.
pval: Numeric. Raw p-values from the statistical test.
DID: Numeric. Difference-in-differences estimate: the phylogenetically corrected group difference (log2 scale) in the within-sample log-ratio between the focal orthologous splice junction and its competing junctions. Positive values indicate higher relative usage in group 2.

For the regularized phyloDE methods (phyloDEOU, phyloDEBM), the underlying model fit returned by phyloDE::phylolmFit (a limma-compatible object, with limma::eBayes already applied) is also attached to S4Vectors::metadata(se)$fit, enabling the downstream phyloDE/limma diagnostics (see the section below). This slot is absent for the other methods.

Details

The analysis proceeds in three steps:

1. Species filtering (two-stage):

Stage 1: Retain orthologous splice junctions present in at least minSpecies species in both groups (AND logic).
Stage 2: Expand to include overlapping orthologous splice junctions present in at least minSpecies species in either group (OR logic).

2. Count filtering: Orthologous splice junctions must have average counts >= minCounts in both splicedIn and splicedOut states across all samples.

3. Statistical testing: The selected method tests for differential splicing by modeling the interaction between splice state (in/out) and group membership. Phylogenetic methods account for non-independence due to shared evolutionary history.

Batch Effect Correction

When samples come from different batches (e.g., sequencing runs, labs, dates), batch effects can confound the analysis. The batch parameter allows correcting for these effects. The batch column must be present in colData(se) and should be a factor or character vector identifying the batch for each sample.

Method Selection

Recommended: phyloDEOU (default)

Accessing the phyloDE model fit

When a regularized phyloDE method is used (phyloDEOU, phyloDEBM), the phyloDE::phylolmFit object (with limma::eBayes already applied) is stored in S4Vectors::metadata(se)$fit. Because it follows the limma MArrayLM interface, it can be passed directly to the downstream phyloDE/limma functions. The coefficient of interest is "group2" (i.e. coef = 2 with the default ~ group design):


se  <- runDS(data, se, methodName = "phyloDEOU")
fit <- S4Vectors::metadata(se)$fit

# limma diagnostics
limma::topTable(fit, coef = "group2")
limma::plotMD(fit)
limma::volcanoplot(fit, coef = "group2")

# phyloDE-specific diagnostics
phyloDE::getParameters(fit)
phyloDE::plotParameters(fit)