Exploration of the gene fusion landscape of glioblastoma using transcriptome sequencing and copy number data

Background RNA-seq has spurred important gene fusion discoveries in a number of different cancers, including lung, prostate, breast, brain, thyroid and bladder carcinomas. Gene fusion discovery can potentially lead to the development of novel treatments that target the underlying genetic abnormalities. Results In this study, we provide comprehensive view of gene fusion landscape in 185 glioblastoma multiforme patients from two independent cohorts. Fusions occur in approximately 30-50% of GBM patient samples. In the Ivy Center cohort of 24 patients, 33% of samples harbored fusions that were validated by qPCR and Sanger sequencing. We were able to identify high-confidence gene fusions from RNA-seq data in 53% of the samples in a TCGA cohort of 161 patients. We identified 13 cases (8%) with fusions retaining a tyrosine kinase domain in the TCGA cohort and one case in the Ivy Center cohort. Ours is the first study to describe recurrent fusions involving non-coding genes. Genomic locations 7p11 and 12q14-15 harbor majority of the fusions. Fusions on 7p11 are formed in focally amplified EGFR locus whereas 12q14-15 fusions are formed by complex genomic rearrangements. All the fusions detected in this study can be further visualized and analyzed using our website: http://ivygap.swedish.org/fusions. Conclusions Our study highlights the prevalence of gene fusions as one of the major genomic abnormalities in GBM. The majority of the fusions are private fusions, and a minority of these recur with low frequency. A small subset of patients with fusions of receptor tyrosine kinases can benefit from existing FDA approved drugs and drugs available in various clinical trials. Due to the low frequency and rarity of clinically relevant fusions, RNA-seq of GBM patient samples will be a vital tool for the identification of patient-specific fusions that can drive personalized therapy. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-14-818) contains supplementary material, which is available to authorized users.


Background
Cancers result from the accumulation of genomic mutations and epigenetic alterations that change gene expression and function. In particular, gene fusions have been recognized as an associated and significant feature of cancer since the characterization of the Philadelphia chromosome [1]. The occurrence of gene fusions in solid tumors has long been noted, but their importance has been appreciated only recently, largely due to high throughput technologies such as transcriptome sequencing (RNA-seq) [2][3][4][5]. RNA-seq permits genome-wide transcription analysis for novel transcript discovery.
RNA-seq has spurred important gene fusion discoveries for a number of different cancers, including lung [6,7], prostate [3,8,9], breast [10][11][12], brain [13], thyroid [14] and bladder carcinomas [15]. One obvious benefit from gene fusion discovery is the potential to develop novel treatments that target these genetic abnormalities. The EML4-ALK translocation fusion is an example in which the fusion causes constitutive kinase activity. Mouse fibroblasts transfected with EML4-ALK formed tumors when this fusion was injected into nude mice, thus demonstrating the oncogenic activity of the resultant protein [6]. Crizotinib, a competitive inhibitor of ALK, has recently been granted FDA approval for the treatment of specific late-stage, non-small cell lung cancers, and presently there are two phase 3 trials in progress [16].
Glioblastoma multiforme (GBM), a grade IV astrocytoma, is the most common form of primary brain cancer, with a median survival of approximately 1 year after multi-modal treatments [17]. Recent studies suggest that nearly 80% of all malignant brain tumors are accounted for by the broad category of gliomas, and 54% of all malignant brain tumors are GBM [18]. The first fusion protein discovered in glioblastoma was the FIG-ROS1 fusion, in which an intra-chromosomal deletion of 240 kb leads to a constitutively active kinase, suggesting oncogenic activity [19]. Two more studies reported fusions of PDGFRA-KDR [20] and LEO1-SLC12A1 [21], each in a single patient sample. The FGFR-TACC fusion is one of the recurrent fusions in GBM and it has been reported in three studies [13,22,23]. The oral administration of an FGFR inhibitor has been shown to prolong the survival of mice harboring intracranial FGFR-TACC-initiated glioma [13]. EGFR-SEPT14 fusions present in about 4% of GBMs were shown to be functional and sensitive to EGFR inhibition in a recent study [24].
In this study, we focus on identification of gene fusion events from GBM transcriptome data. Using our in-house pipeline, we identify and validate 13 fusion events in 24 GBM samples by analyzing SOLiD single-end 50 bp data. We also identify 175 high-confidence gene fusion events in 161 GBM samples by analyzing TCGA Illumina HiSeq paired-end 75 bp transcriptome data. We integrate gene fusion data with copy number data to elucidate fusion mechanisms in GBM.

