P4S

Plant Splicing Site Secondary Structure How does RNA secondary structure stability around splice sites drive splicing types?

Getting splicing data: SpliceScape

🚧🚧🚧

Extracting splice site flanking regions

🚧🚧🚧

Accessing secondary structure stability

🚧🚧🚧

Statistical Analysis

Main objective

Does the distribution of ΔG change depending on the splicing type?

Statistical rationale

To properly address this question, a few considerations are important:

1. We aim to avoid relying solely on non-parametric tests. While robust, these tests do not provide insight into how distributions differ. Parametric approaches allow us to investigate differences in location, dispersion, and shape.

2. Our dataset is relatively large. In such cases, p-values can become overly sensitive, leading to statistically significant results that may not reflect biologically meaningful differences.

Distribution fitting approach

To model the data, we used the gamlss package in R.

The first step is to determine which distribution model best describes the observed data. For this, we used the fitDist function, which evaluates multiple candidate distributions and selects the best one based on an information criterion.

fitDist tests a wide range of distributions (e.g., Normal, Log-Normal, Weibull, Gamma, Skew-Normal), fitting their parameters (such as mean, variance, and skewness) to approximate the empirical data as closely as possible.

Model selection: AIC vs BIC

AIC (Akaike Information Criterion)

AIC penalizes model complexity based only on the number of parameters. It balances goodness-of-fit (via deviance) with model simplicity, but tends to favor more complex models, especially in large datasets where even small fluctuations may appear statistically meaningful.

Penalty term: AIC = Deviance + 2 × (number of parameters)

This penalty is constant and does not depend on sample size (N).

BIC (Bayesian Information Criterion)

BIC introduces a stronger penalty that depends on sample size. As N increases, the penalty for additional parameters becomes larger, favoring simpler and more parsimonious models.

Penalty term: BIC = Deviance + ln(N) × (number of parameters)

Because of this, BIC is generally more conservative than AIC, particularly for large datasets.

Sample Size (n)	Criterion	Rationale
n < 40	AICc (corrected AIC)	AIC and BIC tend to select overly complex models in very small samples. AICc introduces an additional penalty to correct for this bias.
40 < n < 500	AIC	AIC performs well in this range. It has a higher tendency to select models with more parameters, which helps capture important nuances in moderate samples.
500 < n < 1000	Transition	Both criteria usually agree. When they diverge, BIC starts to become more reliable in preventing noise from being interpreted as signal.
n > 1000	BIC	Mandatory. For large n, the AIC penalty (k = 2) becomes negligible. BIC, by using ln(n), maintains stricter control and ensures model consistency.

Defining the baseline model

Based on this framework, the following workflow is applied:

1. Identify the best-fitting model for the reference group (constitutive events) and confirm that it is appropriate through residual diagnostics (see below).

2. If the model is adequate, we then assess whether it can also be applied to alternative splicing events (restricted to groups with more than 100 events), ensuring consistency across analyses:

   2.1. For each type of alternative event, we identify the best-fitting model and compute the ΔAIC. We then evaluate whether models are transferable (see below).

   2.2. We perform residual diagnostics (see below).

Model transferability logic

• If step 2.1 indicates that models are identical or transferable, the constitutive model can be safely applied to the alternative events.

• If models are not clearly transferable, step 2.2 (residual analysis) becomes decisive:

  • If residual diagnostics are satisfactory, the constitutive model may still be applied, with the caveat that it is not the optimal fit.
  • If residual diagnostics are poor, models are considered non-transferable.

In such cases, this result is biologically informative, as it suggests differences in dispersion between groups. However, parametric comparisons are no longer appropriate, since the data cannot be described under a common model. Therefore, differences must be assessed using non-parametric tests.

---

Residual diagnostics

Residual analysis is based on four diagnostic plots:

A. Residuals vs Fitted Values (top-left)

This plot shows residuals as a function of fitted values.

What to look for:

• A uniform cloud of points or a vertically homogeneous spread
• Residuals symmetrically distributed around zero
• No visible patterns (e.g., funnel shape)

---

B. Residuals vs Index (top-right)

This plot shows residuals in the order they appear in the dataset.

What to look for:

