phyloDS Workflow

library(phyloDS)

Introduction

This vignette documents a complete phyloDS analysis, function by function, in the order you would normally call them, from experiment setup to the interpretation of a single differentially used orthologous splice junction (OSJ).

No code is executed here beyond loading the package: every chunk below shows the exact call you would make on your own data, with its arguments and output explained in prose, but none of it runs when this vignette builds. This keeps the vignette static and dependency-free.

For the full statistical model, see the accompanying manuscript.

Step A: Create Experiment Object

The first step is to create an experiment object that contains all metadata and analysis settings.

What does `newExperiment()` do?

This function initializes your analysis by:

Loading sample and species metadata
Loading the phylogenetic tree relating your species
Setting up directory paths for intermediate files (extracted DNA, MSAs)

Parameters explained:

Each argument, and every column expected in its table, is documented in detail on the newExperiment() help page. Briefly:

sampleInfo: one row per RNA-seq sample
speciesInfo: one row per species, including which group ("1" = background/outgroup, "2" = focal lineage) it belongs to
orthogroups: one row per gene, with its orthogroup and genomic coordinates
speciesTree: ultrametric phylogenetic tree (Newick)
dnaDir, msaDir: directories where extracted DNA sequences and alignments will be written

data <- newExperiment(
  sampleInfo  = "samples.txt",
  speciesTree = "speciesTree.nwk",
  speciesInfo = "species.txt",
  orthogroups = "orthogroups.txt",
  dnaDir      = "orthogroupDNA",
  msaDir      = "orthogroupMSA"
)

What’s in the returned object?

The data object is a list containing:

sampleInfo: data frame with sample metadata
speciesInfo: data frame with species metadata
speciesTree: phylogenetic tree (phylo object from ape)
orthogroups: GRanges object of orthologous gene coordinates
dirs: resolved dna/msa directory paths

Understanding group assignment

Group 2 represents the lineage effect you want to test; group 1 is the background/outgroup. In the primate dataset analyzed in the paper, group 2 is the hominoid clade and group 1 the other primates, and phyloDS identifies OSJs that are differentially used along the hominoid branch relative to the outgroup.

Step B: Calculate Splice Junction Coverage

What does `coverageSJ()` do?

Reads BAM alignment files for each sample
Identifies splice junctions from CIGAR strings
Counts reads supporting each junction across samples and genes

Genes are processed independently, so nbThreads sets how many genes are handled in parallel.

cov <- coverageSJ(data, nbThreads = 25)

cov is a data.frame with one row per splice junction observation:

geneID: gene identifier
sampleID: sample where the junction was detected
donor / acceptor: genomic positions of the donor/acceptor splice sites
splicedIn: number of reads supporting this junction

Step C: Extract DNA Sequences and Alignments

orthogroupDNA() extracts, per orthogroup, each species’ genomic DNA sequence and writes it to a FASTA file; orthogroupMSA() then aligns each orthogroup’s FASTA with MAFFT (--auto by default), establishing the multiple sequence alignment (MSA) coordinate system that findOSJ() (Step D) uses to identify orthologous splice junctions across species. Both functions write files to disk rather than returning an R object.

orthogroupDNA() parallelizes across species: each worker loads one species’ genome once and extracts all of that species’ genes, so nbThreads is useful up to the number of species. orthogroupMSA() parallelizes across orthogroups and, within each one, across MAFFT’s own threads: with nbThreads = 15 and nbThreadsByOrthogroup = 3 below, 5 orthogroups are aligned at once, each using 3 MAFFT threads.

orthogroupDNA(data, nbThreads = 2)
orthogroupMSA(data, nbThreads = 15, nbThreadsByOrthogroup = 3)

Scaling up: single node vs. HPC cluster

orthogroupMSA() already parallelizes MAFFT across the cores of a single machine (nbThreads, nbThreadsByOrthogroup), which is enough for small orthogroup sets. For a real, genome-wide dataset (thousands of orthogroups), building all MSAs this way can still take a long time even on a large multi-core node.

The companion repository aLiehrmann/phyloDS_reproducibility, which reproduces the phyloDS manuscript end to end, provides a Nextflow pipeline (code/pipeline/analysis/main.nf) that performs the exact same step (one MAFFT alignment per orthogroup, matching orthogroupMSA()’s logic) but as one independent job per orthogroup, distributed across multiple nodes of an HPC cluster (-profile slurm, pbs, or sge). This is the recommended route for aligning large numbers of orthogroups: it is not a replacement for orthogroupMSA()’s statistical logic, only an alternative execution strategy that trades single-node multi-core parallelism for cluster-wide, per-orthogroup parallelism.

Step D: Counting Orthologous Splice Junctions

What does `findOSJ()` do?

This step identifies orthologous splice junctions (OSJs) across species.

Map junctions to MSA coordinates: for each gene and sample, genomic junction positions (donor/acceptor) are converted to MSA coordinates. Junctions at the same MSA position across species are considered orthologous.
Build a count matrix: rows are orthologous junctions (unique MSA positions), columns are samples x 2 states (splicedIn: reads spanning the junction, splicedOut: reads spanning a competing, overlapping junction instead).

findOSJ() keeps every junction observed in at least one sample, including junctions detected in a single species. Cross-species support is required later, through the minSpecies filter of runDS() (Step G).

These two per-sample counts, $Y_{\text{in}}$ and $Y_{\text{out}}$ , are the raw ingredients of the statistical model used later by runDS() (see the manuscript for the full model): they are combined into a single within-sample log-ratio $\Delta$ per OSJ, which is what phyloDE actually tests for a group effect.

Orthogroups are processed independently, so nbThreads sets how many orthogroups are handled in parallel.

se <- findOSJ(data, cov, nbThreads = 25)

Understanding the output structure

findOSJ() returns a SummarizedExperiment with three main components:

colData: sample x state design. Each sample appears twice (splicedIn and splicedOut), providing the two counts used to compute $\Delta$ for every OSJ.
rowRanges: coordinates of orthologous splice junctions. start/end are MSA coordinates; for each species, <speciesID>.seqnames/.strand/ .donor/.acceptor/.geneID give the genomic position in that species. NA in donor/acceptor marks an alignment gap (junction absent in that species).
assay: the count matrix, read counts per junction and sample-state.

Step E: Estimate Size Factors

What does `estSizeFactors()` do?

estSizeFactors() computes one size factor per sample using DESeq2’s median-ratio method, estimated from OSJs with sufficient splicedIn coverage (minCounts, default 5) and applied identically to that sample’s splicedIn and splicedOut counts.

Because $\Delta$ is itself a within-sample log-ratio, this sample-specific factor cancels out of $\Delta$ and often has little effect on the test itself. It matters mainly when splicedIn or splicedOut is close to zero (e.g. a nearly constitutive or nearly absent OSJ), where that cancellation breaks down.

se <- estSizeFactors(se)

Step F: Simulation Study (Optional)

What does `runSimu()` do?

Evaluates the type~I error rate of different statistical methods under the null hypothesis (no differential usage, $\mathrm{DID}=0$ ) by:

Estimating parameters from your real data (DESeq2 dispersions and mean splicedIn/splicedOut counts)
Generating synthetic null data under a given phylogenetic scenario (BM or OU)
Running one or more methods on the simulated data
Computing false positive rates to compare methods

Parameters explained:

lambda: within-species correlation, i.e. correlation among replicates of the same species (0 = independent replicates, close to 1 = strongly correlated, species-structured replicates). It sets the strength of the correlation.
modelProcess: "OU" (Ornstein-Uhlenbeck, controlled by thalf) or "BM" (Brownian motion)
thalf: OU half-life $t_{1/2} = \log(2)/\alpha$ , relative to tree height. It sets the structure of the between-species correlation: small values give a star-like tree (species blocks, quasi-independent species), and values of 1 or more make the correlation follow the tree topology. Ignored for "BM" or lambda = 0. Note: thalf is relative to tree height only if data$speciesTree has unit height (root-to-tip distance of 1). runSimu() passes the tree to compcodeR as is, so rescale it beforehand (e.g. tree$edge.length <- tree$edge.length / max(ape::node.depth.edgelength(tree))).
rep: number of replicates per parameter combination
methods: which methods to benchmark
skipDataGeneration: reuse previously generated null datasets, useful when comparing additional methods

