Null Simulation Study for Differential Splicing Methods

Generates synthetic datasets under the null hypothesis (no differential splicing, DID=0) with various levels of within- and between-species correlation. Since all data are simulated without effect, any significant orthologous splice junction is a false positive. This allows evaluation of type I error control across methods.

Usage

runSimu(
  data,
  se,
  simDir,
  minSpecies = "max",
  minCounts = 1,
  maxJunctions = NULL,
  lambda = c(0, 0.1, 0.25, 0.5, 0.75, 0.9),
  modelProcess = "OU",
  thalf = c(0.05, 0.1, 0.2, 0.5, 1),
  rep = seq_len(10),
  skipDataGeneration = FALSE,
  methods = c("limma", "limmaCorr", "DESeq2", "EVE", "phylolmBM", "phylolmOU",
    "phyloDEBM", "phyloDEOU"),
  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.

se

SummarizedExperiment. Orthologous splice junction count data from findOSJ with size factors computed by estSizeFactors. Used to estimate realistic simulation parameters (dispersion, mean counts).

simDir

Character. Directory path for saving simulation files. Created if it does not exist.

minSpecies

Named numeric vector or "max". Minimum species per group for filtering. Format: c("1" = n1, "2" = n2). 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). Provide an explicit named vector to override, or set to NULL to skip species filtering.

minCounts

Numeric. Minimum average normalized counts in both states. Default: 1.

maxJunctions

Integer or NULL. Maximum number of orthologous splice junctions to use for simulation. If the number of orthologous splice junctions after filtering exceeds this value, a random subsample is taken. Default: NULL (no limit).

lambda

Numeric vector. Within-species correlation, i.e. correlation among replicates of the same species. Sets the strength of the correlation: values range from 0 (independent replicates) to close to 1 (strongly correlated, species-structured replicates). Default: c(0, 0.1, 0.25, 0.5, 0.75, 0.9).

modelProcess

Character. Evolutionary model for simulation:

"OU": Ornstein-Uhlenbeck process with stabilizing selection (default). Half-life controlled by thalf.
"BM": Brownian Motion (neutral drift). thalf is ignored.

thalf

Numeric vector. OU phylogenetic half-life, $t_{1/2} = \log(2)/\alpha$, expressed relative to tree height. Sets the structure of the between-species correlation:

Small thalf (e.g., 0.05): star-like tree, species blocks (quasi-independent species)
thalf of 1 or more: correlation across species follows the tree topology

Default: c(0.05, 0.1, 0.2, 0.5, 1). Ignored when lambda = 0 or modelProcess = "BM". Note: the tree is passed to compcodeR as is, so thalf is relative to tree height only if data$speciesTree has unit height (root-to-tip distance of 1); rescale it beforehand otherwise.

rep

Integer vector. Replicate IDs for each parameter combination. Default: 1:10.

skipDataGeneration

Logical. If TRUE, skip data generation and load previously generated datasets. Useful for re-running analyses with different methods. Default: FALSE.

methods

Character vector. Statistical methods to evaluate. One or more of:

"limma": Standard limma (ignores phylogeny)
"limmaCorr": limma with duplicate correlation
"DESeq2": Negative binomial GLM (ignores phylogeny)
"EVE": Phylogenetic ANOVA model
"phylolmBM": Phylogenetic linear model with BM
"phylolmOU": Phylogenetic linear model with OU
"phyloDEBM": Phylogenetic limma with BM consensus
"phyloDEOU": Phylogenetic limma with OU consensus (recommended)

Default: all methods.

batch

Character or NULL. Column name in colData(se) for batch correction during analysis. Default: NULL.

nbThreads

Integer. Number of parallel threads. Default: 1.

nbSimuEVE

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

seed

Integer. Random seed used for reproducibility, both when subsampling to maxJunctions and for the EVE null simulations. The same seed is used for both. Default: 42.

...

Additional arguments passed to the underlying method functions. For phylogenetic methods, these are passed to phyloDE::phylolmFit.

Value

A data.frame with simulation results metadata, one row per parameter-method combination:

rep: Replicate ID
modelProcess: Evolution model ("BM" or "OU")
lambda: Within-species correlation
thalf: OU half-life (relative to tree height)
method: Statistical method name
outputResFile: Path to TSV file with pval and DID columns

Details

The simulation workflow has three steps:

1. Parameter Estimation: Estimates realistic parameters from the input data:

Gene-specific dispersions via DESeq2
Mean counts per orthologous splice junction (splicedIn and splicedOut separately)
Sample structure (species, groups, replicates)

2. Data Generation: Generates synthetic counts using compcodeR::generateSyntheticData with:

Phylogenetic Poisson log-normal model with DESeq2-estimated dispersions
Within- and between-species correlation structure (controlled by lambda and thalf, respectively)
No differential splicing effect (effect.size = 1 for all genes)
Independent simulation of splicedIn and splicedOut counts (see Liehrmann et al., 2026 for details on the simulation framework)

3. Analysis: Applies each specified method to simulated datasets. Since DID=0 by design, significant results represent false positives.

Expected behavior:

At lambda ~ 0: All methods should have FPR close to alpha
At lambda > 0: Non-phylogenetic methods (limma, limmaCorr, DESeq2) and unregularized phylogenetic methods (phylolm, EVE) may show inflated FPR

Parameter Grid

The function generates a full factorial design of lambda × thalf × rep, except that when lambda = 0, only the first thalf value is used (since thalf is irrelevant when replicates are independent).