Bayesian Reconstruction and Differential
Testing of Excised mRNA

BREM solves the SME reconstruction problem by representing samples as mixtures of SMEs sampled from a global distribution that is learned across samples (Fig. LABEL:fig:overview). SMEs are sequences of mRNA excisions, typically, but not limited to introns, that can be assigned to the same transcript (i.e., they do not overlap). Briefly, BREM models samples as mixtures of SMEs, which are themselves mixtures of junction reads. BREM learns the structure of SMEs, a global distribution over SMEs, a mapping between junctions reads and SMEs, and a sample-specific distribution of SMEs. A separate model is built for each gene, which includes $i=1\dots N$ samples that are collections of reads overlapping excision junctions. Let the set of unique excisions within a gene be $V$, which is indexed by $v\in\{1,\dots,|V|\}$. The $i^{th}$ sample consists of $1,\dots,J_{i}$ junction reads. The goals are to reconstruct the latent SMEs and assign sample reads to SMEs for subsequent differential testing of SMEs. BREM consists of two major components: (a) combinatorial model for SME structure and (b) a probabilistic model for SME admixture.

We develop a Gibbs sampling algorithm to fit our model that proceeds by first sampling parameter values from their priors. Then, for each latent variable $z$, we sample from their complete conditionals, i.e., the probability of $z$ given all other random variables in the model (see §LABEL:sec:supp_gibbs). To determine Gibbs sampling convergence, we used Relative Fixed-Width Stopping Rules (RFWSR) [13]. RFWSR sequentially checks the width of a confidence interval relative to a threshold (here, $\sigma=0.001$) based on the effective sample size. Most variables yield efficient updates (full derivations are provided in the supplemental materials §LABEL:sec:supp_gibbs); however, special consideration is required for $z$ and $b$.

Here, we define our generalized linear model (GLM) to compute differential SME usage using the fitted model parameters in BREM. We quantify differential SME usage based on the expression profile across all SMEs for $2$ groups of samples in each gene. The $z$ variables represent the mapping of an excision observation to a SME. Let the counts across all unique junction reads for sample $i$ be denoted $z_{i}$. Then, we can express $z_{i}$ as a Dirichlet-multinomial

$$z_{i1},\dots,z_{iJ_{i}}\sim DirMult\left(\sum_{j}z_{ij},\alpha\odot p_{i}\right)$$

where $p_{ij}=\frac{\exp(x_{i}\beta_{j}+\mu_{j})}{\sum_{k}\exp(x_{i}\beta_{k}+\mu_{k})}$. We set $\alpha\sim\gamma(1+10^{-4},10^{-4})$ to stabilize maximum likelihood estimation [58]. Finally, to test differential SMEs between two groups, we construct two models: (a) a DirMult GLM where we set $x_{i}=0$ for one group and $x_{i}=1$ for the other and (b) a DirMult GLM where all $x_{i}=0$. Differential SMEs are quantified by a likelihood ratio test with $K-1$ degrees of freedom, where $K$ is the number of SMEs.

To appropriately evaluate methods that compute transcript segments of varying size we define two measures based on excision matching. If the set of expressed transcripts is known a priori, computed excisions can be evaluated using variants of homogeneity scores and partial precision and recall [51], however that is not the case here. Since we can compute excisions from exons, but not vice versa, we define transcript segments by the set of their component excisions. Reconstructed transcript segments may include exons from disparate transcripts. With a known reference, we compute the number of excisions that appear in any reference transcript and normalize by the total number of excisions. We define the $k^{th}$ computed transcript segment $T_{k}$ as a subset of excisions, or, $T_{k}\subseteq\{1,\dots,V\}$. Let the set of reference transcripts be $T^{t}$ and the $v^{th}$ excision be $e_{v}$. The set $T^{t}$ either represents a simulated baseline or known experimental transcripts from a well characterized cell type. The partial homogeneity score (phs) for transcript $T_{k}$ in sample $i$ can be computed as

$$s_{i}^{phs}(T_{k})=max_{T\in T^{t}}\frac{\sum_{e_{j}\in T_{k}}\mathbbm{1}\left[e_{j}\in T\right]}{|T_{k}|}$$