Results
Gene fusion discovery pipeline for SOLiD single-end 50

bp data
We profiled the transcriptome of 24 GBM samples and 4 non-tumor samples using the SOLiD sequencer. We generated 50 bp single-end RNA-seq reads with sequencing depths ranging from 126 to 205 million reads (details provided in Additional file 1). A variety of software packages are available for gene fusion discovery for Illumina pairedend, Illumina single-end and SOLiD paired-end data [25][26][27][28][29]. We developed an in-house gene fusion discovery pipeline, as no software package was available for singleend SOLiD data (see Figure 1). First, we aligned all the reads and calculated reads per kilobase per million (RPKM) for each exon using Bioscope 1.3 software package [30]. Gene annotations were combined from three databases: Ensembl gene annotation version 66, UCSC and RefSeq genes (the tracks were downloaded on April 4th, 2012, from the UCSC genome browser [31]). Cancer Outlier Profile Analysis (COPA) [3] was performed for each exon to identify exons with substantially higher expression in a small set of samples. We evaluated the expression variation at 5′ and 3′ exons of all genes that had at least one exon with outlier expression levels. We extracted reads that partially mapped to the junction where there was a significant change in expression levels for a given gene. We then constructed a consensus sequence from the partially extracted reads. After converting the consensus sequence from color space to base space format, we used the UCSC  Figure 1 Gene fusion discovery pipeline for the Ivy Center SOLiD single-end data. Reads were aligned to the hg19 assembly using bioscope-1.3 software by Life Technologies. RPKM values were calculated for each exon, followed by modified cancer outlier profile analysis (COPA). If any of the exons of a gene displayed outlier expression in a sample, then the read distribution across that gene was evaluated for that sample. If either the 3′ or 5′ end of the gene had a considerably lower RPKM value compared to the other, the gene was further evaluated for fusion. All partially mapped sequences to a potential fusion breakpoint were extracted. One or more consensus sequence was generated and translated to base space format from the color space format. The consensus sequences were then aligned to the hg19 human genome using UCSC BLAT. If part of the consensus sequence mapped to the known exon and the other part uniquely mapped to the genome, the sequence was considered a potential fusion sequence. All potential fusion sequences were validated with fusion qPCR followed by Sanger sequencing.
BLAT tool [32] to map the consensus sequence to the human genome (hg19 assembly). If part of the consensus sequence mapped to the original gene and the rest mapped uniquely to another region in the genome, then the sequence was considered a fusion sequence. We identified 13 such sequences (see Table 1) in eight samples. We were able to validate all of the 13 fusion transcripts using fusion qPCR followed by Sanger sequencing. Figure 2 illustrates the MON2-MARS gene fusion as one example of a fusion transcript. The outlier expression of the MON2 and MARS exons can be observed with a z-score > 4. Panel A shows the RNA-seq read distribution across all exons for both genes. The MON2 read distribution shows higher 5′ expression relative to its 3′ expression, and the MARS read distribution shows higher 3′ expression relative to its 5′ expression. Partially mapped reads at exon 34 of MON2 and exon 6 of MARS map to the MON2-MARS fusion sequence. Panel B shows the gel image for fusion qPCR. The product can be observed in sample SN214 (TCGA-74-6583) but not in non-tumor brain and MON2-MARS fusion negative GBM samples. Panel C shows the Sanger sequencing trace of the fusion PCR product. Detailed images for the other 12 fusion sequences are available in Additional file 2.

