ASGAL: aligning RNA-Seq data to a splicing graph to detect novel alternative splicing events

Data coming from high-throughput sequencing of RNA (RNA-Seq) can shed light on the diversity of transcripts that results from Alternative Splicing (AS). Computational approaches for transcriptome analysis from RNA-Seq data may be classified according to two primary goals: (i) detection of AS events and (ii) full-length isoform reconstruction. Tools in these two categories may be further classified based on an approach which may be (a) de-novo assembly based or (b) gene annotation guided or reference based. Various tools have been proposed in the literature that fall in the categories listed above. Examples of tools in category (ii.a) that do not require a reference genome are Trinity [1] and ABySS [2], while Cufflinks [3], Scripture [4], and Traph [5], among many others, are known tools of category (ii.b). The first two tools were originally designed for de-novo isoform prediction and can make limited use of existing annotations. While the reconstruction of full-length transcripts (either de-novo or using a reference) is a computationally intensive task, the detection of AS events is computationally feasible and it can be achieved without performing intensive steps related to transcript reconstruction. Observe that given a set of transcripts reconstructed from a sample of RNA-Seq reads, a tool for comparing transcripts is needed to extract AS events. Such a comparison is performed for example by AStalavista [6], a popular tool for the exhaustive extraction and visualization of complex AS events from full-length transcripts. This tool does not use RNA-Seq reads as input but only the gene annotation, and it does not focus on single events (such as exon skipping, alternative splice sites, etc.) but rather uses a flexible coding of AS events [7] to list all the AS events between each pair of transcripts.

Since reconstructing full-length isoforms from RNA-Seq reads is a difficult and computationally expensive problem, one may restrict the task to the direct detection of AS events from RNA-Seq data through an alignment process. Following the latter approach, we propose a computational approach to predict AS events, and we implement this procedure in a tool — ASGAL — belonging to category (i.b). Compared to existing tools, ASGAL has as main goals the splice-aware alignment of RNA-Seq data to a splicing graph and the annotation of the graph with novel splicing events that are supported by such alignments. From this perspective, differently from tools for event detection based on differential analysis, ASGAL is able to detect a novel event in a gene annotation when this event is supported by reads from a single unannotated isoform. Some tools using unannotated splice sites — hence most similar to ASGAL with respect to the goal of predicting AS events — are SpliceGrapher [8] and SplAdder [9] which take as input the spliced alignments of sequencing data (RNA-Seq data for SplAdder, and RNA-Seq data in addition to EST data for SpliceGrapher) against a reference genome, and produce an augmented graph representation of the annotated transcripts, traditionally known as the splicing graph [10], with nodes and edges that may represent novel AS events. The main task of SplAdder is the prediction of AS events that are supported by an input sample, and the quantification of those events by testing the differences between multiple samples. Two other tools whose main goal is differential alternative splicing analysis are SUPPA2 [11] and rMATS [12]. Both SUPPA2 and rMATS analyze RNA-Seq data from different samples (replicates) to obtain the set of differential alternative splicing events between the analyzed conditions. SUPPA2 is only able to detect AS events that are in the annotation, while rMATS only lists novel events that use annotated splice sites. Similarly, MAJIQ [13] analyzes RNA-Seq data and a set of (annotated) transcripts to quantify the relative abundances of a set of Local Splicing Variations which implicitly represent combinations of AS events involving both annotated and novel splice sites, but also changes of these abundances between conditions. Note that both MAJIQ and rMATS do not include an alignment step, but need an external spliced aligner such as STAR [14], while SUPPA2 requires the quantification of the input transcripts, which can be obtained by using a tool like Salmon [15]. In both cases, the identification of AS events stems from an analysis of the expression levels. A most recent tool, LeafCutter [16] analyzes RNA-Seq data and quantifies differential intron usage across samples, allowing the detection of novel introns which model complex splicing events. Like the other cited tools, LeafCutter requires as input the spliced alignments of the RNA-Seq samples of interest. Two crucial computational instruments are usually required by tools of category (i.b): an input file consisting of the alignment of RNA-Seq data to a reference genome, and a gene annotation. The first input may significantly change the performance of such tools, as the accuracy of the alignment may affect the predictions of AS events. In particular, the alignment to a reference genome is usually guided by the annotated transcripts that may be represented by a splicing graph that is then enriched with the information coming from the computed alignments.

