Splicing QTLs From Junction And Isoform Counts

Why Splicing QTLs Need Their Own Phenotypes

A genetic variant can change which transcripts a gene produces without changing how much RNA the gene produces overall. For example, a splice junction may be used more often in one genotype group, while the total number of reads assigned to the gene is similar across groups. A splicing quantitative trait locus (sQTL) analysis targets this difference in transcript processing, not simply a difference in total expression.

The first step is to define a phenotype that represents transcript choice. Junction counts can quantify how often an RNA-sequencing read supports a particular splice junction relative to competing junctions. Transcript-level quantification can instead represent the fraction of a gene’s RNA assigned to each isoform. These measurements are related but not interchangeable: the meaning of a proportion depends on which transcript events and competing outcomes were used to construct it.

Defining a Phenotype From Junction Counts

Suppose a splice event has an inclusion path and an exclusion path. Let \(I\) be the number of informative reads supporting the inclusion path and \(E\) the number supporting the exclusion path for a donor. A common phenotype is percent-spliced-in, usually abbreviated PSI, expressed as a proportion between zero and one:

The denominator \(I+E\) is the informative count for that event; it is not necessarily the total number of reads from the gene. The event definition must identify which junctions count as inclusion and which count as exclusion. For events with more complex structure, the event-specific construction may differ, so the formula should not be applied without specifying the supporting evidence.

A transcript-level phenotype can be defined similarly. If a gene has quantified isoforms with abundances \(A_1,\ldots,A_K\), the proportion assigned to isoform \(k\) is

Provided the abundance estimates are on a comparable scale and the included isoforms form the intended set, these proportions sum to one. In practice, transcript quantification can be uncertain when isoforms share exons or reads cannot distinguish them. A phenotype therefore depends not only on the formula but also on the event or isoform definition, read assignment, and minimum evidence required.

Dosage Regression for a Splicing Phenotype

For a biallelic variant, let \(G_i\) denote donor \(i\)’s allele dosage, coded as 0, 1, or 2 copies of the effect allele. Let \(Y_i\) be that donor’s splicing phenotype, such as PSI. A direct model on the proportion scale is a linear regression:

Here, \(\boldsymbol{C}_i\) contains included covariates, \(\boldsymbol{\gamma}\) their coefficients, and \(\varepsilon_i\) the residual. The coefficient \(\beta_G\) is the expected change in PSI per additional effect allele, conditional on the covariates. If PSI is recorded as a proportion, a coefficient of 0.05 means an estimated increase of five percentage points per allele, not a five-percent relative increase.

Before using ordinary linear-regression inference, the model requires an appropriate mean structure and residual assumptions: conditional on genotype and covariates, the expected phenotype should be represented adequately by the stated linear form; donors should provide independent observations for this analysis; and residuals should have mean zero and suitable variance behavior. Classical small-sample tests and intervals additionally rely on an appropriate residual distribution, commonly normal errors. If donors are related or measurements are repeated, the independence assumption is not appropriate for a simple regression, and the dependence must be modeled. The phenotype’s bounds also matter: a linear fit can predict values below zero or above one, especially when observations cluster near a boundary.

An alternative is to transform an interior proportion with the logit function, \(\operatorname{logit}(Y)=\log\{Y/(1-Y)\}\), and regress the transformed phenotype on dosage and covariates. This keeps the modeled response on an unbounded scale, but changes the effect interpretation: the genotype coefficient is a change in log odds of inclusion relative to exclusion per allele. Exact zeros or ones have no finite logit, so a transformation requires a prespecified, defensible treatment of boundary values rather than an arbitrary adjustment chosen after seeing the association.

Logit-Scale Effects and Isoform Composition

A logit model changes the scale, not the biological target. For an interior PSI value \(Y_i\), a model of the form

makes \(\beta_G\) an additive change in log odds of inclusion relative to exclusion per allele. Exponentiating it, \(\exp(\beta_G)\), gives the multiplicative change in those odds per allele, under this model. It is not an odds ratio for disease, and it is not a change in total gene abundance.

Isoform proportions have a further constraint: they sum to one. If there are three isoforms, increasing the proportion of one or more necessarily reduces the proportion of another. Consequently, separate associations for every isoform are not independent pieces of evidence. A useful way to express this is with a contrast between two isoforms, such as the log ratio \(\log(q_1/q_2)\), when both proportions are positive. A dosage effect on that contrast describes a change in their relative representation. The contrast does not describe a change in total RNA.

The constraint follows directly from the definition. If \(q_1+q_2+q_3=1\), then for any changes \(\Delta q_1,\Delta q_2,\Delta q_3\) that preserve the total, \(\Delta q_1+\Delta q_2+\Delta q_3=0\). A rise of 0.10 in \(q_1\) must therefore be balanced by a combined fall of 0.10 in the other proportions. This is why an isoform-proportion coefficient is inherently relative to the other isoforms in the denominator.

What an sQTL Result Does—and Does Not—Say

A nonzero genotype coefficient supports an association between genotype and the defined splicing phenotype, conditional on the model covariates. Its biological interpretation is tied to the assay and event: for a junction PSI phenotype, the result concerns the relative support for the specified inclusion and exclusion paths; for an isoform proportion, it concerns the isoform’s share of the chosen transcript set. It does not by itself establish a causal mechanism, identify every affected transcript, or show that total gene expression changed.

Phenotype construction and model choice should be reported together. A clear result states the event or isoform, how its proportion was calculated, the genotype dosage coding, the modeled scale, the covariates, and the effect units. For instance, “each effect allele was associated with a 0.04 increase in junction PSI” is interpretable as a four-percentage-point increase in the specified junction-usage proportion. By contrast, “the variant increased expression by 0.04” confuses a splicing phenotype with abundance and leaves the scale unspecified.

Common Mistakes

• Calling PSI total expression. PSI is a fraction among defined informative evidence. A rise in PSI can occur while total gene abundance stays constant, rises, or falls.
• Ignoring the denominator. The same inclusion count can yield different PSI values depending on the competing exclusion evidence. Always state the event definition and denominator.
• Reading a proportion coefficient as a relative percent change. A linear-model coefficient of 0.05 is a five-percentage-point change per allele. It is not automatically a 5% relative change.
• Treating all isoform associations as independent. Isoform proportions sum to one, so a change in one component constrains the others. Interpret results relative to the specified set or contrast.
• Overinterpreting noisy proportions. A proportion based on little informative evidence can be unstable. Read support and phenotype construction should be examined alongside the association estimate.