Gene fusions in the Ivy Center SOLiD dataset
We identified 13 fusion transcripts in eight out of the 24 GBM samples (see Table 1). Two samples, SN214 and SN161, harbored multiple fusions. Both fusion partners in eight of the fusion transcripts are well annotated genes. these eight transcripts are predicted to be in-frame fusions coding for a chimeric protein product. Transcript CREB1 → PARD3B results in a C-terminal truncation of the 5′ fusion gene partner due to a frameshift. In transcript PLEKHA6 → PIK3C2B, the entire PIK3C2B coding sequence is preserved, but the fusion junction is at a novel 5′ UTR exon for PLEKHA6. In transcript LANCL2 → RP11-745C15, the 5′ partner gene fuses with a noncoding RNA resulting in C-terminal truncation. For the other five fusion transcripts, the 5′ partner gene fuses with genomic sequence without gene annotation, denoted as "UAR" in Table 1. Three of these transcripts (SLC35E3 → UAR, ZNF713 → UAR in two samples) result in Cterminal truncation of the 5′ partner genes. One transcript (MDM1 → UAR) is predicted to result in a shorter isoform with an extended 3′ UTR, and one transcript does not have any predicted protein product (SEC61G → UAR). We had tissue available from surgery at recurrence for patient SN159, and we were able to validate the presence of SEC61G → UAR at recurrence. Predicted protein sequences are provided in Additional file 2. All fusions are intra-chromosomal in our cohort, and fusion partners are in close proximity, ranging from a distance of 5.1 million base pairs to 35 kilo base pairs between the two partners. Although there are no recurrent fusions in our small cohort, there are multiple genes, SLC35E3, PIK3C2B and ZNF713, that occurred in more than one fusion transcript. An FGFR3 → TACC3 fusion was recently reported by three independent studies as a recurrent gene fusion [13,22,23]. All our fusion transcripts are highly overexpressed compared to their wildtype gene partners, as is evident in the third column in Table 1, which shows a much higher number of reads spanning the fusion junction compared to the number of reads spanning the known wildtype junctions. One of the fusions, SLC35E3 → UAR, has two isoforms.

Gene fusions in the Illumina HiSeq TCGA dataset
We downloaded RNA-seq data for 169 TCGA samples from CGHub [33] to explore the gene fusion landscape of GBM beyond our cohort. The TCGA transcriptome data are 75 bp paired-end reads generated using Illumina HiSeq with sequencing depths ranging from 54 to 252 million reads per sample. We used the TopHat-Fusion [25] and SnowShoes-FTD [29] software packages to identify fusions because both packages are expected to have a very low false-positive rate. There were 882 and 492 fusion sequences identified by TopHat-Fusion and SnowShoes FTD suggesting a large number of false positives (results are available in Additional files 3 and 4). The number of fusion sequences could be reduced by increasing the threshold for the minimum number of fusion spanning reads, but this modification can lead to the failure to identify some truly important fusions, such as CEP85L → ROS1. In our SOLiD dataset, the FGFR3 → TACC3 fusion has the second lowest number of junction spanning reads. Because our method resulted in a 100% validation rate, we applied filtering steps based on our method to the fusion sequences identified by both packages. We required that at least one of the breakpoints must be a known exon boundary. To reduce the likelihood of identifying passenger fusions [34], we required that at least one of the fusion spanning reads must have a ratio of greater than two compared with its corresponding wild-type exon-exon spanning reads. We discarded gene fusions involving adjacent genes. Exact details are provided in the methods section. After curating fusion sequences from both packages, we obtained a set of 175 high-confidence fusion sequences, which was referred to as the curated set. Curated fusions were present in 53% (85/161) of patients, and 22% (35/161) of these patients harbored more than one fusion. The curated fusion set is available in Additional file 5.