where $\mathbbm{1}[e_{j}\in T]$ is $1$ if $e_{j}$ matches an excision in $T$ and $0$ otherwise. Here, we consider two excisions $e_{v}$ and $e_{w}$ as matching if the donor and acceptor splice sites of $e_{v}$ are at most $6$ bases from the donor and acceptor splice sites of $e_{w}$. The score $s_{i}^{phs}$ enforces that excisions are sampled from the same true transcript and is normalized by the size of the computed transcript segment. We also define $\hat{s}_{i}^{phs}$, which normalizes computed transcript segments by the true transcript length.

$$\hat{s}_{i}^{phs}(T_{k})=max_{T\in T^{t}}\frac{\sum_{e_{j}\in T_{k}}\mathbbm{1}\left[e_{j}\in T\right]}{|T|}$$

Both scores $s_{i}^{phs}$ and $\hat{s}_{i}^{phs}$ are related to the Jaccard index but importantly emphasize different goals. Score $s_{i}^{phs}$ will tend to produce better scores for methods that compute shorter transcript segments; as long as the shorter transcript segments are accurate (they are contained within true transcripts), this score will be close to $1$. In contrast, $\hat{s}_{i}^{phs}$ prefers longer transcript segments and will be close to $1$ if the computed transcript is both accurate and full-length. Either $s_{i}^{phs}$, $\hat{s}_{i}^{phs}$, Jaccard index, or some linear combination thereof can be used depending on the goals of the study.

Finally, let the set of computed transcripts be $T_{(i)}^{c}=\{T_{k}\}$. An overall score for sample $i$ can then be computed as

$$\displaystyle s_{i}^{phs}=\frac{\sum_{T_{k}\in T_{(i)}^{c}}s_{i}^{phs}(T_{k})}{|T_{(i)}^{c}|}\qquad\text{and}\qquad\hat{s}_{i}^{phs}=\frac{\sum_{T_{k}\in T_{(i)}^{c}}\hat{s}_{i}^{phs}(T_{k})}{|T_{(i)}^{c}|}$$

To compute precision and recall, we first match computed transcript segment $T_{k}$ to the true transcript $T\in T^{t}$ with maximum $s_{i}^{phs}$ or $\hat{s}_{i}^{phs}$. Let the matched transcript be $T^{*}$. Then, we label each excision $e_{j}\in T_{k}$ as a true positive (TP) if $1\left[e_{j}\in T^{*}\right]=1$ and a false positive (FP) otherwise. Excisions are labeled as false negatives (FN) if they exist in a true transcript but were not included in any computed transcript. The F1 score is computed as the harmonic mean of precision and recall, which are computed as: $precision=\frac{TP}{TP+FP}$ and $recall=\frac{TP}{TP+FN}$.

The input to our model is the set of junction reads mapped to a reference genome. In this work, we map reads to the reference genome using STAR aligner (V. 2.7.3a). We used gene annotations whenever possible, including in STAR alignments since some methods require gene annotations and STAR highly recommends using them when available. Gene annotations were also used to generate the BAM files for the GEUVADIS data. However, gene annotations are not required to execute BREM.

BREM assumes that the number of SMEs ($K$) is given as input. However, the number of SMEs should not be less than the chromatic number of the excision interval graph, or, equivalently, the maximum independent set ($IS$) of the complement graph. Using the interval graph property, we can compute the chromatic number in polynomial time. Then, for each gene, we trained our model with $K=IS+x$, where $x\in\{0,2,4,6,8,10,12,14,16\}$ and selected the model with the highest predictive likelihood. Predictive likelihood is commonly used to perform model selection on admixture and topic models [47] and is less prone to overfitting than likelihood. To select hyperparameters, we implemented a grid search on held-out genes where $\alpha\in\{0.001,0.01,1,5,10\}$, $\eta\in\{0.01,1,5,10\}$, $r$ and $s\in\{1,5,10\}(r=s)$. Throughout the subsequent experiments, we set $\eta=0.01$ and $\alpha=r=s=1$ for both simulated and experimental data. For convergence, we check RFWSR every 50 iterations after burn-in (500 iterations) and stop sampling after 100 iterations if RFWSR$<\sigma$.

We applied BREM, Cufflinks, LeafCutter, StringTie, and rMATS to reconstruct transcripts, splice events, or SMEs in all $1260(420\times 3)$ simulated genes. Before computing precision, recall, and F1 score, the computed transcript segments must be matched to true transcripts; here, we quantify this using the partial homogeneity scores (Fig. 1).

Interestingly, StringTie does not suffer from the same significant decrease as Cufflinks, though both StringTie and Cufflinks exhibited high variance. In comparison, BREM demonstrated far less variability in $s^{phs}$ and $\hat{s}^{phs}$ than Cufflinks and StringTie, while maintaining high performance.

Next, to evaluate the impact of the parameter that controls the number of SMEs ($K$), we varied $K$ from the chromatic number in the excision interval graph (equivalently, size of the maximum independent set ($IS$) in the complement graph) to $IS+16$. The trend for both $s^{phs}$ (Fig. LABEL:fig:fig1, top) and $\hat{s}^{phs}$ (Fig. LABEL:fig:fig1, bottom) are similar: as $K$ increases, precision increases initially and then remains flat while recall decreases monotonically. This is likely due to two factors. First, BREM collapses SMEs with the same excision configuration after convergence. This means that the effective K is much lower when $K$ is much larger than the number of alternative transcripts. Second, BREM benefits from the flexibility of additional SMEs initially, but eventually when $K>>IS$, BREM learns SMEs that are low abundance and noisy.

Having matched computed transcripts with true transcripts, we next evaluated each method with respect to precision, recall, and F1 score for the top match using $s^{phs}$ (Fig. 2a) and $\hat{s}^{phs}$ (Fig. 2b). First, rMATS is highly selective, exhibiting high precision regardless of the length of the transcript. This, however, is to be expected since rMATS scores consistently high when considering $s^{phs}$, but also consistently low when considering $\hat{s}^{phs}$. Since rMATS is concerned only with singular splicing events, in either case, the task is less difficult. On the other hand, LeafCutter performs some local assembly of splicing events into clusters and thus has a more difficult assembly task than rMATS, though performs similarly in terms of F1 score. Both Cufflinks and StringTie exhibit high variance, but perform considerably better than the local splicing methods in terms of recall. BREM, situated between these extremes, achieves higher precision for most genes than the full-length transcript methods and substantially higher recall and F1 with lower standard errors. We also tested precision, recall, and F1 score as a function of the complexity of the overlap graph (defined by the number of edges). For genes yielding complex graphs ($|E|>200$), BREM achieves the highest recall and F1 Score, while rMATS is the most precise (Fig. LABEL:fig:prf_nf). Importantly, this shows that BREM performs well when there is substantial overlap among the transcripts. As a function of the number of excisions, BREM also achieves the highest recall (Fig LABEL:fig:prf_nodes). The flexibility of our admixture modelling allows BREM to focus on producing high confidence transcript segments rather than fixing the size to be small (e.g., individual splice events) or large (full-length transcript).

Next, we tested the sensitivity of BREM to model parameters; in particular, we tested how the mean posterior SME length (denoted $|SME|$ and defined by the number of excisions in a SME) varied as a function of $K$, $r$, and $s$. We trained the model setting $r,s\in\{0.1,1,10,100\}$ and $K\in\{IS,IS+5,IS+10,IS+15\}$. The parameter $K$ did not correlate with $|SME|$, likely due to BREM collapsing posterior SMEs with the same excision usage. However, as we increased the prior mean of $Beta(r,s)$, $|SME|$ also increased (Fig. 3). This is consistent with the interpretation of $r$ and $s$ in the model: $r$ and $s$ control the prior probability of including an excision in SMEs. As the mean of $Beta(r,s)$ increases, larger SMEs become more likely in the posterior. However, this relationship is not strictly monotonic, as other model parameters, properties of the transcripts, and stochasticity of model inference interact with the effect of $r$ and $s$ on $|SME|$. BREM is also fast, with the running time increasing linearly as a function of $K$, $|C|$, and the average number of junction reads across samples (Fig. LABEL:fig:runtime_2).