• Complete randomness ("white noise")
• No trends, cycles, or clustering

This indicates the absence of ordering bias in the ΔG data.

---

C. Density Estimate (bottom-left)

This represents the distribution of residuals.

What to look for:

• A normal-like distribution centered around zero
• A well-distributed rug plot (no large gaps in sampling)

---

D. Normal Q–Q Plot (bottom-right)

This plot compares observed residuals to a theoretical normal distribution.

What to look for:

• Points lying closely along the diagonal line

Deviations from the line indicate departures from normality.

Interpreting ΔAIC for model transferability

To assess whether models are transferable, we use the scale proposed by Kass and Raftery (1995):

ΔSBC	Evidence	Recommended Action
0–2	Negligible	Models can be unified.
2–6	Positive (weak)	Models can be unified. The information loss is minimal compared to the gain in comparability.
6–10	Strong	Caution. Inspect the Q–Q plot: if diagnostics are good, unify; otherwise, report the difference.
>10	Very strong	Inspect the Q–Q plot, but models are substantially different. This is a biologically meaningful result.

Outcome

This step defines:

• Which model will be used for comparison
• For which splicing types parametric analysis is valid

All remaining cases are analyzed using non-parametric tests.

---

Comparison via simulation

To evaluate whether and how the ΔG distribution changes across splicing types, we use a simulation-based comparison approach. The idea is to assume the simplest possible scenario and test whether the data reject it.

We use the gamlss framework, which fits models describing how the distribution of the data changes as a function of an explanatory variable.

Model comparison strategy

1. First, we assess whether separating the distributions is justified:

• m0: a single distribution for all data
• m1: one distribution per splicing type

We then compare these models using:

Likelihood Ratio Test (LRT)

Based on maximum likelihood:

• Small p-value → strong improvement → distributions differ
  → Evidence that splicing type affects ΔG (stability)

• Large p-value → little improvement
  → No strong evidence of an effect

ΔAIC

Using the same rationale as before:

• ΔAIC > 10 indicates a substantial difference between models

---

Expected outcome

Does including splicing type significantly improve the model fit?

Parameter refinement

If there is significant evidence and ΔAIC indicates that incorporating splicing type improves the model, the next questions are:

• Which splicing types drive this difference?
• Which distribution parameters are responsible for it?

Refinement strategy

We compare:

• Baseline model (no splicing effect):
  μ and σ are constant

• Model with splicing effect:

  • deltaG ~ label → mean (μ) varies across splicing types
  • sigma.fo = ~ label → variability (σ) varies
  • nu.fo = ~ label → skewness (ν) varies
  • tau.fo = ~ label → kurtosis (τ) varies

---

Parameter interpretation

Parameter	Name	Primary Function	What to look for
μ (Mu)	Location	Defines the center of the distribution	Is the value biologically plausible?
σ (Sigma)	Scale	Defines the spread (variance)	Larger σ indicates higher variability
ν (Nu)	Skewness	Controls asymmetry	ν ≠ 1 (in many families) indicates asymmetry
τ (Tau)	Kurtosis	Controls tail behavior	High τ suggests heavy tails (more extreme values)

---

Biological interpretation

• If μ changes → alternative events have different mean stability
• If σ changes → alternative events are more/less heterogeneous
• If ν changes → distributions differ in asymmetry
• If τ changes → distributions differ in frequency of extreme values
• If multiple parameters change → indicates a structurally different distribution profile

---

Model evaluation (per parameter)

For each parameter, we repeat:

• Likelihood Ratio Test (LRT):
  Tests statistical significance
  → Does the more complex model significantly improve the fit?

• AIC (Akaike Information Criterion):
  Evaluates model quality while penalizing complexity
  → Is the added complexity justified?

Only parameters that are both statistically significant and improve model quality are retained.

---

Final step: model summary

For selected models, we use summary() to determine:

• Which parameters are significant
• For which splicing types the inclusion of parameters is justified

---

How to interpret the summary output

• Intercept (μ):
  Represents the estimated mean for the reference group (constitutive events)

• labelXXX (μ):
  Represents the difference relative to the constitutive group (not the absolute mean)

Non-parametric statistical analysis and plot generation

🚧🚧🚧