Gene fusions and copy number changes
The Circos plot of all curated fusions (see Figure 3) shows specific genomic hotspots where fusions occur in GBM. Two major genomic hotspots are on chromosomes 7 (7p11) and 12 (12q14-15). In our SOLiD dataset, 8 of 13 validated fusions were located on 7p11 and 12q14-15. Other regions with higher frequency of fusions are on chromosomes 1, 4, 6 and 19. These genomic hotspots for fusions are the regions that are frequently amplified in GBM, as observed in Figure 3 [35]. Because Affymetrix SNP array data were available for all but two TCGA samples, we looked for associations between fusion points and copy number data. We downloaded level 3 segmented copy number data from TCGA [36]. The start and end points of each segment were considered to be the genomic breakpoints. For the curated set, copy number data were available for 172 fusion sequences, out of which at least one of the partner genes harbored a genomic breakpoint in 135 cases (78%). We also predicted the fusion mechanism for each of the fusion sequences based on the copy number data. Figure 4 shows the distribution of different fusion mechanisms for all curated fusions. We binned fusion mechanisms into six types: 1) No copy number changes -There are no genomic breakpoints around fusion points. These could be either inter-or intra-chromosomal fusions, see Figure 4B. 2) Focal amplifications -Fusion points are within a genomic amplicon, see Figure 4C. 3) Tandem duplications -Fusion points are around the start and end of an amplified genomic segment, see Figure 4D. 4) Deletion-based -Fusion points are around the start and end of a relatively deleted genomic segment, see Figure 4E. 5) Complex genomic rearrangements -Both fusion points are around genomic breakpoints with multiple segments between the two breakpoints, see Figure 4F. 6) Inter-chromosomal -Fusion partners are located on different chromosomes with at least one fusion point near a genomic breakpoint, see Figure 4G.
Only 8% of the fusions are without accompanying copy number changes suggesting that the majority of the fusions in GBM are associated with unbalanced genomic rearrangements. Majority of the fusions in focal amplicons are present on chromosome 7 and restricted to the EGFR locus, see Figure 3 and Figure 4A. Approximately 40% of all the fusions in GBM result from complex genomic rearrangements (CGR), see Figure 4H. Some of the inter-chromosomal rearrangements also display complex fusion mechanisms, see Figure 4G. A recent study analyzed whole genome sequencing data and showed a high incidence of CGR in GBM resulting from chromothripsis-39% in GBM compared to 9% in other tumor types [37]. Fusions generated through CGRs are largely present on chromosomes 12 and 7, see Figure 4A. The distribution of CGR based fusions on chromosome 12 is largely restricted to 12q14-15 amplicon, see Figure 3. Even though partners of fusion sequences formed due to CGRs belong to different copy number segments, they have highly correlated copy number value, see Figure 4F. This suggests co-amplification of segments involving fusion.

