Pegasus: a comprehensive annotation and prediction tool for detection of driver gene fusions in cancer

Abstract

Background

The extraordinary success of imatinib in the treatment of BCR-ABL1 associated cancers underscores the need to identify novel functional gene fusions in cancer. RNA sequencing offers a genome-wide view of expressed transcripts, uncovering biologically functional gene fusions. Although several bioinformatics tools are already available for the detection of putative fusion transcripts, candidate event lists are plagued with non-functional read-through events, reverse transcriptase template switching events, incorrect mapping, and other systematic errors. Such lists lack any indication of oncogenic relevance, and they are too large for exhaustive experimental validation.

Results

We have designed and implemented a pipeline, Pegasus, for the annotation and prediction of biologically functional gene fusion candidates. Pegasus provides a common interface for various gene fusion detection tools, reconstruction of novel fusion proteins, reading-frame-aware annotation of preserved/lost functional domains, and data-driven classification of oncogenic potential. Pegasus dramatically streamlines the search for oncogenic gene fusions, bridging the gap between raw RNA-Seq data and a final, tractable list of candidates for experimental validation.

Conclusion

We show the effectiveness of Pegasus in predicting new driver fusions in 176 RNA-Seq samples of glioblastoma multiforme (GBM) and 23 cases of anaplastic large cell lymphoma (ALCL). Contact: [email protected].

Background

Gene fusions are the result of genetic aberrations (translocations, deletions, amplifications and inversions) involving the juxtaposition of two genes that can generate a single hybrid transcript. Since 1960, gene fusions have been known to play a major role in tumorgenesis. The BCR-ABL1 gene fusion, arising from the Philadelphia chromosome (t(9;22)(q34;q11)), was the first case of a translocation-induced gene fusion associated with the development of a cancer, namely chronic myelogenous leukemia [1]. In this fusion, the N-terminus oligomerization domain of BCR and the tyrosine kinase domain in ABL1 are essential in promoting oncogenic activity [2]. Among the gene fusions associated with tumor development, it is worth mentioning TMPRSS2-ERG, a gene fusion occurring in 40-80% of cases of prostate cancer [3], [4] and fusions involving the ALK gene with different partners in various malignancies [5], such as NPM1-ALK in anaplastic large cell lymphoma (ALCL) [6] and ELM4-ALK in non-small-cell lung cancer [7].

Discovering the relationship between gene fusions and cancer is gaining significant momentum thanks to advances in next generation sequencing (NGS) technology, particularly RNA paired-end sequencing [8]. Recently, the application of this technology allowed the discovery of new chromosomal rearrangements of the CIITA gene with various promiscuous partners in the lymphomagenesis of primary mediastinal B cell lymphomas [9]. In Singh et al. [10], the analysis of RNA-Seq data led to the discovery of the highly oncogenic fusion protein FGFR3-TACC3 in 3% of patients diagnosed with glioblastoma multiforme (GBM). Even though FGFR3-TACC3 occurs at low frequency, the efficacy of FGFR inhibitors in the treatment of these tumors opens the door to personalized therapies for this deadly disease. Moreover, the FGFR3-TACC3 fusion has been found in other cancers such as bladder [11] and lung [12]. These recent discoveries underscore the power of high throughput genomics for the identification of targetable gene fusions, opening the door to personalized cancer therapies.

Several bioinformatics tools are now established for the detection of candidate fusion events from paired-end RNA-Seq data. Generally, the detection of read pairs that discordantly map to two distinct genes generates a first set of gene fusion candidates. Subsequently, the exact fusion junction is determined for each candidate by searching for reads spanning the breakpoint, i.e. reads that partially map to both genes. FusionSeq [13] and deFuse [14] were the earliest examples of software based on this strategy. Detection tools differ in the type and number of cascading filters they apply to reduce the large number of false positive fusions. ChimeraScan [15] implements an algorithm based on trimming reads to increase fusion detection sensitivity. Bellerophontes [16] uses TopHat [17] and Cufflinks [18] to identify gene fusions involving truly expressed genes, and applies a set of modular cascading filters based on an accurate gene fusion model [19]. A comprehensive comparison of fusion detection tools has recently been published [20].

The methods adopted by fusion detection tools to shrink the list of candidates lead to increased specificity but reduced sensitivity. As reported in the comparative analysis performed by Abateet al.[16] , the heterogeneity of filtering strategies often yields poorly overlapping sets of candidate transcripts between algorithms. The union of all candidate fusions reported by different detection tools should be considered for further experimental validation, in order to maximize sensitivity. A problem arises however, since the number of putative gene fusions might be on the order of hundreds of candidates per RNA-Seq sample. This is largely due to the presence of read-through events, reverse transcriptase template switching artifacts, and different systematic errors in the analysis of the reads [21]. The naïve approach of considering all candidates from all detection tools quickly overwhelms the capacity of experimental validation procedures, and highlights the need to focus on a reduced number of select biologically relevant fusions driving the oncogenic progression of disease.

The classification of gene fusions into driver and passenger events is a complex problem that has not been fully explored yet. To address this issue, several databases have collected hundreds of chromosomal translocations involved in cancer cases and reported in the biomedical literature. For instance, Mitelman [22], TICdb [23] and ChimerDB2.0 [24] are manually curated repositories of known gene fusions along with detailed information such as chromosomal breakpoints, reported tissue types, and fusion sequences. New computational approaches to nominate biologically relevant fusions from high-throughput data have been proposed. ConSig assesses driver gene fusions by combining copy number variations (CNV), ontologies and interactomes based on the assumption that fusion events are more likely to arise from genes with similar biological functions [25]. Wu et al. have proposed a network based approach relying on relative co-occurrence of protein domains and domain-domain interactions, and location of the gene fusion in a gene network [26]. Recently, Oncofuse has improved the computational analysis with a machine learning approach based on a Naïve Bayes classifier applied to preserved domains after chromosomal rearrangement [27]. Compared to earlier methods, Oncofuse introduces a new level of detail by considering only the domains that are maintained on the resulting fusion transcripts. The domain analysis should be extended, however, by taking into account all possible transcript isoforms as well as the reading frame, which plays a crucial role since frame-shifted fusions imply a loss of the 3’-gene domains. Moreover, Oncofuse relies on a Naïve Bayes classifier that makes a restrictive assumption on the class conditional independence of all features. Taking the FGFR3-TACC3 gene fusion as an example, however, the acquired coiled-coil domain of the TACC3 gene cooperates with tyrosine kinase functionality of FGFR3 to produce the dramatic oncogenic effect [10]. This example illustrates the limitations of a model assumption that ignores interactions between functional protein domains.

In this paper we aim to discern oncogenic driver fusions from the background of passenger events and artifacts by combining 1) functional domain annotation based on accurate fusion sequence analysis and 2) a binary classification algorithm using gradient tree boosting. The implementation of this methodology is Pegasus, a new framework for the functional characterization of RNA-Seq gene fusion candidates and quantification of their oncogenic potential. Pegasus runs on top of multiple state of the art fusion detection tools in order to maximize detection sensitivity and consider the largest possible set of fusion candidates.

The main innovative steps introduced by Pegasus are as follows:

• Common interface between several fusion detection tools.

• Chimeric transcript sequence reconstruction: a key feature since fusion detection tools do not report whole transcript sequences.

• Reading frame identification and accurate domain annotation, including both preserved and lost protein domains within the assembled chimeric transcript.

• Prediction of fusion oncogenic potential: high performance ensemble learning technique trained on a feature space of protein domain annotations.

• Automated workflow that would otherwise require massive effort if manually executed by the scientist.

We assess the trained Pegasus model’s prediction accuracy by applying it to a set of recently discovered gene fusions where it compares quite favorably with the current state of the art, Oncofuse. Beyond curated datasets, we report the results of Pegasus on real RNA-Seq data from three distinct patient cohorts: public GBM samples from TCGA, non-public GBM samples, and non-public ALCL samples. We successfully identify driver gene fusions in both cancer types and demonstrate the utility of coupling our algorithm with experimental analysis.

Implementation

In order to first motivate our feature engineering, we briefly review the main mechanisms hitherto identified in oncogenic gene fusions (see Figure 1). Fusion transcripts can broadly lead to three scenarios: i) enhanced overexpression of an oncogene ii) deregulation of a tumor suppressor gene iii) formation of a new, aberrant protein.

Figure 1

Common gene fusions in cancer biology. Different mechanisms of gene fusions in cancer are shown. Figure 1 (a) Shows the mechanisms of oncogenic transcription factor activation by means of an endogenous enhancer. Figure 1 (b) Shows an example of disrupting gene fusion where the 5’ gene leads to the deregulation of the 3’ tumor suppressor gene. Finally, Figure 1 (c) and 1 (d) depicts the BCR-ABL1, NPM1-ALK and FGFR3-TACC3 chimeras where a completely new protein is produced.