params <- runSimu(
  data,
  se,
  simDir       = "simulations/",
  lambda       = c(0, 0.5, 0.9),
  modelProcess = c("BM", "OU"),
  minCounts    = 5,
  rep          = seq_len(50),
  methods      = c("limma", "limmaCorr", "phylolmOU", "EVE", "phyloDEOU")
)

Performance evaluation

plotSimu() summarizes the empirical type~I error rate across the $(\lambda, t_{1/2})$ grid at a chosen significance level.

res <- plotSimu(
  params,
  alpha = 0.1
)

Key observation: as shown for the hominoids and mole-rats datasets in the manuscript (Figure 5), phyloDEOU keeps type~I error close to the nominal level across the simulated $(\lambda, t_{1/2})$ grid, while methods that ignore phylogeny (limma, limmaCorr) or estimate OU parameters per-junction without regularization (phylolmOU, EVE) inflate it in the evolutionary regimes estimated from the real data. The full simulation study and its figures are reported in the manuscript.

Step G: Test for Differential OSJ Usage

What does `runDS()` do?

This is the core statistical test:

Filters junctions by species representation (minSpecies, default "max": an OSJ is kept if it is observed in every species, or if it overlaps such an OSJ and is observed in every species of at least one group) and minimum counts (minCounts, default 1: mean normalized count of at least 1 for both splicedIn and splicedOut)
Fits the phylogenetic model described in the manuscript (default method: "phyloDEOU")
Tests the group effect ( $\mathrm{DID}_j$ ) for each OSJ
Adjusts p-values for multiple testing (Benjamini-Hochberg)
Flags significant OSJs based on alpha

se_X <- runDS(
  data,
  se,
  methodName = "phyloDEOU",
  alpha      = 0.1
)

What’s in the results?

mcols(se_X) gains four columns:

significant: TRUE/FALSE, is this OSJ differentially used at the chosen alpha?
padj / pval: adjusted / raw p-value
DID: difference-in-differences, the phylogenetically corrected group difference in $\Delta$ (log2 scale); positive = higher relative usage in group 2, negative = higher relative usage in group 1

Step H: Filter to Differentially Used Genes

What does `filterSignificantGenes()` do?

Once OSJs have been tested, it is often convenient to work with whole genes (orthogroups) that contain at least one significant OSJ, rather than the full set of tested OSJs. filterSignificantGenes() subsets se_X to keep every OSJ belonging to such a gene.

se_X <- filterSignificantGenes(se_X)

Step I: Predict Splice Site Strength with `OpenSpliceAI`

What does `predictDID()` do?

predictDID() provides an orthogonal, sequence-only line of evidence for the OSJs found differentially used in Step G. For each species, it predicts donor and acceptor splice-site probabilities from genomic sequence alone, using pre-trained OpenSpliceAI models, and summarizes them into a predicted DID:

$\mathrm{predDID} = \log_2\!\left(\frac{\bar p_2}{1-\bar p_2}\right) - \log_2\!\left(\frac{\bar p_1}{1-\bar p_1}\right)$

where $\bar p_1$ , $\bar p_2$ are the mean predicted probabilities across species in groups 1 and 2. A positive predDIDdonor/predDIDacceptor means the site is predicted stronger in group 2 from sequence alone; concordance with the RNA-seq-based DID from Step G supports a cis-regulatory basis for the shift (see the manuscript’s TRIM28 example).

Parameters explained:

predDir: directory used to cache per-gene predictions (reused on re-runs)
refSpecies: one of "Homo_sapiens", "Mus_musculus", "Danio_rerio", "Arabidopsis_thaliana" – determines which pre-trained models are used
models: model replicate(s) to average over

Note: this step calls the external openspliceai CLI on the DNA sequences extracted in Step C.

Genes are processed independently, so nbThreads sets how many genes run through openspliceai at once. On the primate dataset, this step was run through the same companion Nextflow pipeline as orthogroupMSA() (see Step C’s Scaling up section), with nbThreads = 20.

se_X <- predictDID(
  data       = data,
  se         = se_X,
  predDir    = "predictions/",
  refSpecies = "Homo_sapiens",
  models     = c("model_10000nt_rs10", "model_10000nt_rs11"),
  nbThreads  = 20
)