Gene fusions and molecular features
We checked to see if samples with at least one curated fusion were enriched in any clinically associated molecular features. We did not find any association with presence of EGFR vIII, mutation/homozygous deletion of PTEN or TP53, mutation of IDH1 or G-cimp status (data obtained from CBio portal [38]). Amplifications of EGFR were more prevalent in samples with at least one fusion compared to samples with no fusions (63% vs. 38%, p = 0.0016, Fisher's exact test). We observed that the samples with classical subtype were more likely to have fusions (72%) and that samples with mesenchymal subtype were less likely to have fusions (39%), see Figure 5A. This result can be explained by the association of genomic fusion hotspots with subtypes. Almost all of the samples with a classical subtype have focal amplification of the EGFR locus, and samples with a mesenchymal subtype have a much lower incidence of focal amplifications on chromosomes 4 and 12 [35]. Figure 5B shows the chromosomal distribution of fusion breakpoints for each subtype. The majority of the fusions in the classical subtype are located on chromosomes 7, 12 and 1. Fusions in the mesenchymal subtype are mostly present on chromosomes 7, 1, 6 and 19, whereas chromosomes 12 and 4 harbor the majority of the fusions with the proneural subtype. Proneural subtype shows enrichment of fusions formed by complex genomic rearrangements, see Figure 4H. Gene fusions in samples with the neural subtype have a broader chromosomal distribution, with the majority of breakpoints on chromosomes 3, 2, 7, 6, 1, 15, 4 and 12.
Predicted structure of the curated fusions set We predicted the amino acid sequence of all curated fusions based on their chimeric nucleotide sequence. A significant portion of the fusions (37%) were predicted to be in-frame fusions with amino acid sequences present from both fusion partner genes. Another 18% were predicted to have C-terminal truncation due to either the out of frame fusion with another gene or fusion with an unannotated region. Approximately 8% of fusions are predicted to result in the same protein product as its 3′ partner gene by borrowing only the promoter from the 5′ partner. In approximately 10% of the fusions, the 5′ gene partner is predicted to contribute only the promoter, but the N-terminal of the 3′ gene is truncated.
We also observed another novel class of fusions that involve non-coding RNA genes. In approximately 14% (25/175) of the fusions, the 5′ partner gene is predicted to have a C-terminal truncation due to fusion with a noncoding RNA. These fusions also result in the expression of non-coding RNAs that are not expressed in other samples. Another important set of fusions involve tyrosine kinases. In 13 cases, the fusion sequences retained the tyrosine

Recurrent gene fusions
Although more than half of GBMs showed evidence of gene fusions, there were very few fusions that were present in more than one sample. One of those fusions is the already reported FGFR3 → TACC3 fusion, which was found in two patients in the TCGA cohort and in one patient in the Ivy Center cohort. EGFR → SEPT14, an in-frame fusion with C-terminal deletion of EGFR, was found in three patients in the TCGA cohort. LANCL2 → SEPT14, an out of frame fusion that leads to the C-terminal truncation of LANCL2, was found in two TCGA patients. Two additional patients, one TCGA and one Ivy Center, had fusions of LANCL2 with non-coding RNA RP11-745C15.2, which also resulted in the C-terminal truncation of LANCL2. The same non-coding RNA RP11-745C15.2 fused with EGFR in two TCGA patients, resulting in the C-terminal truncation of EGFR. There are 27 genes that are fusion partners in more than one patient sample (see Table 2). The majority of these genes (18/27) are on genomic fusion hotspots located on chromosomes 7 and 12.

Discussion
Our study highlights the prevalence of gene fusions as one of the major genomic abnormalities in GBM. Fusions occur in approximately 30-50% of GBM patient samples.
In the Ivy Center cohort of 24 patients, 33% of samples harbored fusions that were validated by qPCR and Sanger sequencing. We were able to identify high-confidence gene fusions from RNA-seq data in 53% of samples in a TCGA cohort of 161 patients. We identified 13 cases (8%) with fusions retaining the tyrosine kinase domain in the TCGA cohort and one case in the Ivy Center cohort. Recent advances in the development of tyrosine kinase inhibitors (TKIs) have demonstrated that these drugs can provide significant benefit to patients whose tumors have a specific genetic abnormality. We also identified a novel class of fusions (14%) that result in the C-terminal truncation of its 5′ partner due to fusion with non-coding RNA genes. One such case was also present in the Ivy Center cohort. This study reveals the diversity of gene fusions in GBM samples. The majority of the fusions are private fusions occurring in one patient. There are a few fusions that recur at low frequency in GBM.
Our study is the first to provide a comprehensive view of the gene fusion landscape in GBM by examining sequences from 185 patients from two independent cohorts. We successfully utilized our in-house pipeline for fusion discovery using SOLiD single-end, 50 bp RNA-seq data with a 100% validation rate. For the TCGA cohort, we used two different gene fusion detection software packages to comprehensively identify fusions from Illumina paired-end, 75 bp RNA-seq data. Ours is the first study to describe recurrent fusions involving non-coding genes. We combined copy number data with gene fusion discovery to elucidate mechanisms of the formation of gene fusions in GBM. All of the fusions detected in this study can be further visualized and analyzed on our website (http://ivygap.swedish.org/fusions).
We were able to validate all of the fusions in our SOLiD single-end RNA-seq data by using strict filtering criteria.  It is likely that we may have underestimated fusions for Ivy Center data. Due to lack of access to the tissue samples, we could not determine the validation rate for our set of curated fusions in the TCGA cohort. The curated fusion set did have a significantly higher percentage of fusions associated with copy number changes relative to the low-confidence set. We applied filters to discard likely passenger fusions [34,39], but the functional significance of these fusions still needs to be evaluated. Singh et al. was the first study to describe multiple fusions of FGFR-TACC in GBM, reporting this phenomenon in 3 of the 97 tumors examined. They showed that the fusion protein has oncogenic activity when introduced into astrocytes and oral administration of an FGFR inhibitor prolongs the survival of mice harboring intracranial FGFR-TACC-initiated glioma [13]. A second study by Parker et al. showed that the fusion gene is overexpressed by escaping miR-99a regulation due to loss of the 3′ UTR of FGFR3 [22]. In their cohort, 4 out of 48 samples harbored the FGFR3 → TACC3 fusion. In our Ivy Center cohort, the FGFR3 → TACC3 fusion was detected in one out of 72 samples. We tested for this fusion in an additional 48 samples in addition to the 24 RNA-seq samples, but did not detect any fusion events. In the TCGA cohort, 2 of 161 samples harbored the FGFR3 → TACC3 fusion. Fusions of FGFR genes are identified in other cancers, including bladder cancer, cholangiocarcinoma, squamous lung cancer, breast cancer, thyroid cancer, oral cancer, head and neck squamous cell carcinoma and prostate cancer [15,23]. Tropomyosin-Receptor Kinases (Trk) are known to play a role in cancer biology. Rearrangements of the NTRK1 gene are consistently observed in a small fraction of papillary thyroid carcinomas [40].
We identified two cases of NTRK1 fusions in the TCGA cohort. Frattini et al. [24] screened 248 samples for NFASC-NTRK1 fusion but did not find any. We identified a single case of a CEP85L-ROS1 fusion in the TCGA patient samples. In a recent study by Giacomini et al. [41], a CEP85L-ROS1 fusion was detected for angiosarcoma. There have been two more reported cases, one angiosarcoma and one epithelioid hemangioendothelioma, with ROS1 rearrangements. ROS1 rearrangements also define a unique molecular subclass of lung cancer that may respond to an ALK inhibitor [42]. We identified fusions of EGFR in nine patient samples from the TCGA cohort, out of which six retained the tyrosine kinase domain and resulted in a carboxyl-terminal truncation. A study by Cho et al. has shown that cetuximab prolonged the survival of intracranially xenografted mice with oncogenic EGFR carboxyl-terminal deletion mutants compared with untreated control mice [43]. It is likely that patients with fusions of EGFR leading to carboxyl-terminal truncation will show sensitivity to EGFR inhibitors. Frattini et al. [24] showed that EGFR-SEPT14 fusions which occur in about 4% of GBMs was a functional gene fusion in GBM and confers mitogen independence and sensitivity to EGFR inhibition. A total of 13 cases from both cohorts have fusions of genes involved in chromatin remodeling and modification. These genes include ARID1A, ARID1B, ASH1L, CHD4, HDAC1, HMGA2, JMJD1C, KDM4B, RERE, SETD1B and YEATS4. ARID1A-MAST2 fusion has been shown to be a critical driver fusion in an MDA-MB-468 breast cancer cell line [10]. In 27 samples, the 5′ partner gene fuses with non-coding RNA. These fusions are predicted to have a C-terminal truncation. These cases also have highly expressed non-coding RNAs Table 2 Fusions involving genes that partner in more than one fusion (Continued)   TCGA-28-5209-01A  LANCL2  chr7  55433921  RP11-310H4  chr7  55714590   TCGA-14-1829-01A  SEC61G  chr7  54823471  RP11-310H4  chr7  There is now growing evidence of an oncogenic and tumor suppressive role for long, non-coding RNAs in tumor biology [45]. Their identification in gene fusion events has thus far been neglected, as most studies focus on fusions of the coding genes. Even though gene fusion events in GBM are abundant with scarce recurrent events, they are not random events. Majority of the fusion events occur at 7p11, 12q14-15, 1q32 and 4q12 which are also recurrently amplified regions in GBM. These fusion hotspots are consistent in both Ivy and TCGA cohorts. Also majority of the fusion events are due to unbalanced genomic rearrangements. Analysis of whole genome sequencing data also showed that 88% of genic rearrangements in GBM are associated with copy number alterations [46]. Some of the key genes implicated in GBM biology within these hotspots are EGFR, MDM2, CDK4, PIK3C2B, MDM4 and PDGFRA. A recent study [46] identified a dense breakpoint pattern on 12q14-15 indicative of local chromosome instability and defined this region as "breakpoint enriched region" (BER). They showed that patients with BER pattern had poor survival and this pattern was associated with MDM2/CDK4 coamplification. There are two other cancers, dedifferentiated liposarcomas and lung adenocarcinomas that also show MDM2/CDK4 co-amplification in 90% and 4% of cases respectively [38,47]. All three types of cancer display distinct genomic aberration patterns in 12q14-15 region in spite of having MDM2/CDK4 co-amplification. GBM samples show shattering of the region with alternate high level deletions and gains, lung adenocarcinomas mostly contain large amplified segments and dedifferentiated liposarcomas contain multiple amplified segments (see Additional file 6). Whole genome sequencing, copy number and RNA-seq datasets show that GBMs contain deletion bridges that connect these amplified segments and generate a large number of fusion transcripts. Such complex genomic rearrangements are more prevalent on chromosome 12 but not limited to as shown in the study by Malhotra et al. [37] where they analyzed whole genome sequencing data of 18 GBMs. About 40% of fusion transcripts are formed due to such complex genomic rearrangements. With the advent of RNA-seq technology the list of fusion sequences in solid tumors is growing exponentially but little is known about the mechanisms that facilitate fusion events. The formation of the TMPRSS2-ERG gene fusion that occurs in about 50% of prostate cancers has been shown to be facilitated by androgen signaling which induces proximity of the TMPRSS2 and ERG genomic loci and then exposure to gamma irradiation which causes DNA double-strand breaks [48]. The overview of the fusion landscape in GBM leads to questions about what mechanisms are responsible for generating highly site specific DNA double-strand breaks and then joining of these breaks that result in complex genomic rearrangements.

Conclusions
Gene fusions are frequent genomic abnormalities in GBM. The majority of the fusions are private fusions, with a minority recurring in multiple patients. Complex genomic rearrangements are the major mechanism by which fusions are formed in GBM. Due to the low frequency and rarity of clinically relevant fusions, RNA-seq of GBM patient samples is an essential tool for the identification of patient specific fusions that can drive personalized therapy.

Ethics statement
This study was reviewed and approved by Western IRB (IRB00000533) in compliance with the ethical principles set forth in the report of the National Commission for the Protection of Human Subjects of Biomedical and Behavioral Research, titled "Ethical Principles and Guidelines for the Protection of Human Subjects of Research (Belmont Report)". The research protocol was also approved by the Swedish Neuroscience Institute research steering committee. All participants provided written informed consent according to IRB guidelines prior to their participation in this study.

Patient samples
Tumors were obtained from surgeries performed in the years 2009 through 2011 at the Swedish Medical Center (Seattle, WA) according to institutional guidelines. Patient samples used in this study had a histopathology diagnosis of WHO grade IV glioblastoma multiforme.
Transcriptome sequencing on SOLiD 5500 RNA isolation and purification Total RNA was extracted from human brain tumor tissues with Trizol (Life Technologies, CA) and then purified using the MEGAclear kit (Life Technologies) as per the manufacturer's instructions. The integrity and quantity of RNA was assessed on the Agilent 2100 Bioanalyzer (Agilent, CA) as per the manufacturer's recommendations.
Ribosomal RNA depletion from total RNA Qualified total RNA was subjected to depletion of ribosomal RNA by using the Ribo-Zero rRNA removal Kit (Epicentre, IL). A total of 5 μg of purified total RNA was mixed with rRNA removal reagents for 25 minutes, added to prepared Ribo-Zero microspheres according to the manufacturer's instructions, and then incubated for 20 minutes. The mixture was applied to a spin-filter column and centrifuged for 2 minutes to remove the microspheres. rRNA-depleted total RNA was concentrated using the Ribominus concentration module (Life Technologies, CA) and assessed on the Agilent 2100 Bioanalyzer for the confirmation of rRNA removal.

RNA fragmentation
A total of 500 ng of rRNA-depleted total RNA was subjected to fragmentation by chemical hydrolysis using the SOLiD Total RNA-Seq kit (Life Technologies, CA) according to the manufacturer's instruction and was assessed on the Agilent 2100 Bioanalyzer for fragment yield and size distribution.

Construction of an amplified whole transcriptome library
The fragmented rRNA-depleted total RNA samples were used for the construction of an amplified library using the SOLiD Total RNA-Seq kit (Life Technologies, CA) according to the manufacturer's instruction. Briefly, 100 ng of fragmented RNA was hybridized with SOLiD adaptor mix and followed by ligation of the fragments. Reverse transcription was performed with SOLiD RT primers to generate the cDNA library. The cDNA library was then purified and size selected using AMPure XP reagent (Agencourt, CA) as per the manufacturer's instruction. Amplification of the cDNA library was performed for multiplex SOLiD sequencing using barcoded 3′ primers. Purification of the amplified DNA was performed using the PureLink PCR micro kit (Life Technologies, CA). Purified DNA was assessed on the Agilent Bioanalyzer 2100 for yield and size distribution.

Sequencing
The bar-coded libraries were quantified by using the SOLiD Library TaqMan Quantitation kit (Life Technologies, CA), and four bar-coded libraries were pooled together in equal concentrations into one pool. The pooled libraries were used as the template for the next step of emulsion PCR and were followed by enrichment. Emulsion PCR and enrichment were performed at the E120 scale in SOLiD EZ Bead System (Life Technologies, CA) according to the manufacturer's instructions. Each pool was sequenced in a SOLiD FlowChip on the SOLiD 5500 (Life Technologies, CA) according to the manufacturer's instructions.
TCGA transcriptome and copy number data TCGA transcriptome data were downloaded from CGHub [33]. The level 3 copy number data were obtained from the TCGA data portal [36].
Gene fusion discovery process for SOLiD 5500 data Reads were aligned to hg19 assembly using bioscope 1.3 software by Life Technologies [30]. The gene annotation file was obtained by combining annotations from Ensembl gene annotation version 66, UCSC and RefSeq genes (the tracks were downloaded on April 4th, 2012, from the UCSC genome browser [31]). RPKM values were calculated for each exon, followed by a modified cancer outlier profile analysis (COPA) [3]. If any of the exons of a gene displayed outlier expression in a sample, then the read distribution across that gene was evaluated for that sample. If either the 3′ or 5′ end of the gene had a considerably lower RPKM value compared to the other end, the gene was further evaluated for fusion events.
Alignments Bioscope 1.3 was run using its default settings. The RPKM values were calculated using the "Count Known Exons" tool with quality cutoffs minMapq = 10 and scoreClear-Zone = 5. Exons that have an RPKM value greater than 20 in at least one of the samples were evaluated for outlier expression.

Cancer outlier profile analysis (COPA)
For each exon, RPKM values are sorted in ascending order. We calculate z-score z i in sample i as Z i = (x i − μ)/σ where average and standard deviation are calculated as follows: where n = number of samples and k = index to the array of sorted RPKM values. An exon is considered to have an outlier expression if the z-score is greater than 4.

Exon-walking RNA-seq expression pattern
Earlier studies have utilized exon-walk PCR to identify fusion breakpoints [3]. We used a similar approach using RNA-seq RPKM data for each exon. For each exon number j of gene i and sample k RPKM is normalized by the 7th quantile RPKM values for exon number j of gene i as follows: When walking from j th exon to the (j + 1) th exon, if there is a two-fold drop or rise in normalized RPKM value, then the exon-exon boundary is considered a potential fusion breakpoint.

Consensus sequence
All partially mapped sequences to a potential fusion breakpoint were extracted. These sequences have less than 35 matches to the known exon. One or more consensus sequences were generated and translated to base space format from the color space format. At least two sequences were used to generate a consensus sequence.

Blat
The consensus sequences were then aligned to the hg19 human genome using BLAT [32].
If the part of the consensus sequence mapped to the known exon and the other part uniquely mapped to the genome, the sequence was considered a fusion sequence.
Fusion qPCR cDNAs were synthesized by using the High Capacity cDNA Reverse Transcription Kit (Life Technologies) with 1 μg of purified total RNA. Primers specific for fusion genes that were used in RT-qPCR are listed in Additional file 7. GUSB was used as the internal reference gene. For each fusion sequence three samples were used: the GBM sample containing the fusion, the GBM sample without that fusion and the non-tumor brain sample.

Sanger sequencing
The RT-PCR products were selectively extracted from an agarose gel and cloned into the pCR2.1-TOPO cloning vector (Life Technologies). All clones were confirmed by sequencing using 3130 Genetic Analyzer (Life Technologies).

Curated fusions
For TopHat fusions, we considered all of the potential fusions in the output file, potential_fusion.txt, and not just the fusions reported in result.txt. We used the output file final_fusion_report_RNA.txt for SnowShoes-Ftd fusions. For all exons of the genes involved in fusion transcripts, we calculated z-scores as described in the above section. Exon RPKM for the TCGA data was calculated using script coverageBed in package BEDTools-Version-2.16.2 [49]. The following conditions were met by the fusion transcripts in the curated set: