Fusion-sq: detecting tumor-specific gene fusions with high confidence
To resolve gene fusions and investigate their relevance to pediatric cancer, we combined RNA-seq and paired tumor/normal WGS samples of 129 patients across 53 pediatric cancer types (Fig. 1a). Known clinically relevant fusions were identified with RNA-seq in 30 patients (Fig. 1a, gray overlay). We investigated combining WGS structural variant analysis with RNA-seq data to increase detection specificity and identify potentially pathogenic gene fusions in the remainder of the cohort. Hereto, we developed Fusion-sq which integrates (“fuses”) predicted gene fusions from RNA-seq data with SVs from WGS data (Methods, Fig. 1b,c). For every predicted gene fusion, Fusion-sq first derives genomic intervals to match RNA and DNA breakpoints based on intron-exon gene structure. Next, DNA breakpoints falling within these intervals are used to identify SVs that link the fusion 5’ and 3’ partner genes. To optimize both recall and precision, gene fusion predictions from STAR-Fusion and FusionCatcher are integrated with SVs detected by Manta [23], DELLY [24] and GRIDSS [25]. We selected high confidence fusions that are 1) predicted by both RNA gene fusion tools and 2) supported by SVs detected by at least two callers (Fig. 1b,d, Additional file 1: Figure S1). Fusions are then classified as tumor-specific, (likely) germline or low allele fraction (AF) based on the SV’s variant supporting reads in the paired tumor/normal samples, as reflected in the AFs of SVs in both samples (Fig. 1b,d). In total, 2906 fusions were predicted in 129 patients using RNA-seq data alone (Fig. 1b). By combining these fusion predictions with WGS data, Fusion-sq identified 232 high confidence fusions in 83 patients (Fig. 1b,e). To investigate potentially pathogenic gene fusions, we further analyzed the 155 high confidence tumor-specific fusions (hcTSF) by combining RNA evidence with properties of the underlying SVs, such as SV type and AF.
a Distribution of 129 pediatric cancer patients colored according to ICD-O-3 primary group. Gray overlay indicates patients for which a clinically relevant fusion was identified by RNA-seq. b Schematic overview of the number of fusions at different steps throughout the Fusion-sq pipeline. Subsets discussed extensively in the main text are highlighted in red and available in Table 1. Note that the 35 clinically relevant gene fusions include reciprocal orientations. In addition, slight differences in numbers with the main text arise from excluding patients without known clinically relevant fusions, e.g the 27/30 high confidence tumor-specific fusions (hcTSFs) involving oncogenes (onco) or tumor suppressor genes (TSG).
c Schematic overview of the Fusion-sq algorithm to find SVs that support fusion predictions by linking the upstream (5’, blue) and downstream (3’, red) partner genes. First, genomic intervals to match RNA-DNA breakpoints are derived based on the intron-exon gene structure and the RNA breakpoints. Next, DNA breakpoints located in these matching intervals are used to identify SVs that link together the partner genes (Supplementary Methods). d Upper panel: schematic representation of high confidence fusion detection based on SVs identified by two or more tools (> 50% reciprocal overlap). Lower panel: scatter plot of tumor and normal AF, resulting in classification of fusions in tumor-specific (blue), likely germline (red) and low AF (green). e Number of fusions in individual patients: hcTSFs (circles, grouped and color-coded by primary cancer type group) and fusion predictions by both RNA fusion tools with FFPM > 0.1 (gray bars). Fusion status is indicated by a black square for clinically relevant fusions and a circle for oncogene/tumor suppressor gene fusions.
All clinically relevant fusions match to tumor-specific SVs
We first focused on the subset of 30 patients known to carry a clinically relevant gene fusion. In all patients, we resolved the tumor-specific SVs for the predicted gene fusions (Table 1) which helped to reconcile variability in RNA-based predictions. To better understand gene fusions at the genomic level, we investigated the underlying SVs supporting the chimeric transcripts. For all patients but one, the associated duplications (DUP), deletions (DEL) and inversions (INV) were almost identical (> 99% overlap) and breakpoints of intra-chromosomal translocations (CTX) were resolved within 10 base pairs (bp) showing a high agreement between SV callers (Fig. 2a). The distances between breakpoints in RNA from the chimeric transcripts and corresponding breakpoints in WGS data from the underlying SVs are highly variable amongst gene fusions and patients (Fig. 2a). Notably, for nine patients the distance between corresponding RNA-DNA breakpoints is larger than 10 kilobase pairs (kb) (Fig. 2a, red lines), illustrating the advantage of using intron-exon gene structure to define genomic intervals for matching chimeric transcripts with SVs.
a Distance between RNA breakpoints and matching SVs as resolved by DELLY (circle), Manta (square) and GRIDSS (triangle) for the 5' and 3' partner genes in 30 patients carrying clinically relevant fusions. Colors represent specific gene fusions. Red lines indicate 10 kb fixed-size matching intervals and would fail to match SVs for nine gene fusions (labeled). Note that all clinically relevant fusions are detected by at least two SV tools at nucleotide resolution (< 10bp, overlapping symbols) except for ASPSCR1–TFE3. b Gene fusion predictions supported by both RNA fusion tools (all circles, same patients as in a) with their read support in fusion fragments per million total RNA-seq fragments (FFPM). The red line indicates the default cutoff below which fusions are usually discarded (low expression FFPM < 0.1). Supporting high confidence tumor-specific SVs were identified for all clinically relevant fusions (colored circles; colors same as in a) and seven additional fusions (black) of unknown significance, but not for the remaining RNA-only fusion predictions (gray). Colored bar at the x-axis indicates the primary cancer type group.
One clinically relevant fusion that proved difficult to resolve was ASPSCR1–TFE3, as the underlying translocation t(7;X) was identified differently by Manta, DELLY and GRIDSS (Table 1). Manta resolved it as a composite fusion of a CTX + INV in chrX, whilst GRIDSS and DELLY had unusually broad CTX footprints of 393 bp and 1181 bp respectively. The associated SV breakpoints of all three tools overlapped SINE elements, which may indicate difficulties in resolving this gene fusion to base pair accuracy in DNA. Despite this potentially complex variant, this translocation could still be mapped to the ASPSCR1–TFE3 chimeric transcript by Fusion-sq.
Simple SV events underlie the remaining (n = 29) clinically relevant fusions [34] and represent all major SV types, such as deletions, duplications, inversions and translocations. As expected, inter-chromosomal translocations (CTX, n = 18) were the most common SV type underlying clinically relevant gene fusions [1]. In some cases, these translocations support both the canonical and reciprocal transcripts. As an example, we identified a reciprocal gene fusion product for three out of four patients with an ETV6–RUNX1 fusion, and for one out of four patients with an EWSR1–FLI1 fusion. In nine patients, we identified duplications resulting in KIAA1549–BRAF fusions. Interestingly, the KIAA1549–BRAF fusion of patient M218AAA was identified as an inversion by all three SV callers, and it is likely an inverted duplication as read depth is also increased (+ 0.42 copy number log2 fold change, CN l2fc). The two remaining clinically relevant fusions are FGFR1–TACC1 caused by a 420 kb INV, and LMNA–NTRK1 caused by a 740 kb DEL, showing that a variety of SV types can result in clinically relevant fusion events.
To reduce false positives in RNA-seq fusion detection, filtering based on read support for chimeric transcripts is often implemented with a default minimum of 0.1 fusion fragments per million total RNA-seq fragments (FFPM) [2, 35]. However, by relying only on RNA evidence, four clinically relevant fusions would be missed due to their low expression (Fig. 2b, Table 1), two KIAA1549–BRAF fusions (0.08–0.09 FFPM), one ETV6–RUNX1 (0.07 FFPM) and one PAX3–FOXO1 (0.03 FFPM). These lowly expressed fusions can be discerned from false positives by integration with WGS. In total, 715 fusions were predicted by both STAR-Fusion and FusionCatcher in these 30 patients alone, of which 147 passed the read support threshold (FFPM > 0.1). In contrast, for 24 patients their clinically relevant gene fusions are the only hcTSFs indicating a high specificity of Fusion-sq. In the remaining six patients, Fusion-sq resolved an additional seven hcTSFs with similar support but of unknown clinical relevance (Fig. 2b, black dots). This shows that integration of RNA-seq and WGS by Fusion-sq can accurately resolve tumor-specific fusions, effectively removing the need for extensive manual filtering to select fusions for follow-up.
Underlying SVs distinguish tumor-specific fusions from healthy chimera
The 232 high confidence fusions (hcFs) identified in 83 patients were classified as tumor-specific, germline or low AF (Fig. 1b). Both the clinically relevant and additionally detected hcTSFs have a higher tumor AF relative to the normal AF, whereas germline fusions and low AF fusions show similar tumor and normal AFs (Additional file 1: Figure S2). While it is counter-intuitive to have high confidence variants with these low AFs, we reasoned that a high number of variant and reference reads could explain this. Indeed, 81% of these low AF hcFs originated from amplified regions (CN l2fc > 1.58). Next, we evaluated the efficacy of identifying hcTSFs by 1) assessing the underlying SV properties, 2) annotating with databases of chimera and SVs, and 3) annotating with cancer-related genes.
We then investigated whether the underlying SVs of additionally detected gene fusions resemble those of known clinically relevant fusions. Hereto, we mapped the high confidence fusions resolved in individual patients to distinct fusions to account for recurrent fusions. In total the 232 hcFs mapped to 189 distinct fusions of which 134 tumor-specific, 18 germline and 37 low AF (Table 2, Additional file 1: Figure S3). Like the clinically relevant fusions, interchromosomal translocations (CTX, 39 of 119) are the most common SV type underlying the additionally identified tumor-specific fusions. The remaining tumor-specific intra-chromosomal SVs are distributed over DUP, DEL and INV (Fig. 3a). In contrast, germline fusions are depleted in CTX events and generally caused by shorter intra-chromosomal SVs. The tumor-specific SVs are approximately 3.5x larger than germline SVs (Fig. 3b). Only a single germline fusion (CCDC32–CBX3) has an underlying CTX. CCDC32–CBX3 is a known healthy chimera and the SV breakpoints overlap ALU repeat elements which indicates either potential mapping difficulties or a mechanism involved in forming this event. In conclusion, these results show that the types and sizes of SVs underlying tumor-specific gene fusions are distinct from germline events and resemble the SVs of known clinically relevant fusions.
Number of (a) and average SV length in base pairs (bp) for (b) distinct high confidence fusions categorized either as tumor-specific (upper panel), likely germline (middle panel) and low AF (lower panel). Each category is further subdivided according to SV type: inversion (INV), duplication (DUP), deletion (DEL) or intra-chromosomal translocation (CTX). Fusions are colored based on known clinical relevance (clinical; purple), involving a cancer-related gene or cancer chimera (cancer; red), population SV or healthy chimera (common; blue), or both cancer and common (both; orange).
Secondly, to compare hcTSFs with population variants, the gene fusions were annotated with databases cataloging healthy chimeric transcripts or databases of SVs occurring in the general population as normal variation (Table 2), such as NCBI Common SV database [26], gnomAD [27] and DGV [28]high confidence fusions, ten are flagged as a healthy chimera by the gene fusion caller, of which seven were also classified as germline fusions. Further comparison of SVs underlying fusions with SVs occurring in the general population showed overlap for eight germline fusions and five tumor-specific fusions. Overall, 56% of the distinct germline hcFs either overlap a population SV or occur in a healthy chimera database, compared to only 4% of tumor-specific fusions (Fig. 3a). These results indicate that the hcTSFs are depleted for germline population variants.
Finally, to identify potentially pathogenic fusions, we compared the results of Fusion-sq with chimera previously detected in cancer and annotated fusions with cancer-related genes. 44 of the 189 distinct hcFs are present in ChimerDB [30] or the Mitelman database [1]. Most of these are classified as tumor-specific (28), but also germline (8) and low AF (8) fusions occur in these databases. Six of the eight germline fusions annotated as cancer chimera are also annotated as either healthy chimera or have underlying intra-chromosomal SVs that overlap with general population SVs (Fig. 3a). This corresponds to previous observations that cancer chimera databases can include passenger fusions [36, 37]. After accounting for recurrence, 128 tumor-specific fusions were identified which have no evidence of occurring in the normal population. Next, we annotated the fusion partner genes as proto-oncogene or tumor-suppressor gene (TSG). None of the germline hcFs affect an oncogene or TSG, further substantiating the classification based on AF of the underlying SVs (Fig. 3a). In summary, we resolved 27 distinct hcTSFs involving an oncogene and/or TSG with similar high confidence as the known clinically relevant fusions, which therefore provide sufficient evidence for further investigation.
Resolving tumor-specific fusions in individual genomes
Focusing on patients without a known clinically relevant fusion, we identified 113 hcTSFs of unknown significance to investigate further. These fusions are detected in 32 patients across different cancer types. Generally, each individual patient carries few to none hcTSFs (Fig. 1e), whereas a high burden of gene fusions is associated with copy number instability (Fig. 4. Additional file 1: Figure S4). High fusion burden ( > = 5 hcFs, ~ 95th percentile) is associated with a high fraction of genome altered (FGA) by copy number alterations (CNAs) (median 63% vs 3.8%, mean 48% vs 11%, wilcox p < 0.01) in tumor types prone to copy number instability such as osteosarcoma (4x), embryonal sarcoma, embryonal rhabdomyosarcoma, neuroblastoma and ependymoma [33].
Circos plots of osteosarcoma patient M691AAA (a), neuroblastoma patient M787AAA (b) and embryonal rhabdomyosarcoma patient M002AAB (c). The plots are annotated with high confidence tumor-specific gene fusions (multi-colored links) and copy number gains (red) and losses (blue). d Relationship between gene expression changes of the 5' (left panel) and 3' (right panel) partner gene relative to the cancer type supergroup (FPKM z-score, color scale displayed on the right) and copy number (CN) log-2 fold changes in read depth. Genes located inside amplified regions (CN l2fc > 1.58) are displayed with larger circles. Cancer-related genes are labeled according to their patient of origin.
To further investigate this association between copy number instability and high fusion burden, we closely studied four osteosarcoma patients carrying TP53 and ATRX fusions (Additional file 1: Figure S5). These patients’ tumors have many dispersed CNAs (median FGA 67%) and many hcFs distributed across their genomes (range 5–15, Fig. 1e, Additional file 1: Figure S5a). In three of these osteosarcoma patients, we resolved fusions involving the first exon of TP53 and in each case a different 3’ partner gene resulting from an inversion and two translocations; t(17;6), t(17;20) (Additional file 1: Figure S5b) (Table 1). Translocations with the first exon of TP53 have previously been identified as cancer driver events in osteosarcoma [38]. Therefore, these fusions are potentially also pathogenic in these patients, especially since a driver mutation had yet to be identified for these tumors. Apart from the TP53 gene fusions and associated SVs, we did not observe additional CN losses or deleterious SNVs in the TP53 loci of these patients. In the fourth osteosarcoma patient (M691AAA, Fig. 4a), we detected an ATRX fusion which has been suggested previously as a potential driver mutation for osteosarcoma [38]. For the ATRX fusion, gene expression concomitantly was reduced relative to the cancer type supergroup (-3.3 z-score FPKM (zfpkm), p < 0.01). Similarly, the group of patients with a TP53 fusion showed reduced expression relative to the cancer type supergroup (0.88 vs 2.0 log2 FPKM, p < 0.05) (Additional file 1: Figure S5c). However, this was not clear for the individual patients, illustrating that the underlying SV provides additional evidence for a disruptive SV event that could not be easily derived from RNA-seq alone.
Multiple fusions originating from highly amplified regions were resolved in two patients with neuroblastoma (M787AAA, Fig. 4b) and embryonal rhabdomyosarcoma (M002AAB, Fig. 4c). The neuroblastoma patient has a focal amplification in chromosome 12q13-15 involving the oncogenes MDM2 and CDK4 (4–6 copy number log2 fold change, CN l2fc). In this region, we identified 18 high confidence, low AF gene fusions. Including fusions previously identified as cancer-related chimera (i.e. FRS2–MDM1, FRS2–PTPRR, CCT2–BEST3, RAB3IP–BEST3) of which both partner genes are overexpressed due to the amplification (Fig. 4d). Similarly, patient M002AAB carries 11 fusions originating from a focal amplification in chr11q22 (5–7 CN l2fc) which are detected with high confidence but some with a low AF. Here also, the fusion partner genes are overexpressed which is consistent with the amplification (Fig. 4d). Although chr11q22 amplification itself is not known to be clinically relevant, we resolved multiple fusion combinations with oncogenes BIRC3, PGR and in particular YAP1. YAP1–CFAP300 involves exons 1–5 of the YAP1 oncogene which is highly amplified (7.2 CN l2fc) and overexpressed (5.3 zfpkm, p < 0.1). Also, this exon 1–5 fragment has previously been identified as pathogenic when fused to other 3’ partner genes by activating the TEAD pathway [39, 40].
These case studies illustrate that Fusion-sq can confidently resolve fusions in unstable genomes with a high FGA or complex alterations. For some tumor-specific fusions, the underlying SVs are potentially pathogenic (e.g. the TP53 and ATRX fusions) while in other cases the fusions seem to be the result of copy number instability. In addition, gene fusions in amplified regions can exhibit high expression of their partner genes due to the underlying CN gain. Therefore, characteristics of underlying SVs are key for interpreting potential functional effects of individual gene fusions, especially in patients with unstable genomes.
Tumor-specific fusions affecting gene expression
Gene fusions can activate oncogenes or disrupt tumor-suppressor genes (TSG) and the resulting pathogenic effects can be reflected in dysregulation of gene expression. To identify potentially pathogenic fusions, we assessed the functional effects of the 27 distinct hcTSFs involving an oncogene or TSG (in bold) (Table 1). Their SVs represent all the simple SV types: CTX (10), DUP(7), DEL(5) and INV(5). As a proxy for functional effect, we combined expression data, underlying SVs and gene annotation.
Overexpression of oncogenes resulting from gene fusions is often suggested as an activation mechanism [2]. However, we did not observe an enrichment of oncogenes amongst overexpressed fusion partner genes (> 1.96 zfpkm) compared to the cancer type supergroup or the full cohort. Instead, our data suggests an association with copy number gain (CN l2fc > 0.58) irrespective of gene annotation, as reflected by approximately half of the fusions with CN gain found to be overexpressed (24 of 42). We identified fusions in individual patients for which the 3’ oncogenes are significantly overexpressed relative to the cancer type supergroup (p < 0.1): PAX3–WWTR1 (2.1 and 3.0 zfpkm) and SYMPK–MEF2B (2.1 zfpkm) (Fig. 5). In addition, we identified a MED14–HOXA9 fusion and associated overexpression of HOXA9 (2.2 zfpkm, p = 0.15) in a pre-T-cell lymphoblastic leukemia patient (M385AAA) (Supplementary Results).
Gene expression changes for the 5' (left) and 3' (right) partner genes of tumor-specific gene fusions involving oncogenes and/or tumor-suppressor genes (TSG). Changes are relative to the cancer type supergroup (FPKM z-score; color scale on the right). Individual 5' or 3' genes are marked as either oncogene (circle), TSG (triangle), both TSG and oncogene (square) or kinase (k). Fusions were divided into ‘activating’ or ‘disruptive’ functional effect categories based on partner gene annotation, and marked with an asterisk if they originate from patients with a high gene fusion burden.
Gene fusions can also be pathogenic through activation of kinases due to e.g. loss of an auto-inhibitory domain or increased dimerization. We identified eight hcTSFs involving kinases, of which three in patients with unstable genomes (QKI–MAPK3K4, PRKD3–LTBP1, STK40–OSCP1) and three fusions (ERBB4–LINC01807, FBXW7–DCLK2 and IGSF3–AKT3) with unknown relevance since the kinase domains are not part of the chimeric product. However two fusions resulted in a chimera containing the kinase domains, MEF2A–IGF1R and TNK1–GPS2. The TNK1–GPS2 fusion was found in a patient with anaplastic large cell lymphoma resulting from a 74 kb inversion (0.40 tumor AF). Resolving the underlying SV confirmed the presence of this gene fusion which was predicted with different breakpoints and in reciprocal orientations by the two RNA-seq tools (Supplementary Results).
Finally, we identified three gene fusions involving tumor-suppressor genes (TSGs) which showed a significant decrease in gene expression (p < 0.1) relative to the cancer type supergroup: ZBTB20–LSAMP (-4.0 zfpkm, p < 0.05), NF1–RAB11FIP4 (-2.4 zfpkm) and the previously mentioned ATRX–LINC01280 (-3.5 zfpkm) in a patient with osteosarcoma (Supplementary results).
These examples illustrate that combining WGS structural variant identification with RNA-seq based gene fusions results in the identification of novel high confidence gene fusions with strong biological relevance.