Preliminaries

Let \(\Sigma\) be an alphabet of fixed size \(\sigma\). Here we always assume \(\Sigma =\{A,C,T,G\}\). Given a sequence (string) \(s \in \Sigma ^*\), let |s| denote its length, s[i] the ith element of s, and s[i, j] the substring \(s[i] s[i+1] \ldots s[j]\) for any \(1 \le i<j \le |s|\).

A k-mer is a sequence \(s \in \Sigma ^k\). Given an integer k and a set S of sequences each of length \(n \ge k\), we define span(S, k) as the set of all distinct k-mers that appear as a substring in S.

Given a directed graph \(G = (V,A)\) and a vertex \(v \in V\), we denote its out-neighbourhood (resp. in-neighbourhood) by \(N^+(v)=\{ u \in V \mid (v,u) \in A \}\) (resp. \(N^-(v)=\{ u \in V \mid (u,v) \in A \}\)), and its out-degree (resp. in-degree by \(d^+(v)=|N^+(v)|\) (\(d^-(v)=|N^-(v)|\)). A (simple) path \(\pi = s \leadsto t\) in G is a sequence of distinct vertices \(s = v_0, \ldots , v_l = t\) such that, for each \(0 \le i < l\), \((v_i, v_{i+1})\) is an arc of G. If the graph is weighted, i.e. there is a function \(w : A \rightarrow Q_{\ge 0}\) associating a weight to every arc in the graph, then the length of a path \(\pi\) is the sum of the weights of the traversed arcs, and is denoted by \(|\pi |\).

An arc \((u,v) \in A\) is called compressible if \(d^+(u)=1\) and \(d^-(v)=1\). The intuition behind this definition comes from the fact that every path passing through u should also pass through v. It should therefore be possible to “compress” or contract this arc without losing any information. Note that the compressed de Bruijn graph [2, 3] commonly used by transcriptomic assemblers is obtained from a de Bruijn graph by replacing, for each compressible arc (u, v), the vertices u, v by a new vertex x, where \(N^-(x) = N^-(u)\), \(N^+(x) = N^+(v)\) and the label is the concatenation of the k-mer of u and the k-mer of v without the overlapping part (see Fig. 1).

Simple uniform model for repeats

We now present the model we adopted for representing high copy-number repeats, e.g. transposable elements, in a genome or transcriptome. First, we would like to clarify that our model is a simple one and, as such, should be seen as only a first approximation, yet realistic enough, of what may happen in reality. We consider here that sequencing errors can be successfully removed. Indeed, there are several techniques to remove the big majority of the sequencing errors in RNA-seq data. In KisSplice, for example, we prune the de Bruijn graph using an absolute and a relative cut-off based on the k-mer coverage. The absolute cut-off enables us to remove sequencing errors in general, and the relative one is tailored to deal with highly-expressed genes (more details can be found in [14]). Furthermore, while we realise that there is room for improvement, in practice, the sequencing-error-removal procedure in KisSplice seems to be effective, as most sequencing errors are removed at the expense of losing some rare genomic variants [14].

Basically, our model consists of several “similar” sequences, each generated by uniformly mutating a fixed initial sequence. In particular, it enables to model well recent invasions of transposable elements which often involve high copy-number and low divergence rate (i.e. divergence from their consensus sequence). Consider indeed as an example the recent subfamilies AluYa5 and AluYb8 with 2640 and 1852 copies respectively, which both present a divergence rate below \(1\%\) [15] (see [16] for other subfamilies with high copy-number and low divergence).

The model is as follows. First, due to mutations, the sequences \(s_1, \ldots , s_m\) that represent the repeats are not identical. However, provided that the number of such mutations is not high (otherwise the concept of repeats would not apply), the repeats are considered “similar” in the sense of having a small pairwise Hamming distance between them. We recall that, given two equal length sequences \(s\) and \(s'\) in \(\Sigma ^n\), their Hamming distance, denoted by \(d_H(s,s')\), is the number of positions i for which \(s[i] \ne s'[i]\). Indels are thus not considered in this model.

The model has then the following parameters: \(\Sigma\), the length n of the repeat, the number m of copies of the repeat, an integer k (for the length of the k-mers considered), and the mutation rate, \(\alpha\), i.e. the probability that a mutation happens in a particular position. The sequences \(s_1, \ldots , s_m\) are then generated by the following process. We first choose uniformly at random a sequence \(s_0 \in \Sigma ^n\). At step \(i \le m\), we create a sequence \(s_i\) as follows: for each position j, \(s_i[j]=s_0[j]\) with probability \(1-\alpha\), whereas with probability \(\alpha\) a value different from s[j] is chosen uniformly at random for \(s_i[j]\). We repeat the whole process m times and thus create a set \(S(m,n,\alpha )\) of m such sequences from \(s_0\) (see Fig. 2 for a small example). The generated sequences thus have an expected Hamming distance of \(\alpha n\) from \(s_0\).

Topological characterisation of the subgraphs generated by repeats

Given a de Bruijn graph \(G_k(R)\), if a is a compressible arc labelled by the sequence \(s=s_1 \ldots s_{k+1}\) then, by definition, a is the only outgoing arc of the vertex labelled by the sequence s[1, k] and the only incoming arc of the vertex labelled by the sequence \(s[2,k+1]\). Hence the \((k-1)\)-mer s[2, k] appears as a substring in R, always preceded by the symbol s[1] and followed by the symbol \(s[k+1]\). We refer to such \((k-1)\)-mers as being boundary rigid. It is not difficult to see that the set of compressible arcs in a de Bruijn graph \(G_k(R)\) stands in a one-to-one correspondence with the set of boundary rigid \((k-1)\)-mers in R.

We now calculate and compare among them the expected number of compressible arcs in \(G=G_k(R)\) when R corresponds to a set of sequences that are generated: (1) uniformly at random, and (2) according to our model. We show that \(\gamma\) is “small” in the cases where the induced graph corresponds to similar sequences, which provides evidence for the relevance of this parameter.

We consider now \(\gamma (G_k(R))\) for \(R=S(m,n,\alpha )\). We upper bound the expected number of compressible arcs by upper bounding the number of boundary rigid \((k-1)\)-mers.

The previous result shows that the number of compressible arcs is a good parameter for characterising a repeat-associated subgraph.

Identifying a repeat-associated subgraph

As we showed, a subgraph due to repeated elements has a distinctive feature: it contains few compressible arcs. Based on this, a natural formulation to the repeat identification problem in RNA-seq data is to search for large enough subgraphs that do not contain many compressible arcs. This is formally stated in Problem 1. In order to disregard trivial solutions, it is necessary to require a large enough connected subgraph, otherwise any set of disconnected vertices or any small subgraph would be a solution. Unfortunately, we show that this problem is NP-complete, so an efficient algorithm for the repeat identification problem based on this formulation is unlikely.

In Theorem 2, we prove that this problem is NP-complete for all directed graphs with (total) degree, i.e. sum of in and out-degree bounded by 3. The reduction is from the Steiner tree problem which requires finding a minimum weight subgraph spanning a given subset of vertices. It remains NP-hard even when all arc weights are 1 or 2 (see [17]). This version of the problem is denoted by STEINER(1, 2). More formally, given a complete undirected graph \(G = (V,E)\) with arc weights in \(\{1,2\}\), a set of terminal vertices \(N \subseteq V\) and an integer B, it is NP-complete to decide if there exists a subgraph of G spanning N with weight at most B, i.e. a connected subgraph of G containing all vertices of N.

We specify next a family of directed graphs that we use in the reduction. Given an integer x, we define the directed graph R(x) as a cycle on 2x vertices numbered in a clockwise order and where the arcs have alternating directions, i.e. for any \(i \le x\), \((v_{2i},v_{2i+1})\) is an arc. Observe that in R(x), all vertices in even positions, i.e. all vertices \(v_{2i}\), have out-degree 2 and in-degree 0, while all vertices \(v_{2i+1}\) have out-degree 0 and in-degree 2. Clearly, none of the arcs of R(x) is compressible.

We can obtain the same result for the specific case of subgraphs of de Bruijn graphs. The reduction is more technical but follows similarly.

Enumerating bubbles avoiding repeats

In this section, we modify the algorithm of [1] to enumerate all bubbles with at most b branching vertices in each path. Given a weighted directed graph \(G = (V,E)\) and a vertex \(s \in V\), let \(\mathcal {B}_s(G)\) denote the set of \((s,*,\alpha _1,\alpha _2,b)\)-bubbles of G. The algorithm recursively partitions the solution space \(\mathcal {B}_s(G)\) at every call until the considered subspace is a singleton (contains only one solution), and in that case it outputs the corresponding solution. In order to avoid unnecessary recursive calls, it maintains the invariant that the current partition contains at least one solution. The algorithm proceeds as follows.

In order to maintain the invariant (INV), we only perform the recursive calls when \(\mathcal {B}_{s}(\pi _1 \cdot e, \pi _2, G' - u)\) or \(\mathcal {B}_{s}(\pi _1,\pi _2, G'')\) are non-empty. In both cases, we have to decide if there exists a pair of (internally) vertex-disjoint paths \(\bar{\pi }_1 = u_1 \leadsto t_1\) and \(\bar{\pi }_2 = u_2 \leadsto t_2\), such that \(|\bar{\pi }_1| \le \alpha _1'\), \(|\bar{\pi }_2| \le \alpha _2'\), and \(\bar{\pi }_1,\bar{\pi }_2\) have at most \(b_1,b_2\) branching vertices, respectively. Since both the length and the number of branching vertices are monotonic properties, i.e. both are smaller for a prefix instead of for the full path, we can drop the vertex-disjoint condition. Indeed, let \(\bar{\pi }_1\) and \(\bar{\pi }_2\) be a pair of paths satisfying all conditions but the vertex-disjoint one. The prefixes \(\bar{\pi }^*_1 = u_1 \leadsto t^*\) and \(\bar{\pi }^*_2 = u_2 \leadsto t^*\), where \(t^*\) is the first intersection of the paths, satisfy all conditions and are internally vertex-disjoint.

Moreover, using a dynamic programming algorithm, we can obtain the following result.

As a corollary of Lemma 1, we can decide if \(\mathcal {B}_{s}(\pi _1,\pi _2, G)\) is non-empty in O(b|E|) time. Now, using an argument similar to [1], i.e. the leaves of the recursion tree and the solutions are in one-to-one correspondence and the height of the recursion tree is bounded by 4b, we obtain the following theorem.

Local assembly: experimental setup

To evaluate the performance of our method, we simulated RNA-seq data using the FluxSimulator version 1.2.1 [20]. We generated 100 million reads of 75 bp using its default error model. We used the RefSeq annotated Human transcriptome (hg19 coordinates) as a reference and we performed a two-step pipeline to obtain a mixture of mRNA and pre-mRNA (i.e. with introns not yet spliced). To achieve this, we first ran the FluxSimulator with the Refseq annotations. We then modified the annotations to include the introns and re-ran it on this modified version. In this second run, we additionally constrained the expression values of the pre-mRNAs to be correlated to the expression values of their corresponding mRNAs, as simulated in the first run. Finally, we mixed the two sets of reads to obtain a total of 100M reads. We tested two values, namely 5 and 15% for the proportion of reads from pre-mRNAs. Those values were chosen so as to correspond to realistic ones (see “Measuring the confidence of a transcript in full-length transcriptome assemblers” section).

On these simulated datasets, we ran KisSplice [12] versions 2.1.0 (Ks_2.1.0) and 2.2.0 (Ks_2.2.0, with a maximum number of branching vertices set to 5) and obtained lists of detected bubbles that are putative alternative splicing (AS) events. We also ran the full-length transcriptome assemblers Trinity version r2013_08_14 and Oases version 0.2.09 on both datasets, obtaining a list of predicted transcripts, from which we then extracted a list of putative AS events. For Oases, there was one locus in each dataset for which we were not able to extract the putative AS events. A manual inspection revealed that they were mostly composed of subparts of introns or UTRs flanked by repeats, usually copies of ALUs. The presence of such high copy-number repeats in these transcripts induced the clustering of all these unrelated sequences into one complex locus. More precisely, in the dataset containing 5% of the reads from pre-mRNAs, the largest locus that Oases predicted comprised 20,769 transcripts, while the second largest contained only 139 transcripts. In the other simulated dataset, the largest locus comprised 39,389 transcripts, and the second largest contained only 205 transcripts. This indicates that Oases faces similar issues to Ks_2.1.0. For fairness of comparison, we did not post-process these complex loci, and we were then unable to retrieve the potential AS events that they could describe. It is worth mentioning though, that the majority of the transcripts inside these loci corresponded to subparts of introns or UTRs, and not to full introns or exons, and therefore could not describe AS events.

In order to assess the precision and the sensitivity of our method, we compared our set of found AS events to the set of true AS events. Following the definition of Astalavista, an AS event is composed of two sets of transcripts, the inclusion/exclusion isoforms respectively. We consider that an AS event is true if at least one transcript among the inclusion isoforms and one among the exclusion isoforms is present in the simulated dataset with at least 5 reads per kilobase (RPK). The rationale for adding this threshold is that very rare events are considerably hard, or even impossible, to detect by all methods.

To compare the results of KisSplice with the true AS events, we propose that a true AS event is a true positive (TP) if there is a bubble such that one path matches the inclusion isoform and the other the exclusion isoform. If there is no such bubble among the results of KisSplice, the event is counted as a false negative (FN). If a bubble does not correspond to any true AS event, it is counted as a false positive (FP). To align the paths of the bubbles to transcript sequences, we used the Blat aligner [21] with 95% identity and a constraint of 95% of each bubble path length to be aligned (to account for the sequencing errors simulated by FluxSimulator). We computed the sensitivity TP/(TP+FN) and precision TP/(TP+FP) for each simulation case and we report their values for various classes of expression of the minor isoform. Expression values are measured in RPK.

Local assembly: results

The overall sensitivity and precision of Ks_2.2.0, Ks_2.1.0, Trinity and Oases can be found in Fig. 7a, b, respectively. A first look reveals that Ks_2.2.0 outperforms the other three methods in both measures and datasets.

A closer look at Fig. 7a shows that both versions of KisSplice had better sensitivity than both transcriptome assemblers in the 5% pre-mRNA dataset. However, due to its inability to deal with repeat-associated regions, the performance of Ks_2.1.0 drops substantially, from 46 to 33%, when a higher rate of 15% of pre-mRNA is present in the data. The same happened with Oases. Ks_2.2.0 and Trinity, on the other hand, were able to slightly improve their sensitivity from the 5 to the 15% pre-mRNA dataset. It is however worth mentioning that the sensitivity of Ks_2.2.0 is substantially higher than the one of Trinity in the 15% pre-mRNA dataset. In summary, we can say that Ks_2.2.0 is substantially more sensitive than all the other three methods. This reflects the fact that most problematic repeats are in introns. A small unspliced mRNA rate leads to few repeat-associated subgraphs, so there are not many AS events drowned in them. In this case, the advantage of using Ks_2.2.0 is less obvious, whereas a large proportion of pre-mRNA leads to more AS events drowned in repeat-associated subgraphs which are identified by Ks_2.2.0 and usually missed by the other methods.

In Fig. 7b we can see that Ks_2.2.0 and Trinity presented the highest precision (98%) of all methods in the 5% pre-mRNA dataset. Although Ks_2.1.0 is ranked third, it still presents a very high precision (95%), while Oases presented a moderate value (80%). Nevertheless, the most important aspect to be observed in Fig. 7b is that Ks_2.2.0 kept the same high precision even when a higher rate of 15% of pre-mRNA is present in the data. Trinity, on the other hand, dropped significantly from 98 to 79%. This drop in precision is actually mostly due to the prediction of a large number of intron retentions, since Trinity assembles both the mRNA and the pre-mRNA. Ks_2.1.0 dropped slightly from 95 to 91%, and Oases dropped moderately, from 80 to 70%. In summary, we can say that both versions of KisSplice are more precise than both transcriptome assemblers, except that Trinity shows comparable precision if a small rate of pre-mRNA is present in the data, and, more specifically, that Ks_2.2.0 outperformed all the other three methods. The high precision we obtain indicates that we have very few false positives. Those mostly correspond to repeat-induced bubbles mistakenly identified as AS events.

Finally, Fig. 7c, d present the detailed sensitivity by expression levels of the four methods on both datasets, allowing for a better understanding of their performance. As we can see, Ks_2.2.0 presented the best sensitivity in practically all expression levels in both datasets, while the other three methods were worse, but comparable between themselves. We can also observe that the gap between the sensitivity of Ks_2.2.0 and the sensitivity of the other methods tends to increase as the expression levels of the genes increase, especially in the 15% pre-mRNA dataset. Since highly-expressed genes tend to present higher levels of pre-mRNA content, they bring several repeat copies in their introns, and thus some parts of their associated graphs are complex and repeat-induced. Therefore, AS events inside such genes tend to be trapped in troublesome regions of the graph, making them harder to find. As Ks_2.2.0 is able to avoid such complex regions and retrieve the AS events deeply plunged into them, it presents better sensitivity than the other methods, especially in highly-expressed genes and datasets with high rate of pre-mRNAs.

As was already reported in [12], KisSplice (i.e. both Ks_2.2.0 and Ks_2.1.0) is faster and uses considerably less memory than Trinity and Oases. For instance, on these datasets, KisSplice uses around 5 GB of RAM, while Trinity uses more than 20 GB, and Oases, around 18 GB. However, it should be noted that both these latter methods try to solve a more general problem than KisSplice, that is reconstructing the full-length transcripts.

To conclude, our results show that Ks_2.2.0 is significantly more sensitive and precise than Ks_2.1.0, Trinity and Oases for the task of identifying AS events. The advantage of using Ks_2.2.0 over the other three methods is more evident when the input data contains high pre-mRNA content or the AS events of interest stem from highly-expressed genes.

On the usefulness of Ks_2.2.0 on real data

In order to give an indication of the usefulness of our repeat-avoiding bubble enumeration algorithm with real data, we also ran Ks_2.2.0 and Ks_2.1.0 on the SK-N-SH Human neuroblastoma cell line RNA-seq dataset (wgEncodeEH000169, total RNA). In Fig. 8, we have an example of a non-annotated exon skipping event not found by Ks_2.1.0. Observe that the intronic region contains several transposable elements (many of which are Alu sequences), while the exons contain none. This is a good example of a bubble (exon skipping event) drowned in a complex region of the de Bruijn graph. The bubble (composed by the two alternative paths) itself contains no repeated elements, but is surrounded by them. In other words, this is a bubble with few branching vertices that is surrounded by repeat-associated subgraphs. Since Ks_2.1.0 is unable to differentiate between repeat-associated subgraphs and the bubble, it spends a prohibitive amount of time in the repeat-associated subgraph and fails to find the bubble.

Global assembly

To test our hypothesis that the Branching Measure is able to identify problematic transcripts, we evaluated it on the transcripts assembled by Trinity on the two simulated RNA-seq datasets described in “Local assembly: results” section, and on two other real RNA-seq datasets: one from the GEUVADIS project [18]¹ and one from a neuroblastoma SK-N-SH cell line (ENCODE) differentiated or not using retinoic acid.² Even though our method is reference-free, we have chosen to evaluate it under a model species so that we could make use of annotated reference genomes to assess if our predictions are correct. We compared our measure against two state-of-the-art methods for de novo transcriptome evaluation, Rsem-Eval [4] and TransRate [5], on the specific task of identifying chimeric transcripts in Trinity’s assemblies on all four described datasets. In all our tests, we used the contig impact score of Rsem-Eval as a measure of contig correctness. Formally, the contig impact score is a statistical measure that compares the hypothesis that a particular contig (i.e. transcript) is a true contig with the null hypothesis that the reads composing the contig actually represent the background noise [4]. In other words, the contig impact score determines the relative contribution of each transcript to explaining the assembly. Well-assembled transcripts should therefore have a high contig impact score, and badly assembled transcripts, including chimeras, should have a low contig impact score. TransRate [5], on the other hand, allowed us to work with a specific metric representing the probability that a contig is derived from a single transcript. This metric denotes the probability that the read coverage of a transcript is best modelled by a single Dirichlet distribution, rather than two or more distributions, and it corresponds to the field sCseq of TransRate’s output file contigs.csv.

As was shown before, one of the main errors that transcriptome assemblers do is to build chimeric transcripts. We compared the performances of the Branching Measure, Rsem-Eval, and TransRate on identifying chimeric transcripts. In order to have our ground truth, we first identified which assembled transcripts are chimeric with respect to a reference genome by using Algorithm 1. In total, 253 out of 18,706 transcripts (1.3%) in the 5% pre-mRNA dataset, 376 out of 26,407 transcripts (1.4%) in the 15% pre-mRNA dataset, 375 out of 99,591 transcripts (0.3%) in the GEUVADIS dataset, and 2830 out of 457,383 transcripts (0.6%) in the SKNSH dataset were classified as chimeric. Figure 9 depicts four ROC curves showing the performance of the three methods on all datasets. We can observe that the Branching Measure outperforms both Rsem-Eval and TransRate by a large margin in all tests and, with a high-value threshold, is also able to flag a majority of the chimeric transcripts while keeping a low false positive rate. These experiments show that, in the provided datasets, chimeric transcripts could be well captured by the Branching Measure. Our false positives include well-assembled transcripts traversing high copy-number low divergence repeats, and our false negatives include chimeric transcripts that did not go through a complex region. The main issue with Rsem-Eval and TransRate, on the other hand, is that both methods failed to find chimeric transcripts assembled from genes with similar expression levels. These chimeras had low scores and corresponded to the false negatives at the end of the ROC curves for Rsem-Eval and TransRate. As a side effect of this misclassification, many well-assembled transcripts had higher scores than several real chimeras, and were mistakenly flagged as chimeras.