Enhanced overexpression of an oncogene is exemplified by the famous IgH-MYC fusion (Figure 1a), and is the main reason for our explicit annotation of oncogene status and interactions with known oncogenes in our feature space representation of fusion transcripts. In other cases, deregulating properties can be associated with the fused transcript, such as insertion of one or two nucleotides across the junction breakpoint introducing a shift of the reading frame. This scenario is illustrated in the PPP2RA-CHEK2 fusion [28] (Figure 1b) where the introduced frame-shifted sequence prevents the formation of the CHEK2 protein that is a known tumor suppressor gene. Here we see the motivation for our explicit annotation of tumor suppressor status and interactions with known tumor suppressors in the feature space, as well as the need for computing reading frame of each candidate fusion. Finally, fusion transcripts can also yield a completely new chimeric protein. BCR-ABL1 [1] and NPM1-ALK [6] are well studied examples of such in-frame fusions. The new protein is generally larger than the kinase involved and causes an increase of the tyrosine kinase activity (Figure 1c). Moreover, in the recently discovered FGFR3-TACC3 gene fusion, the acquired coiled-coil domain of TACC3 gene drives the localization of the fusion protein to the mitotic spindle through a mechanism that is dependent on tyrosine kinase functionality [10] (Figure 1d). It seems that a reasonable feature space representation for predicting the oncogenic properties of such novel chimeric proteins should maintain knowledge of both preserved and lost functional domains in the partner genes.

The methodology of Pegasus is composed of three phases (see Figure 2): a) integration of candidates from fusion detection tools b) chimeric transcript sequence reconstruction and domain annotation c) classifier training and driver prediction. The first phase, Fusion Detection Tools Integration, involves pooling the entire set of unique gene fusion candidates detected by any of the fusion detection tools. The second phase, Chimeric Transcript Sequence Reconstruction and Domain Annotation, includes two steps: i) reconstruction of the chimeric transcript using the genomic breakpoint coordinates and the partner gene annotations ii) annotation of the assembled sequence to provide information on the fusion frame and to generate a report of all the protein domains conserved or lost in the gene fusion. The final phase frames Driver Fusion Prediction as a binary classification task and fits an ensemble of decision tress via the gradient boosting algorithm.

Figure 2

Pegasus fusion annotation flow. For each phase, the figure shows how the feature vector is constructed on the left side of the panel. In Fusion Detection Tools Candidates Integration step, report files from several fusion detection tools are loaded in a unique fusion database. In Chimeric Transcript Sequence Reconstruction and Functional Analysis phase, the fusion transcript is assembled according to the fusion breakpoint coordinates, the reading frame is checked and the protein domain annotation is performed on the resulting fused sequence. Finally, the Driver Fusion Prediction applies machine learning techniques to determine prediction scores.

Fusion detection tools integration

The Fusion Detection Tools Integration is the repository of the entire set of fusion candidates detected by any of the fusion detection tools. Several fusion detection algorithms are supported in Pegasus: Bellerophontes, deFuse, and ChimeraScan. Each tool adopts a private formalism for reporting fusion information with different levels of detail. However, some chimeric fusion features are common to all the fusion detection tool reports (e.g. genes involved in the fusion, genomic breakpoint coordinates, number of reads encompassing and spanning the fusion breakpoint, etc.). Thus, the internal database structure of Pegasus provides a unique point of access for all the information needed to fully describe a gene fusion candidate. Furthermore, experimental analysis might involve the comparison of several RNA-Seq samples per case study. To this end, the common repository embedded in Pegasus provides an organized overview of all the fusions occurring in the entire sample set. This feature allows comparison and the recurrence analysis of the fusion candidates within both the same experimental dataset (samples of the same disease) and within different experimental datasets (samples across different diseases).

Chimeric transcript sequence reconstruction and domain annotation

For each gene fusion candidate, the entire chimeric transcript sequence is first assembled according to publicly available gene annotations and the fused gene breakpoint coordinates. This is the most computationally intensive step in the methodology. For each gene fusion candidate, Pegasus assembles the chimeric sequence based on the possible isoforms and splicing junctions of each gene, as well as the genomic breakpoint coordinates (Figure 3). It is worth specifying that Pegasus reconstructs the fusion sequence exclusively on the basis of gene annotation and fusion breakpoint, and it does not exploit the sequenced reads because they are not an input to the program. Therefore the reconstructed sequence may not reflect the actual sequences especially in case of alternative splicing events.

Figure 3

Chimeric transcript sequence reconstruction. Both for the 5’ and the 3’ gene the annotated isoform sequences are retrieved. All the possible combinations between the isoforms of the 5’ and 3’ genes are considered as putative fusion transcripts.

We adopt the annotation file from ENSEMBL database [29]. Since several distinct isoforms might be available for a specific gene, Pegasus considers the combinations of all possible isoforms reported in the annotations of those genes involved in the fusion (Figure 3). The chimeric transcript sequence is therefore reconstructed combining the 5’ gene isoform sequence (from the isoform start codon to the genomic breakpoint) and the 3’ gene isoform sequence (from the genomic breakpoint to the isoform stop codon). Different gene isoforms allow for different protein domains to be retained or disrupted during the fusion. If this scenario occurs, Pegasus considers the union of all possible domains that are retained and lost and as input features for downstream classification. Furthermore, the fusion breakpoint can fall in either the coding region (exon-exon junction boundaries), or in non-coding regions (exon-intron or intron-intron junction boundaries). Pegasus takes the latter scenario into account and if the fusion breakpoint falls in an intron, the intronic sequence is retained.

After sequence reconstruction we assess the preservation of the reading frame in the chimeric transcript, which enables a great deal of our downstream feature engineering (see Figure 4a). If the gene fusion introduces or deletes a nucleotide in one of the codons, the entire reading frame is shifted and the corresponding amino acid sequence changes. Consequently, the resulting protein sequence is different from the one encoded by the gene involved in the fusion. The gene fusion encodes a protein sequence that either corresponds to a completely unknown protein or contains a premature stop codon (the presence of a premature stop codon in the chimeric sequence interrupts the protein translation resulting in the truncation of the protein encoded by the 5’ fused gene). This class of mutations is functionally similar to nonsense point mutations that play a role in many cancers and might imply the loss of functionality of the 5’ fused gene. The sequence is labeled as in-frame if the number of nucleotides composing the fusion sequenced is a multiple of three and no premature stop codons are introduced in the chimeric sequence.

Figure 4

Frame check and domain annotation. Figure 4 (a) The length of the fusion transcript, from the start to the stop codon, must be multiple of three (three nucleotides per single encoded codon). If the length of the sequence module three is non-zero, the fusion sequence if frame-shifted. A premature stop codon can be introduced in the protein sequence. Figure 4 (b) The nucleotide sequence resulting from the fusion of the 5’ and 3’ gene is translated into amino acid sequence. Similarly, the genomic breakpoint coordinates are translated into protein amino acid coordinates. UniProt Web Service is queried and the list of the available domains for both the gene is retrieved. On the basis of the protein domain sequence and protein breakpoint, the list of both conserved and lost domains is reported.

The annotation of the preserved and lost protein domains is essential in order to capture the oncogenic potential of a translated chimeric transcript. The nucleotide fusion sequence assembled in the previous step is translated into an amino acid sequence. Subsequently, the UniProt web service [30] is queried for all available annotations of the putative protein encoded by the two genes involved in the fusion (Figure 4b). Leveraging the reading frame information and fusion breakpoint, Pegasus determines the conserved and lost domains associated with both 5’ and 3’ genes. It is worth emphasizing that both conserved and the lost domains are valuable features of a fusion transcript, with the former more likely to discriminate oncogene related fusions and the latter more likely to discriminate tumor suppressor related fusions.

The domain annotation permits the creation of a detailed feature space for the fusion transcripts, a prerequisite step for posing the ensuing machine learning task. In Pegasus the feature space is composed of:

• Binary information about reading frame and breakpoint region (if the breakpoint falls in coding regions, introns and UTRs);

• Presence or absence of ~1000 protein domains from UniProt. Our selection was based on the domains occurring in the training set from ChimerDB2.0.

• Number of oncogenic or tumor suppressor domains, as defined by association with the keywords “tumor suppressor” or “oncogene” in the UniProt database.

• Number of protein-protein oncogenic interacting domains. We check if one or more domains of the fusion interact with both oncogenic and tumor suppressor domains.

Driver fusion prediction as a binary classification task

We aim to fit a model that can identify oncogenic fusions from the background of passenger events and artifacts. More precisely, we aim to learn a mapping f : X→y from the fusion transcript feature space X to a label yϵ{0, 1} representing oncogenic driver status. Since we desire a biologically interpretable model that is also capable of capturing interactions between features, the decision tree is a natural choice. On the other hand, high dimensional feature spaces predispose to overfitting, and previous driver fusion prediction studies [27] focused a great deal on upfront dimensionality reduction for this reason. A single, large decision tree classifier is not likely to generalize beyond some training depth. An ensemble of shallower decision trees, if learned in a boosting framework, can guard against overfitting because of the iterative nature of the learning and the additive structure of the model [31]. Therefore the balance we strike between the expressive power of decision trees and robustness to overfitting comes in the form of stagewise additive modeling. In Pegasus we employ an additive model $f\left(x\right)={\displaystyle \underset{m}{\Sigma }{\alpha }_m{h}_m\left(x\right)}$ composed of weighted decision trees, an instance of which is shown in Figure 5, and fit via gradient boosting [32]. There is no manual reduction of features or feature space dimension in this strategy, unlike the manual selection of 6 enriched functional categories in the Oncofuse framework [27].

Figure 5

Decision tree base learner. The base learner in the boosted classification model is the decision tree, which defines a recursive partitioning of the feature space into disjoint regions each modeled by a constant. Decision nodes are colored brown while leaf nodes are colored green. At the leaf nodes, the samples have been partitioned into disjoint sets with a higher degree of label homogeneity than at the root.

Thus we require neither upfront dimensionality reduction schemes nor the restrictive assumption of class conditional feature independence in the Naïve Bayes model. Gradient tree boosting is an ensemble learning technique, wherein decision trees are used as base learners and the final model is expressed as an expansion in these basis functions. Figure 5 depicts a sample regression tree that would be added to the ensemble in a single round of boosting. Although the base learners are performing regression, appropriate choice of loss function for gradient boosting yields a classification task. Here we use the binomial deviance loss, and enforce a maximum depth of 5 nodes in the individual decision trees. The gradient tree boosting algorithm, originally published in 2000, is outlined below [33] and the implementation we use can be found in the scikit-learn python library [34]. We denote the number of training examples by N, the number of boosting rounds by M, and the loss function by ℒ.

$$\begin{array}{c}\hfill \left(a\right)\mathrm{For}i=1,2,\dots ,N\mathrm{compute}\hfill \\ \hfill \mathrm{Let}{y}_i\equiv \mathrm{class}\mathrm{label}\mathrm{of}\mathrm{transcript}i\hfill \\ \hfill \mathrm{Let}{f}_m\left(x\right)\equiv \mathrm{classification}\mathrm{function}\mathrm{at}\mathrm{boosting}\mathrm{round}m\hfill \\ \hfill \mathrm{Initialize}{f}_0\left(x\right)={\mathrm{argmin}}_{\mathrm{\gamma }}{\displaystyle \underset{i=1:N}{\Sigma }\mathrm{ℒ}\left(y_i,\gamma \right)}\hfill \\ \hfill \mathrm{For}m=1\mathrm{to}M:\hfill \end{array}$$

(1)

$$\begin{array}{c}\hfill \left(a\right)\mathrm{For}i=1,2,\dots ,N\mathrm{compute}\hfill \\ \hfill r_{\mathit{im}}=-{\left[\frac{\partial \mathrm{ℒ}(y_i,f\left(x_i\right)}{\partial f\left(x_i\right)}\right]}_{f=f_{m-1}}\hfill \end{array}$$

(2)

• (b) Fit a regression tree to r _im producing regions R _jm , j = 1, 2,…, J _m

• (c) $\begin{array}{c}\hfill \mathrm{For}j=1,2,\dots ,J_m\mathrm{compute}\hfill \\ \hfill {\gamma }_{\mathit{jm}}={\mathrm{argmin}}_{\mathrm{\gamma }}{\displaystyle \underset{x_iϵR_{\mathit{jm}}}{\Sigma }\mathrm{ℒ}\left(y_i,f_{m-1}\left(x_i\right)+\gamma \right)}\hfill \end{array}$

• (d) $\begin{array}{c}\hfill \mathrm{For}j=1,2,\dots ,J_m\mathrm{compute}\hfill \\ \hfill {\gamma }_{\mathit{jm}}={\mathrm{argmin}}_{\mathrm{\gamma }}{\displaystyle \underset{x_iϵR_{\mathit{jm}}}{\Sigma }\mathrm{ℒ}\left(y_i,f_{m-1}\left(x_i\right)+\gamma \right)}\hfill \end{array}$

An alternative ensemble classification strategy, the random forest algorithm, demonstrated comparable performance to gradient tree boosting in our experiments. There is recent precedent in the machine learning literature for initializing gradient tree boosting models with rankings learned via random forests for achieving superior performance to either algorithm alone [35]. Interestingly, those authors found that posing web-search rankings as a classification task rather than a regression task increased the performance of a gradient boosted regression tree model, confirming our hypothesis in constructing the current Pegasus classifier.

Results and Discussion

This section highlights the performance of Pegasus in detecting driver gene fusions. First, we examine the performance of the classifier on the training data and compare its effectiveness to a recently published tool, Oncofuse, on a separately curated validation dataset. Next, we run Pegasus on two experimental datasets and demonstrate its role in reducing the search space of potential oncogenic drivers by accurately ranking fusion transcripts from a vast set of putative candidates. The first is the publicly available RNA-Seq data of GBM from TCGA. The second is a non-public set of 23 RNA-Seq samples from a cohort of patients with ALCL, with 2 out of the 23 samples reporting the NPM1-ALK fusion.

We analyze these datasets with ChimeraScan or deFuse and apply Pegasus to the entire set of detected fusions. It is worth specifying that in the reported results about chimeric transcript annotations, if two or more fusions share the same junction breakpoint coordinates, they are counted as a single fusion. The rationale is that according to the Pegasus fusion domain analysis, if two genes fuse in different samples with the same breakpoint they also share exactly the same domain. Conversely, if two genes occur in different samples with different junction breakpoint coordinates, the domain analysis accordingly changes.

Classifier performance on training corpus and independent validation set

The corpus of labeled data used to train the classifier comes from two sources. Positive examples, meaning true oncogenic driver fusions, are drawn from ChimerDB2.0, which contains 501 curated driver fusions. 1500 negative examples are then drawn from an internal collection of reactive lymph node tissue in patients with no clinical history of malignancy. The negative examples contain passenger fusions as well as read-through transcripts. We also supplement the negative training data with 416 deliberately frame-shifted transcripts from ChimerDB2.0 such that the necessary driver domains are lost. In total there are 501 positive examples and 1916 negative examples in the training corpus. The rationale for augmenting the negative set with 416 frame-shifted fusions from ChimerDB2.0 is to include the scenario of chimeric transcripts containing an oncogene at the 3’ position that is frame-shifted. Since such events occur at low frequency in normal lymph node tissue, this design choice improves the performance of the classifier (for a detailed discussion please refer to Additional file 1). In summary, the 501 fusions from ChimerDB2.0 form the positive training set and provide mostly in-frame fusions involving oncogenes. The 1500 fusions from normal tissue contribute to the negative set and provide both in-frame and frame-shifted fusions. The 416 deliberately frame-shifted fusions from ChimerDB2.0 complete the negative set and provide frame-shifted gene fusions mostly with an oncogene at the 3’ position.

The classifier is trained for 100 rounds of boosting under 10-fold stratified cross validation (CV) and achieves a mean test split AUROC of 0.96. As expected, in Figure 6b the loss on the train split monotonically decreases with increasing model complexity, while we see no sign of overfitting in the form of rising loss on the test split. Since each decision tree base learner implies a hierarchy of informative features, we can average over the boosting rounds to produce an aggregate view of the most important features in the classification task. Specifically, the relative feature importances in Figure 6a are computed via an ensemble average of the single decision tree feature importances as defined in Breiman et al.[36]:

$$\begin{array}{l}\mathrm{Let}{I}_k^2\equiv \mathrm{squared}\mathrm{importance}\mathrm{of}\mathrm{feature}k\hfill \\ \mathrm{Let}M\equiv \mathrm{number}\mathrm{of}\mathrm{boosting}\mathrm{rounds}.\hfill \\ \mathrm{Let}T\mathrm{be}\mathrm{a}\mathrm{decision}\mathrm{tree}\mathrm{with}J–1\mathrm{pairs}\left(t,{î}_{t}2\right)\mathrm{o}f\hfill \\ \left(\mathrm{internal}\mathrm{decision}\mathrm{node},\mathrm{estimated}\mathrm{improvement}\mathrm{in}\mathrm{squared}\mathrm{error}\mathrm{risk}\right).\hfill \\ I_k^2\left(T\right)={\displaystyle \underset{t=1:J–1}{\Sigma }{î}_{t}2\mathbb{I}\left(t=k\right).}\hfill \\ \Rightarrow I_k^2=\frac{1}{M}{\displaystyle \underset{m=1:M}{\Sigma }{I}_k^2\left(T_m\right)\text{.}}\hfill \end{array}$$

(3)

Figure 6

Classifier training under stratified 10-fold cross validation. Because of the imbalance in positive/negative training examples we stratify on the label and ensure that each 10% of data held out as a test set contains equal proportions of both labels. Figure 6 (a) Relative feature importances averaged over the 10 CV folds. Figure 6 (b) Binomial deviance loss function plotted against model complexity for both the train and test splits of the data.

We observe that the computationally expensive step of computing the fusion transcript reading frame is justified in the eyes of the classifier, as it is the single most informative feature. Looking a little further down the list we learn other transcript features that are highly informative of driver events, such as having breakpoints in the CDS and conserving domains shared with or interacting with known oncogenes.

Despite strong performance of the model under 10-fold CV on the training corpus, we are interested to see whether the classifier can generalize to new fusion transcripts that were unseen during the training phase. A list of 39 driver fusions, the majority of which are more recent than ChimerDB2.0, and corresponding to the validation set used in [27], is adopted as the positive validation set examples. To balance the label frequencies, we also select 39 transcripts from benign, reactive lymph node tissue as the negative validation set examples. None of the 78 validation examples are included in the training data. The negative examples are selected to contain at least one oncogene or tumor suppressor domain with the rationale that such transcripts more closely resemble driver fusions and would be most challenging for a classification function. In Figure 7a we demonstrate the favorable performance of the trained Pegasus classifier versus Oncofuse, the current state of the art in data-driven prediction of driver fusions. Since the ROC curve does not necessarily reveal how well separated the Pegasus scores are for the two class labels, we include Figure 7b to illustrate the remarkable resolution the classifier achieves between positive and negative examples. We also verify that Pegasus outperforms Oncofuse on randomly drawn sets of 39 non-oncogenic transcripts, though by a smaller margin (ROC curve in Additional file 2: Figure S1). This is to be expected because the majority of non-oncogenic fusions are very easily classified, whereas our curated subset represents a more challenging task. Such robust performance on manually curated data sets naturally leads to the next proving ground, applying Pegasus to the enormous candidate lists generated from real RNA-Seq samples.

Figure 7

Classification performance of trained model. A curated set of 39 recently reported fusions that were not present in the ChimerDB2.0 training corpus are now used as an independent validation set. Figure 7 (a) Superior classification performance of Pegasus compared to Oncofuse on the validation dataset. Figure 7 (b) Boxplot demonstrating the high resolution of the Pegasus driver score in discerning the class boundary. Positive examples are consistently scored near 1 and negative examples consistently near 0.

Pegasus driver fusion predictions in non-public GBM data

In order to demonstrate the effectiveness of Pegasus in predicting driver fusions, we analyze 15 samples of short-term glioblastoma stem cells freshly isolated from individuals with GBM. RNA-Seq samples were first analyzed with ChimeraScan and deFuse [14],[15] for fusion detection. Next, we apply Pegasus to the set of gene fusion candidates and consider as driver events all those fusions having a number of supporting reads greater than 10 and a Pegasus Driver Score (PDS) greater than 0.8. As shown in Figure 7(b), a threshold of PDS > 0.8 promises a good trade-off between specificity and sensitivity. Table 1 reports the 4 detected driver fusions. All fusions have been validated with RT-PCR (see Additional file 3: Figure S2) yielding a 100% rate of transcript validation. And while recurrence is often the surrogate measure of functional importance, the four unique driver candidates from this 15 patient sample still contain features associated with oncogenic gene fusions. In fact, CAND1-EGFR has been reported in [37] and the EGFR gene has been demonstrated to have an oncogenic role in GBM. Moreover, fusions involving MAPK1 and VOPP1 are reported as frequent in GBM with different gene partners [37],[38]. These results show that Pegasus can successfully detect relevant driver fusion candidates from RNA-Seq data and that a threshold of PDS > 0.8 and number of supporting reads greater than 10 provide a strong transcript validation rate.

Table 1 Pegasus predictions on 15 private GBM RNA-Seq data

Pegasus driver fusion predictions in public TCGA GBM data

As the most common and deadly primary brain cancer, GBM has recently undergone a deep investigation by the multi-institutional consortium, the cancer genome atlas (TCGA). TCGA makes its collected RNA-Seq data available to the larger scientific community, and we analyze a set of 161 samples from their GBM cohort. We first analyze the 161 RNA-Seq samples with ChimeraScan (default parameters) [15], detecting a total of 9349 unique fusions across the entire dataset. Next, we apply Pegasus to the set of candidates and consider as driver events all fusion transcripts having a number of supporting reads greater than 10, Pegasus Driver Score (PDS) greater than 0.8 and recurrence greater than 1 (for the complete list see Additional file 4: Table S1). As shown in the non-public RNA-Seq data, these filtering thresholds provide a good validation rate by RT-PCR. The application of these filters reduces the original list of 9349 candidates down to 13 high confidence fusions, making further functional analysis and validation tractable. Pegasus computes a score greater than 0.8 for both FGFR3-TACC3 and EGFR-SEPT14 gene fusions, which are already reported as driver translocation events in GBM [10],[37]. However, since TCGA biological material is not available, we are unable to perform further functional analysis of all predicted driver fusions with experimental procedures. Nonetheless, in order to validate Pegasus performance we compare PDS values with the frequencies reported in both Frattini et al. and Brennar et al. (see Table 2). Of the 13 high confidence Pegasus predictions, 6 are recurrent in TCGA data suggesting a potential functional driver role in GBM tumorgenesis [39]. Some of the recurrent fusions involve the EGFR gene that is usually amplified in GBM, where it is known to activate STAT3 signaling and is thus a drug target. Particularly interesting is also the BCAN-NTRK1 gene fusion. In fact, NTRK1 is often translocated with different partners in cancers beyond just GBM [37].

Table 2 Pegasus predictions on GBM RNA-Seq data

Functional validation of new recurrent driver fusion in anaplastic large cell lymphoma

Anaplastic large cell lymphoma (ALCL) is a form of peripheral T-cell lymphoma that is often associated with translocations of the ALK gene. In 23 non-public ALCL samples (~450 million properly mapped reads) we detect a total of 5201 candidate fusion transcripts by means of deFuse [14] and ChimeraScan [15]. Beyond the two NPM1-ALK fusion transcripts (PDS = 0.98) that are already reported, Pegasus properly annotates and reveals 16 new biologically relevant fusions in these 23 samples. All 16 candidate driver fusions have been validated with RT-PCR, and 4 gene fusions have successfully undergone functional assays and in vivo validation. An example of Pegasus’ effectiveness in functional domain analysis lies in the oncogenic role of TRAF1-ALK [40], a novel fusion in ALCL that Pegasus reports as driver.

The TRAF1-ALK fusion has been reported in three cases of ALCL (one in [40] and two in Abate et al., manuscript under review) suggesting a driver role. Pegasus accurately assembles and annotates the in-frame fusion sequence and correctly detects that the ALK protein kinase domain is completely conserved. Various fusions involving the ALK gene have been reported in literature and the oncogenic effect is generally promoted by ALK signaling. Interestingly, TRAF1 is also known to be involved in both the canonical and non-canonical NFkB pathway. Pegasus correctly annotates that the meprin and TRAF-C homology (MATH) domain is conserved, a domain that ubiquinates the IKK complex, activating NFkB transcription factors. As depicted in Figure 8, Pegasus properly identifies both the presence of an oncogenic protein domain (ALK) and an interacting oncogenic domain (TRAF1 to NFKB activation). Experimental work demonstrates and validates this Pegasus prediction showing the oncogenic effect of TRAF1-ALK in vivo, with activation of both ALK and NFkB signaling (Abate et al., manuscript under review).

Figure 8

Novel driver TRAF1-ALK gene fusion in ALCL. A graphical representation of the Pegasus annotation on the TRAF1-ALK fusion in ALCL. Conserved domains are reported according to the junction breakpoint.

Conclusion

Since the first application of whole transcriptome sequencing to gene fusion discovery in 2009 [8], many new aberrant events have been reported, opening an exciting frontier for molecular understanding of cancer biology and targeted therapies. The unprecedented sensitivity of NGS technology, however, often yields numbers of fusion candidates too large to be experimentally validated. The frontier is made ever more exciting (and challenging) by large consortia such as TCGA and ICGC (International Cancer Genome Consortium) who are making available large sets of RNA-Seq samples spanning the spectrum of human malignancies. The analysis of this data has been revealing the limits of the theory that associates driver events with recurrent events. In fact, out of 161 RNA-seq GBM samples, the most frequent fusion (EGFR-SEPT14) occurs in only 6 samples, and the highly expressed FGFR3-TACC3 fusion is recurrent in only the 3% of GBM cases. Thus, this data suggest that in order to select relevant driver fusion candidates for biological validation, a functional analysis of the putative gene fusion candidate is necessary.

Here we present Pegasus, an accurate prediction tool for the discovery of new driver gene fusions in cancer studies. The proposed methodology is based on a computational model of the features that make chimeric transcript a driver oncogenic event. The framework provides a common interface for several fusion detection tools and it predicts driver events by properly analyzing the detected gene fusion candidates according to the assembled fused sequence.

The application of ensemble learning techniques reveals the most informative features in discriminating oncogenic gene fusions. The data confirm our intuition that an accurate analysis of fusion transcript sequence is necessary. The reading frame in particular is a dominant features in the discrimination of passenger and driver fusions. Similarly, the molecular characterization of the main reported oncogenic domains accurately increase the sensitivity and the PDS computation.

The problem of computationally assessing the biological and clinical relevance of a gene fusion is still very much an open question. However, some driver prediction tools have been recently proposed. To better determine Pegasus performance and accuracy, we compare our predicted results with Oncofuse. The data confirm that an approximate analysis of the fusion transcript sequence negatively impacts the performance of the algorithm. Using a set of known driver fusions as positive examples and a set of passenger fusions from normal tissue as negative examples, we observe the superior performance of Pegasus in ROC space where its AUC is 0.97.

Moreover, we demonstrate the practical role of the Pegasus framework in computing PDS scores that allow for triaging lists of gene fusion candidates for experimental validation in two actual case studies. In the first, we apply Pegasus to public GBM TCGA data and almost 50% of the detected driver fusions turn out to have been reported as a biologically relevant in independent studies [37],[38]. In the second, we compute Pegasus scores for internal ALCL data and we successfully detect novel oncogenic and targetable driver fusions that have undergone complete functional and experimental validation. In this work, we extensively report the driver fusion TRAF1-ALK that has been correctly detected and highly ranked by Pegasus.

Finally, the accuracy of Pegasus in detecting driver fusions in both the curated validation dataset and the real biological cases demonstrates the efficacy of the framework in supporting biological analysis and cancer research. We believe that the Pegasus prediction score, as well as the accurate annotations provided via our feature engineering, will be of great use to other investigators searching for biologically relevant gene fusions in NGS data.

Availability and Requirements

Project name: Pegasus

Project home page: http://sourceforge.net/p/pegasus-fus.

Operating system: UNIX

Programming language: Java, Perl, Python, BASH.