Each simulated gene consisted of $8$ groups of $100$ samples with fold changes $1$, $1$, $1$, $1.1$, $1.25$, $1.5$, $3$, and $5$. We computed differential expression for each method and all pairwise groupings of the samples ($28$ in total). Pairwise comparisons between the first three groups enabled estimation of false discoveries. We randomized the processing order of genes and allocated each method a full week on a $128$ core computer to processes the simulated data (Table 1). StringTie, rMATS, LeafCutter, and BREM all finished in less than a day, while Cufflinks only finished $21.4\%$ of configurations. Additionally, recent comparisons have shown higher ability to detect differential splicing for rMATS, StringTie, and LeafCutter when compared to Cufflinks [58, 63, 61]; thus, we excluded Cufflinks from the comparison. BREM achieved the highest sensitivity and accuracy of identifying differential usage (of SMEs), though StringTie achieved relatively high sensitivity with an impressively high specificity ($0.996$).

We applied BREM and our differential SME model to both GEUVADIS and EGA datasets. After filtering genes expressed at low levels and those without conflicts in the excision interval graph, we applied BREM to infer SMEs in $3983$ and $4278$ genes in the GEUVADIS and EGA data respectively. Using our results on the precision and recall for BREM with varying $K$ (Fig. LABEL:fig:fig1), we set $K=IS+4$. We then applied our Dirichlet Multinomial model to compute differential SME usage across the two data sets. We used the super population (African vs. European) to group samples in GEUVADIS and treatment status in the EGA dataset. After multiple comparisons correction using Benjamini-Hochberg [6], p-values were well-calibrated (Fig. LABEL:fig:exp) and we observed $2105$ and $1961$ genes with significant differential SME usage in the GEUVADIS and EGA data (FDR corrected $p<0.05$).

We conducted a gene ontology (GO) analysis using genes with differential SME expression as the target list and all genes input into BREM as the background list [12]. In the GEUVADIS data, the top $12$ GO terms in the Biological Process ontology ranked by p-value (p $<2.97\times 10^{-6}$) referenced regulation of biomolecular processes. This is consistent with a growing body of evidence that suggests splicing plays a major role in regulating gene expression [16, 17] and metabolism [23, 3, 37]. In the Molecular Function ontology, alternative splicing plays an integral role in the top $18$ GO terms, which reference ATP, DNA, drug, and other molecular binding (p $<5.69\times 10^{-5}$) [39, 19]. In the EGA data, both molecular function and biomolecular processes exhibited significant associations with regulation of and binding to kinase proteins (GO:0046330, GO:0043507, GO:0046328, GO:0019901; p $<9.7\times 10^{-4}$). Alternative splicing is known to (a) regulate the binding of kinase proteins [20] and (b) increase kinase protein diversity [2].

We quantified the total number and percentages of novel versus known splice junctions and SMEs or transcripts in both GEUVADIS and EGA datasets. Since our processing pipeline focuses on excisions, and is thus similar to LeafCutter, and we are testing reconstruction, we compared our results to only Cufflinks and StringTie. We only consider SMEs that are expressed in $10$ or more samples, where the $k^{th}$ SME is considered expressed in sample $i$ if there exists $10$ or more $z_{ij}=k$. We allowed inferred junction locations to differ by at most $6$ nucleotide bases from the reference to be considered matching. We followed the recommended pipelines for StringTie and Cufflinks and merged per-sample assemblies. We considered an SME or computed transcript as novel if it was not a subset of an annotated transcript (and known otherwise). Splice junctions are considered novel if they do not exist in the reference (and known otherwise).

All methods produce far fewer novel SMEs, transcripts, and splice junctions in GEUVADIS compared to the EGA data (Fig. LABEL:fig:novel, a and b). This may be due to the reference transcripts of lymphoblastoid cell lines being more well characterized than monocytes or due to the differences in sequencing platforms (bulk versus single cell RNA-seq). BREM and Cufflinks produced larger proportions of known to novel splice junctions and SMEs or transcripts in GEUVADIS whereas StringTie produces far more novel transcripts (Fig. LABEL:fig:novel, c and d). This discrepancy was larger in the EGA data, where StringTie produced much higher proportions of novel transcripts ($<5\%$ were observed in the reference).