With the main goal of enriching a gene annotation with novel AS events supported by a RNA-Seq sample, we investigated an alternative approach that directly aligns the input reads against a splicing graph representing a gene annotation. The main motivation of our proposal is that, by using the splicing graph during the alignment phase, we are able to obtain an alignment focused on enriching a gene annotation with AS events that produce novel isoforms by using annotated or unannotated splice sites with respect to the actual graph. For this purpose, we implemented ASGAL (Alternative Splicing Graph ALigner), a tool that consists of two parts: (i) a splice-aware aligner of RNA-Seq reads to a splicing graph, and (ii) a predictor of AS events supported by the RNA-Seq mappings. Currently, there are several tools for the spliced alignment of RNA-Seq reads against a reference genome or a collection of transcripts but, to the best of our knowledge, ASGAL is the first tool specifically designed for mapping RNA-Seq data directly to a splicing graph. Differently from SplAdder, which enriches a splicing graph representing the gene annotation using the splicing information contained in the input spliced alignments, and then analyzes this enriched graph to detect the AS events differentially expressed in the input samples, ASGAL directly aligns the input sample to the splicing graph of the gene of interest and then detects the AS events which are novel with respect to the input gene annotation, comparing the obtained alignments with it. More precisely, ASGAL extracts the introns supported by the alignments of reads against the splicing graph, then compares them against the input annotation to detect whether novel events may be predicted from the input reads. This allows ASGAL to detect novel event types even when the input RNA-Seq sample consists only of reads that are not consistent with the input splicing graph, because of the AS event, provided that the number of alignments confirming the AS event is above a certain threshold. Instead, SplAdder needs to receive reads originating from both transcripts that induce the AS event.

The approach of inferring AS events directly from RNA-Seq reads, without assembling isoforms, is also proposed in [17], where the main idea is to perform a de-novo prediction of some AS events from the De Brujin graph assembly of RNA-Seq data, i.e. without using any gene annotation. An investigation of the de-novo prediction of AS events directly from RNA-Seq data is also given in [18], where a characterization of the splicing graph that may be detected in absence of a gene annotation (either given as a reference or as a list of transcripts) is provided.

The ASGAL mapping algorithm improves a previous solution to the approximate pattern matching to a hypertext problem (an open problem faced in [19]). The approximate matching of a string to a graph with labeled vertices is a computational problem first introduced by Manber and Wu [20] and attacked by many researchers [21–23]. Navarro [24] improved all previous results in both time and space complexity, proposing an algorithm which requires \(\mathcal {O}{m(n+e)}\) time, where m is the length of the pattern, n is the length of the concatenation of all vertex labels, and e is the total number of edges. The method in [19] improves the latest result by Thachuk [25]: an algorithm with time complexity \(\mathcal {O}{m + \gamma ^{2}}\) using succinct data structures to solve the exact version of matching a pattern to a graph — i.e. without errors — where γ is the number of occurrences of the node texts as substrings of the pattern. The algorithm in [19] is based on the concept of Maximal Exact Match and it uses a succinct data structure to solve the approximate matching of a pattern to a hypertext in \(\mathcal {O}{m+\eta ^{2}}\) time, where η is the number of Maximal Exact Matches between the pattern and the concatenation of all vertex labels. In this paper, we improve the results in [19] by extending the algorithm to implement a RNA-Seq data aligner for detecting general AS event types from the splicing graph.

An experimental analysis on real and simulated data was performed with the purpose of assessing the quality of ASGAL in detecting AS event types that are annotated or novel with respect to a gene annotation. We note that the current implementation of ASGAL is not able to detect the insertion of novel exons inside an intron and intron retention events caused by the union of two exons. In the first part of our experimental analysis, we compared the alignment step of ASGAL with STAR, one of the best-known spliced aligner. The results show a good accuracy of ASGAL in producing correct alignments by directly mapping the RNA-Seq reads against the splicing graph of a gene. Although ASGAL works under different assumptions than other existing tools, we decided to compare ASGAL with SplAdder, rMATS, and SUPPA2. For this purpose we first ran an experimental analysis on simulated data and compared the event identification step of ASGAL. We performed two distinct analysis. In the first one, we evaluated the accuracy of the tools in predicting novel AS event types, i.e. events that are not already contained in the annotation. Instead in the second analysis, all the tools were compared to assess their accuracy in detecting AS events that are already present in the input annotation and are supported by the RNA-Seq experiments. We also ran an experimental analysis on real data with the main goal of evaluating ASGAL, SplAdder, rMATS, and SUPPA2 in identifying RT-PCR validated alternative splicing events. We performed this last experiment also to test the ability of ASGAL in detecting such events as novel ones, that is by removing the events from the input annotation and keeping their evidence only in the RNA-Seq data.

The results in the simulated scenario show that ASGAL achieved the best values of precision, recall and F-measure in predicting alternative splicing events supported by the reads that are novel compared to the annotation specified by a splicing graph. The results on real data show the ability of ASGAL to detect RT-PCR validated alternative splicing events when they are simulated as novel events with respect to the annotated splicing graph.

ASGAL (Alternative Splicing Graph ALigner) is a tool for performing a mapping of RNA-Seq data in a sample against the splicing graph of a gene with the main goal of detecting novel alternative splicing events supported by the reads of the sample with respect to the annotation of the gene. More precisely, ASGAL takes as input the annotation of a gene together with the related reference sequence, and a set of RNA-Seq reads, to output (i) the spliced alignments of each read in the sample and (ii) the alternative splicing events supported by the sample which are novel with respect to the annotation. We point out that ASGAL uses the input reference sequence only for building the splicing graph as well as for refining the alignments computed against it, with the specific goal of improving the precision in the AS event type detection. Each identified event is described by its type, i.e. exon skipping, intron retention, alternative acceptor splice site, alternative donor splice site, its genomic location, and a measure of its quantification, i.e. the number of alignments that support the identified event.

This section is organized as follows. We first introduce the basic definitions and notions that we will use in the section spliced graph-alignment, and finally we describe the steps of our method. For the sake of clarity, we will describe our method considering as input the splicing graph of a single gene: it can be easily generalized to manage more than a gene at a time. However, the current version of ASGAL tool cannot manage more than a limited set of genes. At the end of this section, we will propose a possible procedure an user can adopt to use our tool in a genome-wide analysis.

Definitions

From a computational point of view, a genome is a sequence of characters, i.e. a string, drawn from an alphabet of size 4 (A, C, G, and T). A gene is a locus of the genome, that is, a gene is a substring of the genome. Exons and introns of a gene locus will be uniquely identified by their starting and ending positions on the genome. A transcript T of gene G is a sequence 〈[a₁,b₁],[a₂,b₂],…,[a_n,b_n]〉 of exons on the genome, where a_i and b_i are respectively the start and the end positions of the i-th exon of the transcript. Observe that a₁ and b_n are the starting and ending positions of transcript T on the genome, and each [b_i+1,a_i+1−1] is an intron represented as a pair of positions on the genome. In the following, we denote by \(\mathcal {E}_{G}\) the set of all the exons of the transcripts of gene G, that is \(\mathcal {E}_{G}=\cup _{T \in \mathcal {T}}\mathcal {E}(T)\), where \(\mathcal {E}(T)\) is the set of exons of transcript T and \(\mathcal {T}\) is the set of transcripts of G, called the annotation of G. Given two exons e_i=[ a_i,b_i] and e_j=[ a_j,b_j] of \(\mathcal {E}_{G}\), we say that e_iprecedes e_j if b_i<a_j and we denote this by e_i≺e_j. Moreover, we say that e_i and e_j are consecutive if there exists a transcript \(T \in \mathcal T\) and an index k such that e_k=e_i and e_k+1=e_j, and e_i, e_j in \(\mathcal {E}(T)\).

The splicing graph of a gene G is the directed acyclic graph \(\mathcal {S}_{G} = (\mathcal {E}_{G},E)\), i.e. the vertex set is the set of the exons of G, and the edge set E is the set of pairs (v_i,v_j) such that v_i and v_j are consecutive in at least one transcript. For each vertex v, we denote by seq(v), the genomic sequence of the exon associated to v. Finally, we say that \(\mathcal {S}_{G}^{\star }\) is the graph obtained by adding to \(\mathcal {S}_{G}\) all the edges (v_i,v_j)∉E such that v_i≺v_j. We call these edges novel edges. Note that the novel edges represent putative novel junctions between two existing exons (that are not consecutive in any transcript of G). Figure 1 shows an example of the definitions of gene, exon, annotation, and splicing graph.

In the following, we will use the notion of Maximal Exact Match (MEM) to perform the spliced graph-alignment of a RNA-Seq read to \(\mathcal {S}_{G}\). Given two strings R and Z, a MEM is a triple m=(i_Z,i_R,ℓ) representing the common substring of length ℓ between the two strings that starts at position i_Z in Z, at position i_R in R, and that cannot be extended in either direction without introducing a mismatch. Computing the MEMs between a string R and a splicing graph \(\mathcal {S}_{G}\) can be done by concatenating the labels of all the vertices and placing the special symbol ϕ before each label and after the last one, obtaining a string \(Z = \phi {\texttt {seq}}(v_{1}) \phi {\texttt {seq}}(v_{2}) \phi \ldots \phi {\texttt {seq}}(v_{|\mathcal {E}_{G}|}) \phi \) that we call the linearization of the splicing graph (see Fig. 1 for an example). It is immediate to see that, given a vertex v of \(\mathcal {S}_{G}\), the label seq(v) is a particular substring of the linearization Z. For the sake of clarity, let us denote this substring, which is the one related to seq(v), as Z[ i_v,j_v]. Then, by employing the algorithm by Ohlebusch et al. [26], all the MEMs longer than a constant L between R and Z, thus between R and \(\mathcal {S}_{G}\), can be computed in linear time with respect to the length of the reads and the number of MEMs. Thanks to the special character ϕ which occurs in Z and not in R, each MEM occurs inside a single vertex label and cannot span two different labels. In the following, given a read R and the linearization Z of \(\mathcal {S}_{G}\), we say that a MEM m=(i_Z,i_R,l) belongs to vertex v if i_v≤i_Z≤j_v where [ i_v,j_v] is the interval on Z related to the vertex label seq(v) (that is, seq(v)=Z[ i_v,j_v]). We say that a MEM m=(i_Z,i_R,l) precedes another MEM \(m^{\prime }=\left (i^{\prime }_{Z}, i^{\prime }_{R}, l^{\prime }\right)\) in R if \(i_{R} < i^{\prime }_{R}\) and \(i_{R} + l < i^{\prime }_{R} + l^{\prime }\), and we denote this by m≺_Rm^′. Similarly, when m precedes m^′ in Z, we denote it by m≺_Zm^′, if the previous properties hold on Z and the two MEMs belong to the same vertex label seq(v). When m precedes m^′ in R (in Z, respectively), we say that \(lgap_{R} = i^{\prime }_{R} - (i_{R} + l)\) (\({lgap}_{Z} = i^{\prime }_{Z} - (i_{Z} + l)\), respectively) is the length of the gap between the two MEMs. If lgap_R or lgap_Z (or both) are positive, we refer to the gap strings as sgap_R and sgap_Z, while when they are negative, we say that m and m^′overlap either in R or Z (or both). Given a MEM m belonging to the vertex labeled seq(v), we denote as PREF_Z(m) and SUFF_Z(m) the prefix and the suffix of seq(v) upstream and downstream from the start and the end of m, respectively. Figure 2 summarizes the definitions of precedence between MEMs, gap, overlap, PREF_Z, and SUFF_Z.

Spliced graph-alignment

We are now able to define the fundamental concepts that will be used in our method. In particular, we first define a general notion of gap graph-alignment and then we introduce specific constraints on the use of gaps to formalize a splice-aware graph-alignment that is fundamental for the detection of alternative splicing events in ASGAL.

A gap graph-alignment of R to graph \(\mathcal {S}_{G}\) is a pair (A,π) where π=〈v₁,…,v_k〉 is a path of the graph \(\mathcal {S}_{G}^{\star }\) and

$$A = \left\langle (p_{1}, r_{1}), \left(p_{1}^{\prime}, r_{1}^{\prime}\right), \ldots, \left(p^{\prime}_{n-1}, r^{\prime}_{n-1}\right),(p_{n}, r_{n}) \right\rangle $$

is a sequence of pairs of strings, with n≥k, such that seq(v₁)=x·p₁ and seq(v_k)=p_n·y, for x,y possibly empty strings and \(P= p_{1} \cdot p_{1}^{\prime } \cdot p_{2} \cdot p_{2}^{\prime } \cdot p_{3} \cdots p_{n-1}^{\prime } \cdot p_{n} \) is the string labeling the path π and \(R = r_{1} \cdot r_{1}^{\prime } \cdot r_{2} \cdots r_{n-1}^{\prime } \cdot r_{n} \).

The pair (p_i,r_i), called a factor of the alignment A, consists of a non-empty substring r_i of R and a non-empty substring p_i of the label of a vertex in π. On the other hand, the pair \(\left (p_{i}^{\prime }, r_{i}^{\prime }\right)\) is called a gap-factor of the alignment A if at least one of \(p_{i}^{\prime }\) and \(r_{i}^{\prime }\) is an empty substring ε. Moreover, either \(p_{i}^{\prime }\) is empty or \(|p_{i}^{\prime }| > \alpha \), and either \(r_{i}^{\prime }\) is empty or \(|r_{i}^{\prime }| > \alpha \), for a fixed value α which represents the maximum alignment indel size allowed. When an insertion (or a deletion) is smaller than α, we consider it an alignment indel and we incorporate it into a factor; otherwise, we consider it as a clue of the possible presence of an AS event and we represent it as a gap-factor. We note that an “alignment indel” is a small insertion or deletion which occurs in the alignment, due to a sequencing error in the input data or a genomic insertion/deletion. Intuitively, in a gap graph-alignment, factors correspond to portions of exons covered (possibly with errors) by portions of the read, while gap-factors correspond to introns, which can be already annotated or novel, and which can be used to infer the possible presence of AS events. We note that to allow the detection of alternative splice site events known as NAGNAG resulting in a difference of 3bps, if an alignment indel occurs at the beginning or at the end of an exon, we consider it during the detection of the events, even though it is not modeled as a gap-factor since in these cases the insertion may be smaller than α.

We associate to each factor (p_i,r_i) the cost δ(p_i,r_i), and to each gap-factor \(\left (p_{i}^{\prime },r_{i}^{\prime }\right)\) the cost \(\delta \left (p_{i}^{\prime },r_{i}^{\prime }\right)\), by using a function δ(·,·) with positive values. Then the cost of the alignment (A,π) is given by the expression:

$$ {\mathtt{cost}}(A, \pi) = \sum_{i=1}^{n} \delta(p_{i}, r_{i}) + \sum_{i=1}^{n-1} \delta\left(p_{i}^{\prime},r_{i}^{\prime}\right). $$

Moreover, we define the error of a gap graph-alignment as the sum of the edit distance of each factor (but not of gap-factors). Formally, the error of the alignment (A,π) is:

$$ {\mathtt{Err}}(A, \pi) = \sum_{i=1}^{n} d(p_{i}, r_{i}), $$

where d(·,·) is the edit distance between two strings.

To define a splice-aware alignment, that we call spliced graph-alignment, we need to classify each gap-factor and to assign it a cost. Our primary goal is to compute a gap graph-alignment of the read to the splicing graph that possibly reconciles to the gene annotation; if this is not possible, then we want to minimize the number of novel events. For this reason we distinguish three types of gap-factors: annotated, novel, and uninformative. Intuitively, an annotated gap-factor models an annotated intron, a novel gap-factor represents a novel intron, while an uninformative gap-factor does not represent any intron.

Formally, we classify a gap-factor \(\left (p_{i}^{\prime }, r_{i}^{\prime }\right)\) as annotated if and only if \(p^{\prime }_{i} = r^{\prime }_{i} = \epsilon \) and the two strings p_i, p_i+1 are on two different vertices that are linked by an edge in \(\mathcal {S}_{G}\). We classify a gap-factor \(\left (p_{i}^{\prime }, r_{i}^{\prime }\right)\) as novel in the following cases:

1. 1
   \(r_{i}^{\prime } = \epsilon \) and \(p_{i}^{\prime }= \epsilon \) occurs between the strings p_i and p_i+1 which belong to two distinct vertices linked by an edge in \(\mathcal {S}_{G}^{\star }\) and not in \(\mathcal {S}_{G}\) (i.e. this gap-factor represents an exon skipping event — Fig. 3a).

   

    
2. 2

   \(r_{i}^{\prime } = \epsilon \) and \(p_{i}^{\prime }\neq \epsilon \) occurs between the strings p_i and p_i+1 which belong to the same vertex of \(\mathcal {S}_{G}^{\star }\) (i.e. this gap-factor represents an intron retention event — Fig. 3b). Actually, we note here that this type of gap-factor may represent also a genomic deletion: currently, our program does not distinguish between intron retentions and genomic deletions that are entirely contained in an exon, therefore we might overpredict intron retentions.

    
3. 3

   \(r_{i}^{\prime } = \epsilon \) and \(p_{i}^{\prime }\neq \epsilon \) occurs between the strings p_i and p_i+1 which belong to two distinct vertices linked by an edge in \(\mathcal {S}_{G}^{\star }\) (i.e. this gap-factor represents an alternative splice site event — Fig. 3c).

    
4. 4

   \(r_{i}^{\prime } \ne \epsilon \) and \(p_{i}^{\prime }= \epsilon \) occurs between the strings p_i and p_i+1 which belong to two distinct vertices linked by an edge in \(\mathcal {S}_{G}^{\star }\) (i.e. this gap-factor represents an alternative splice site extending an exon or a new exon event — Fig. 3d-e).

    

Note that Case 1 allows to detect a novel intron whose splice sites are both annotated (see Fig. 3a). Case 2 supports a genomic deletion or an intron retention (see Fig. 3b), and in case of intron retention, ASGAL finds the two novel splice sites inside the annotated exon. Case 3 gives an evidence of a novel alternative splice event shortening an annotated exon (see Fig. 3c) and ASGAL finds the novel splice site supported by this case. Finally, in Case 4, ASGAL is able to detect a novel alternative splice site (extending an annotated exon) or a novel exon (see Fig. 3d), but only in the first case (alternative splice site) ASGAL is able to find the novel splice site induced by the gap-factor.

For ease of presentation, Fig. 3 shows only “classic” AS event types and not their combination as those modeled with the notion of Local Splicing Variations (LSV) [13]. We note here that our formalization takes into account combinations of AS event types as those given by an exon skipping combined with an alternative splice site (see definition of gap-factor in cases 3 and 4). However, the actual version of the tool is designed only to detect the AS event types shown in Fig. 3. For completeness, in Fig. 4 we show the same AS event types (shown in Fig. 3) with respect to the annotated case, i.e. when the gap-factor is annotated and it represents an already known AS event.

Finally, we classify a gap-factor \(\left (p^{\prime }_{i}, r^{\prime }_{i}\right)\) as uninformative in the two remaining cases, which are (i) \(r_{i}^{\prime } = \epsilon \) and \(p_{i}^{\prime }= \epsilon \) occurs between strings p_i and p_i+1 which belong to the same vertex, and (ii) \(r_{i}^{\prime } \neq \epsilon \) and \(p_{i}^{\prime }= \epsilon \) occurs between strings p_i and p_i+1 which belong to the same vertex. We notice that in the former case, factors (p_i,r_i) and (p_i+1,r_i+1) can be joined into a unique factor.

Let \(\mathcal {G_{F}}\) be the set of novel gap-factors of a gap graph-alignment A. Then a spliced graph-alignment (A,π) of R to \(\mathcal {S}_{G}\) is a gap graph-alignment in which uninformative gap-factors are not allowed, whose cost is defined as the number of novel gap-factors, and whose error is at most β, for a given constant β which models any type of error that can occur in an alignment (sequencing errors, indels, etc). In other words, in a spliced graph-alignment (A,π), we cannot have uninformative gap-factors, and the δ function assigns a cost 1 to each novel gap-factor and a cost 0 to all other factors and annotated gap-factors: thus \({\texttt {cost}}(A,\pi) = |\mathcal {G_{F}}|\) and Err(A,π)≤β. We focus on a bi-criteria version of the computational problem of computing the optimal spliced graph-alignment (A,π) of R to a graph \(\mathcal {S}_{G}\), where first we minimize the cost, then we minimize the error. The intuition is that we want a spliced graph-alignment of a read that is consistent with the fewest novel splicing events that are not in the annotation. Moreover, among all such alignments we look for the alignment that has the smallest edit distance (which is likely due to sequencing errors and polymorphisms) in the non-empty regions that are aligned (i.e. the factors). Figure 5 shows an example of spliced graph-alignment of error value 2, and cost 2 — since it has two novel gap-factors.

In this paper we propose an algorithm that, given a read R, a splicing graph \(\mathcal {S}_{G}\), and three constants, which are L (the minimum length of a MEM), α (the maximum alignment indel size), and β (the maximum number of allowed errors), computes an optimal spliced graph-alignment — that is, among all spliced graph-alignments with minimum cost, the alignment with minimum error. The next section details how ASGAL computes the optimal spliced graph-alignments of a RNA-Seq sample to the splicing graph \(\mathcal {S}_{G}\), and how it exploits novel gap-factors to detect AS events.

ASGAL approach

We now describe the algorithm employed by ASGAL to compute the optimal spliced graph-alignments of a sample of RNA-Seq reads to the splicing graph of a gene, to be used in order to provide the alternative splicing events supported by the sample and a measure of their quantification (i.e. the number of reads supporting the event).

The ASGAL tool implements a pipeline consisting of the following steps: (1) construction of the splicing graph of the gene, (2) computation of the spliced graph-alignments of the RNA-Seq reads, (3) remapping of the alignments from the splicing graph to the genome, and (4) detection of the novel alternative splicing events. Figure 6 depicts the ASGAL pipeline.

In the first step, ASGAL builds the splicing graph \(\mathcal S_{G}\) of the input gene using the reference genome and the gene annotation, and adds the novel edges to obtain the graph \(\mathcal S^{\star }_{G}\) which will be used in the next steps.

The second step of ASGAL computes the spliced graph-alignments of each read R in the input RNA-Seq sample by combining MEMs into factors and gap-factors. For this purpose, we extend the approximate pattern matching algorithm of Beretta et al. [19] to obtain the spliced graph-alignments of the reads, which will be used in the following steps to detect novel alternative splicing events. As described before, we use the approach proposed by Ohlebusch et al. in [26] to compute, for each input read R, the set of MEMs between Z, the linearization of the splicing graph \(\mathcal {S}_{G}\), and R with minimum length L, a user-defined parameter (we note that the approach of [26] allows to specify the minimum length of MEMs). We recall that the string Z is obtained by concatenating the strings seq(v) and ϕ for each vertex v of the splicing graph (recall that ϕ is the special character used to separate the vertex labels in the linearization Z of the splicing graph). We point out that the concatenation order does not affect the resulting alignment and that the splicing graph linearization is performed only once before aligning the input reads to the splicing graph.

Once the set M of MEMs between R and Z is computed, we build a weighted graph G_M=(M, E_M) based on the parameter α, representing the maximum alignment indel size allowed, and the two precedence relations between MEMs, ≺_R and ≺_Z, respectively. Then we use such graph to extract the spliced graph-alignment. Intuitively, each node of this graph represents a perfect match between a portion of the input read and a portion of an annotated exon whereas each edge models the alignment error, the gap-factor of the spliced graph-alignment, or both. More precisely, there exists an edge from m to m^′, with m,m^′∈M, if and only if m≺_Rm^′ and one of the following six conditions (depicted in Fig. 7) holds:

1. 1
   m and m^′ are inside the same vertex label of Z, m≺_Zm^′, and either (i) lgap_R>0 and lgap_Z>0, or (ii) lgap_R=0 and 0<lgap_Z≤α. The weight of the edge (m,m^′) is set to the edit distance between sgap_R and sgap_Z (Fig. 7a).

   

    
2. 2

   m and m^′ are inside the same vertex label of Z, m≺_Zm^′, lgap_R≤0, and lgap_Z≤0. The weight of the edge (m,m^′) is set to |lgap_R−lgap_Z| (Fig. 7b).

    
3. 3

   m and m^′ are inside the same vertex label of Z, m≺_Zm^′, lgap_R≤0 and lgap_Z>α. The weight of the edge (m,m^′) is set to 0 (Fig. 7c).

    
4. 4

   m and m^′ are on two different vertex labels seq(v₁) and seq(v₂), with v₁≺v₂, and lgap_R≤0. The weight of the edge (m,m^′) is set to 0 (Fig. 7d).

    
5. 5

   m and m^′ are on two different vertex labels seq(v₁) and seq(v₂), with v₁≺v₂, lgap_R>0, and SUFF_Z(m)=PREF_Z(m^′)=ε. The weight of the edge (m,m^′) is set to 0 if lgap_R>α, and to lgap_R otherwise (Fig. 7e).

    
6. 6

   m and m^′ are on two different vertex labels seq(v₁) and seq(v₂), with v₁≺v₂, lgap_R>0, at least one between SUFF_Z(m) and PREF_Z(m^′) is not ε. The weight of the edge (m,m^′) is set to the edit distance between sgap_R and the concatenation of SUFF_Z(m) and PREF_Z(m^′) (Fig. 7f).

    

Note that the aforementioned conditions do not cover all of the possible situations that can occur between two MEMs, but they represent those that are relevant for computing the spliced graph-alignments of the considered read. Intuitively, m and m^′ contribute to the same factor (p_i,r_i) in cases 1 and 2 and the non-zero weight of the edge (m,m^′) concurs to the spliced graph-alignment error. In cases 3-5, the edge (m,m^′) models the presence of a novel gap-factor. More precisely, m contributes to the end of a factor (p_i,r_i) and m^′ contributes to the start of the consecutive factor (p_i+1,r_i+1) and the novel gap-factor in between models an intron retention or a genomic deletion on an annotated exon (case 3), an alternative splice site shortening an annotated exon (case 4), and an alternative splice site extending an annotated exon or a new exon (case 5). Finally in case 6, m contributes to the end of a factor (p_i,r_i) and m^′ contributes to the start of the consecutive factor (p_i+1,r_i+1) whereas the gap-factor in between can identify either a novel exon skipping event or an already annotated intron. In both these cases, the non-zero weight of the edge contributes to the spliced graph-alignment error.

The spliced graph-alignment of the read R is computed by a visit of the graph G_M. More precisely, each path π_M of this graph represents a spliced graph-alignment and the weight of the path is the number of differences between the pair of strings in R and Z covered by π_M. For this reason, for read R, we select the lightest path in G_M, with weight less than β (the given error threshold) which also contains the minimum number of novel gap-factors, i.e. we select an optimal spliced graph-alignment.

The third step of ASGAL computes the spliced alignments of each input read with respect to the reference genome starting from the spliced graph-alignments computed in the previous step. Exploiting the annotation of the gene, we convert the coordinates of factors and gap-factors in the spliced graph-alignment to positions on the reference genome. In fact, observe that factors map to coding regions of the genome whereas gap-factors identify the skipped regions of the reference, i.e. the introns induced by the alignment, modeling the possible presence of AS events (see Fig. 3 for details). We note here that converting the coordinates of factors and gap-factors to positions on the reference genome is pretty trivial except when factors p_i and p_i+1 are on two different vertices and only \(p^{\prime }_{i}\) is ε (case d-e of Fig. 3). In this case, the portion \(r^{\prime }_{i}\) must be aligned to the intron between the two exons whose labels contains p_i and p_i+1 as a suffix and prefix, respectively. If \(r^{\prime }_{i}\) aligns to a prefix or a suffix of this intron (taking into account possible errors within the total error bound α), then the left or right coordinate of the examined intron is modified according to the length of \(r^{\prime }_{i}\) (Fig. 3d). In the other case (Fig. 3e), the portion \(r^{\prime }_{i}\) is not aligned to the intron and it is represented as an insertion in the alignment. Moreover, the third step of our approach performs a further refinement of the splice sites of the introns in the obtained spliced alignment since it searches for the splice sites (in a maximum range of 3 bases with respect to the detected ones) determining the best intron pattern (firstly GT-AG, secondly GC-AG if GT-AG has not been found).

In the fourth step, ASGAL uses the set \(\mathcal {I}\) of introns supported by the spliced alignments computed in the previous step, i.e. the set of introns associated to each gap-factor, to detect the novel alternative splicing events supported by the given RNA-Seq sample with respect to the given annotation. Let \(\mathcal {I}_{n}\) be the subset of \(\mathcal {I}\) composed of the introns which are not present in the annotation, that is, the novel introns. For each novel intron \(\left [p_{s}, p_{e}\right ] \in \mathcal {I}_{n}\) which is supported by at least ω alignments, ASGAL identifies one of the following events, which can be considered one of the relevant events supported by the input sample:

• exon skipping, if there exists an annotated transcript containing two non-consecutive exons [ a_i,b_i] and [ a_j,b_j], such that b_i=p_s−1 and a_j=p_e+1.

• intron retention, if there exists an annotated transcript containing an exon [ a_i,b_i] such that (i) a_i<p_s<p_e<b_i, (ii) there exists an intron in \(\mathcal {I}\) ending at a_i−1 or a_i is the start of the transcript and (iii) there exists another intron in \(\mathcal {I}\) starting at b_i+1 or b_i is the end of the transcript.

• alternative acceptor site, if there exists an annotated transcript containing two consecutive exons [ a_i,p_s−1] and [ a_j,b_j] such that p_e<b_j, and there exists an intron in \(\mathcal {I}\) starting at b_j+1 or b_j is the end of the transcript.

• alternative donor site, if there exists an annotated transcript containing two consecutive exons [ a_i,b_i] and [ p_e+1,b_j] such that p_s>a_i, and there exists an intron in \(\mathcal {I}\) ending at a_i−1 or a_i is the start of the transcript.

We note here that these definitions are accurately designed to minimize the chances of mistaking a complex AS event as those modeled with the notion of LSV for an AS event. For example, if we remove conditions (ii) and (iii) from the definition of intron retention, we could confuse the situation shown in Fig. 8 with an intron retention event.