predDIDdonor/predDIDacceptor are picked up automatically by showOSJ() (Step K, below) when present, and by exportResults() (Step J) to render OpenSpliceAI tracks in the exported IGV session.

Step J: Annotate Junctions with Exons and Export Results

What does `annotateOSJ()` do?

Projects OSJs to a reference species’ genome coordinates
Identifies flanking exons from that species’ GTF annotation (donor and acceptor), and any fully skipped exon

annots_hs <- rtracklayer::import("genomes/Homo_sapiens.GRCh38.113.gtf")

X_se_hs <- annotateOSJ(
  data       = data,
  se         = se_X,
  refSpecies = "Homo_sapiens",
  refAnnots  = annots_hs
)

mcols(X_se_hs) gains:

donorExon / acceptorExon: exon ID(s) whose boundary matches the donor/acceptor splice site
skippedExon: exon(s) fully contained within the OSJ span
Multiple exon IDs (separated by ;) indicate alternative splicing complexity; NA means a novel junction, not matching any annotated exon boundary

`exportResults()`: IGV-ready output

For interactive exploration beyond this vignette, exportResults() writes, per orthogroup, the MSA as a BAM track, OpenSpliceAI predictions as BigWig tracks (if predDir from Step I is supplied), a color-coded GTF of DID results, and species-specific annotation tracks, plus an IGV session file that loads them all with a single click.

Orthogroups are processed in parallel (nbThreads), and within each orthogroup the per-species tracks are themselves processed in parallel (nbThreadsBySpecies).

Important: call exportResults() on se_X, not on X_se_hs from annotateOSJ() above. annotateOSJ() projects OSJ coordinates onto a single reference species’ genome, so X_se_hs no longer carries the per-species genomic coordinates that exportResults() needs to build its tracks.

exportResults(
  data                = data,
  se                  = se_X,
  outputDir           = "results/",
  predDir             = "predictions/",   # from predictDID(), Step I; NULL to skip OpenSpliceAI tracks
  nbThreads           = 16,
  nbThreadsBySpecies  = 8
)

See the manuscript’s TRIM28 example (Figure 3) for what such a session looks like once opened in IGV: OSJ-span, OpenSpliceAI strength, MSA, and per-species annotation tracks aligned on the same evolutionary coordinate system.

Step K: Explore a Single Result with `showOSJ()`

Rather than reproducing an IGV session, showOSJ() prints a compact, text-based report for a single OSJ, gathering everything computed so far (reference coordinates, runDS statistics, exon annotation, and predictDID scores when available) – a cheap way to inspect individual candidates.

id <- rownames(X_se_hs)[which(mcols(X_se_hs)$significant)[1]]
showOSJ(id, X_se_hs)

The “Predicted DID” section is shown only when predDIDdonor/ predDIDacceptor columns are present in se (Step I).

Introduction

Step A: Create Experiment Object

What does newExperiment() do?

Parameters explained:

What’s in the returned object?

Understanding group assignment

Step B: Calculate Splice Junction Coverage

What does coverageSJ() do?

Step C: Extract DNA Sequences and Alignments

Scaling up: single node vs. HPC cluster

Step D: Counting Orthologous Splice Junctions

What does findOSJ() do?

Understanding the output structure

Step E: Estimate Size Factors

What does estSizeFactors() do?

Step F: Simulation Study (Optional)

What does runSimu() do?

Parameters explained:

Performance evaluation

Step G: Test for Differential OSJ Usage

What does runDS() do?

What’s in the results?

Step H: Filter to Differentially Used Genes

What does filterSignificantGenes() do?

Step I: Predict Splice Site Strength with OpenSpliceAI

What does predictDID() do?

Parameters explained:

Step J: Annotate Junctions with Exons and Export Results

What does annotateOSJ() do?

exportResults(): IGV-ready output

Step K: Explore a Single Result with showOSJ()

What does `newExperiment()` do?

What does `coverageSJ()` do?

What does `findOSJ()` do?

What does `estSizeFactors()` do?

What does `runSimu()` do?

What does `runDS()` do?

What does `filterSignificantGenes()` do?

Step I: Predict Splice Site Strength with `OpenSpliceAI`

What does `predictDID()` do?

What does `annotateOSJ()` do?

`exportResults()`: IGV-ready output

Step K: Explore a Single Result with `showOSJ()`