Illumina Genome Analyzer sequencing

Each prepared cDNA library was subjected to high-throughput sequencing using the Illumina Genome Analyzer IIx platform, following the standard RNA sequencing protocol. Samples were sequenced at various sites and times and are described as follows: HCC1937 was sequenced by Illumina (50 bp SE). SUM149PT, SEC1 and SEC2 were sequenced at the Institute of Cancer Research (36 bp, 54 bp and 54 bp PE, respectively). SUM1315O2, HCC3153, HCC2337, T92, T50, T160 and MCF10A were sequenced at the McGill University and Genome Quebec Innovation Centre (76 bp PE per sample). All reads were processed using ELAND (AJ Cox, unpublished) to generate FASTQ files. Due to poor base quality at the 3' end of SUM1315O2 and HCC3153, reads in these samples were trimmed by 40 bp from the 3' end to produce 36 bp high quality PE reads. Additional read statistics are summarized in Additional file 1.

RNA-Seq post-processing

Illumina FASTQ reads were first converted to Sanger FASTQ format using the ill2sanger script (patch by D Cittaro) from MAQ (http://maq.sourceforge.net/). Reads were then aligned to the human reference genome (hg18) using BWA [25]. Reads that span an exon-exon splice junction will not map to the reference since it will be separated by an intron. Thus, the remainder of the unmapped reads were aligned to a custom library of junction sequences based on known RefSeq annotation. For each splice junction and read length x, we joined x-5 bp (or the length of the exon if less than x-5) of sequence from the upstream exon with x-5 bp (or the length of the exon if less than x-5) of the downstream exon. This ensured that reads mapping to a splice junction overlapped to at least 5 bp of either end of the junction. The results were then processed using SAMTools [26] for downstream analysis. Visualization of the mapped reads was carried out using Integrated Genome Viewer (http://www.broadinstitute.org/igv).

Identification of gene fusion transcripts from RNA-Seq

In this study, we implemented two strategies adapted from methods described by Maher et al. [3,4] to facilitate gene fusion discovery with SE reads and PE reads, respectively.

In the SE approach, gene fusions are identified by checking for reads which map across the exon-exon fusion junction between the two fused genes. To simplify the analysis, we reasoned that only the subset of reads which did not fully align to the reference genome or splice junctions would be interesting (unmapped reads), as they may harbour fusion junction spanning reads. All unmapped reads were realigned to RefSeq mRNA sequences using Blat [27]. Reads that partially aligned to a RefSeq sequence (20-70% identity) and did not fully map to anything else were retained. Ambiguous reads that were partially aligned to more than five genes were also discarded. Based on the position of the alignment, we further filtered for reads that partially aligned within 5 bp of an exon boundary using RefSeq annotation. This is because we expected that potential reads mapping across a gene fusion junction would be aligning at the exon boundaries. The resulting set of partially aligning reads was denoted as group A. Next, the unaligned portions of these reads were extracted and again subjected to alignment against RefSeq mRNA exon sequences using Vmatch [28], for effective alignment of shorter reads (<20 bp). Reads that fully mapped to an exon and localized within 5 bp of 5' end or 3' were retained in a similar fashion as before (denote as group B). Finally, we combined the results from group A with group B by finding reads that had a partial alignment in group A and a full alignment in group B, thus accounting for the entire read sequence and ultimately representing a gene fusion candidate read (GFCR). All the GFCRs were tallied and summarized. Fusion partners with greater than or equal to three unique GFCRs were considered a gene fusion candidate. Each candidate was manually examined to determine its potential biological significance and to check whether it could be a false positive. This included disqualifying reads that mapped to too many loci, which may represent a repetitive region, or gene fusion candidates which mapped to two homologs that were part of the same gene family.

For the PE approach, we searched for two features which may define a gene fusion. First, we limited our selection to reads that fully and uniquely mapped to the reference genome. Second, we reasoned that reads which map to different chromosomes may reflect an interchromosomal translocation and have the potential to form a gene fusion. Alternatively, reads that map to the same chromosome but at a larger than expected distance apart may represent an intrachromosomal translocation. However, this may be confounded by large introns situated between two exons that the reads map to. To account for this, we used an arbitrary minimum intrachomosomal mapping distance of 1 Mb. After identifying such read pairs, we realigned the reads to RefSeq mRNA sequences using Blat. Ensuring that each read in a pair mapped to a different gene, we summarized a list of candidate gene fusions in a similar fashion as for SE reads and generated a paired end fusion score (PEFS). The scoring scheme was defined as follows: +2 for each identified read pair, -1 for each duplicate read, and -0.5 for each mismatch in a read. We prioritized those that had a minimum of three read pairs that flanked a fusion junction and had a minimum PEFS score of 5. For further supporting evidence of each finding, we searched for additional reads that spanned across the fusion junction, analogous to the SE approach but restricted to a search space containing only the two partner genes. Here, we generated a library of all possible splice junction sequences between the 5' fusion partner and the 3' fusion partner, which includes the predicted fusion junction. All remaining unmapped reads were aligned to these junction sequences using Blat.

Quantification of exon expression levels

Exon expression levels were quantified by first determining the number of reads that uniquely mapped to each exonic region using BEDTools [29]. The genomic coordinates for each exon were obtained from the UCSC Table Browser [30]. The number of reads was normalized and expressed as reads per kilobase of exon model per million mapped reads (RPKM) [31].

Sanger sequencing

Primer pairs for PCR amplification and sequencing of each coding exon were generated using Primer3 [32]. Primers (see Additional file 2 for sequences) were designed to span the predicted exons forming the fusion breakpoint and to generate a maximum product size of 300 bp, which was considered optimal for amplification, purification, and sequencing. To minimize amplification of homologous genomic sequences, primer pairs were filtered using the UCSC In Silico PCR software, and only pairs yielding a single product were used. PCR reactions were performed on the cell lines' cDNA at least twice, using the following thermocycling parameters: 95°C × 15 min, (95°C × 20 s, 60°C × 20 s, 72°C × 20 s) for 30 cycles, 72°C × 10 min. PCR products were purified as recommended by the manufacturer (QIAGEN). Products were sequenced by conventional Sanger methods and compared to the reference sequence to validate the gene fusion. Sequence products were obtained from MUGQIC. Sequence chromatograms were aligned and analyzed with the Staden package and Mutation surveyor software version 3.24.

RT-PCR

Total RNA (100 ng) was reverse transcribed into single-stranded cDNAs using SuperScript III reverse transcriptase (Invitrogen) and Oligo (dT)_12-18 (Invitrogen) in 20 μl reaction at 50°C for 50 min, 85°C for 5 min, 37°C for 20 min. 2 μl of cDNA was used for a subsequent 20 μl PCR amplification. To detect fusion transcripts, we design the forward primer targeting the 5' partner gene and reverse primer targeting the 3' partner. Primer pairs (see Additional file 2 for sequences) for the coding exons of the fusion genes were generated using Primer3 [32]. GAPDH was amplified simultaneously. PCR was performed at an annealing temperature of 60°C for 26 cycles. PCR products were separated by gel electrophoresis in a 3% agarose gel (Metaphor, MRC) and visualized by ethidium bromide staining.

Proof of concept

To evaluate our implementation, we attempted to identify known gene fusions within our samples. We applied our SE approach on HCC1937. This cell line has been previously found to have a translocation involving exon 2 of NFIA on chromosome 1 joining with exon 5 of EHF on chromosome 11 using massively parallel DNA sequencing [24]. We were able to successfully confirm the expression of NFIA-EHF in the transcriptome of this cell line. Several 50 bp reads were mapped across the junction between exon 2 of NFIA and exon 5 of EHF (Figure ​(Figure1A1A).

Open in a separate window

Figure 1

RNA-Seq evidence of previously described gene fusions. We first tested our SE approach on (A) the HCC1937 cell line that harbours the fusion NFIA-EHF. SE reads are shown to map across the exon-exon fusion junction between exon 2 of NFIA and exon 5 of EHF, as illustrated in the schematic. Next, we tested our PE approach on (B) two primary tumors that contain ETV6-NTRK3. Results from the sample SEC2 are shown. Paired reads (indicated by two solid lines joined by a dotted line) as well as single reads (red lines) are shown to map across the fusion junction between exon 5 of ETV6 and exon 14 of NTRK3.

Next, we tested the PE approach on two primary tumors (SEC1 and SEC2) with secretory breast ductal carcinoma, a rare subtype of breast cancer. These samples were tested by other molecular methods to harbour the gene fusion ETV6-NTRK3 [34]. This fusion is caused by a balanced translocation, t(12;15)(p13;q25), joining exon 5 of ETV6 and exon 14 of NTRK3 and is recurrent in this cancer [18]. Furthermore, functional studies have established the fusion to have a causative role in the pathogenesis of this breast cancer subtype [35,36]. We successfully confirmed the gene fusion transcript in SEC2 by identifying read pairs that flanked the fusion junction as well as individual reads mapping across the fusion junction (Figure ​(Figure1B).1B). The reciprocal fusion, involving exon 13 of NTRK3 as the 5' partner gene and exon 6 of ETV6 as the 3' gene, was also identified in the same sample. On the other hand, we could not find the fusion in SEC1. Since no supporting PE reads were found, we checked to see if there were only individual reads that mapped across the ETV6-NTRK3 fusion junction. To test this, we mapped our reads to a custom-generated ETV6-NTRK3 junction sequence using Blat [27] and did not find any supporting reads. Subsequent inspection of starting sample from which SEC1 was prepared indicated that the tumor content was less than 50%. Hence, there may have been insufficient content of this gene fusion to be detected by RNA-Seq.

Discovery of gene fusions

We analyzed the remaining samples for novel gene fusions. From the initial list of candidates, we found that many of the candidates with fewer than three supporting paired reads were due to duplicated reads, which may represent false positives. Samples T50 and T160, in particular, were sequenced at a higher depth than the other samples, yielding more sequencing data. A summary of the number of fusions found in each sample is shown in Table ​Table22 In total, we found four candidates in samples with BRCA1 mutations that were supported by both misaligning PE reads and individual fusion junction-spanning reads. No fusions were identified in MCF10A. We manually inspected each of the four candidates and found that three of them formed in-frame transcripts, as described below and summarized in Table ​Table33.

We identified a novel in-frame interchromosomal fusion transcript in SUM149PT involving exon 6 of MTAP on chromosome 9 joining with exon 3 of PCDH7 on chromosome 4. We found five read pairs flanking the fusion junction and two individual reads spanning the junction (Additional file 3A). The 5' gene, MTAP, encodes for a methylthioadenosine phosphorylase that has been frequently observed to be co-deleted with tumor suppressor gene encoding p16 in numerous cancers [37]. Interestingly, a deletion resulting in a fusion protein between MTAP and tumor suppressor gene encoding p15^INK5B was reported in a glioma xenograft as well as other malignant cell lines [38]. The 3' gene, PCHD7, encodes for an extracellular protocadherin protein involved in cell-cell recognition and adhesion. We were unable to amplify the MTAP-PCDH7 fusion junction using RT-PCR of cDNA in SUM149PT. One possible explanation for this failure could be due the low expression of both genes in this sample. We investigated this by examining their gene expression profiles in the RNA-Seq data, as will be discussed later.

A second in-frame interchromosomal fusion transcript was found in HCC3153 involving exon 19 of WWC1 on chromosome 5 joining with exon 10 of ADRBK2 on chromosome 22. We found three read pairs flanking the fusion junction and seven individual reads spanning the junction (Figure ​(Figure2A).2A). The 5' gene, WWC1, encodes for a protein KIBRA, a cytoplasmic phosphoprotein that regulates the Hippo/SWH signaling pathway and has been shown to be involved in tumor suppression [39]. The 3' gene, ADRBK2, encodes for a beta-adrenergic receptor kinase involved in the phosphorylation of G protein-coupled receptors. The ATP binding and kinase domains of this gene are conserved in the predicted fusion. To verify our finding, we performed PCR amplification of the WWC1-ADRBK2 cDNA at the fusion junction followed by Sanger sequencing (Figure ​(Figure2B).2B). Furthermore, in a separate experiment not described in this study, HCC3153 was subsequently re-sequenced at the Institute of Cancer Research using RNA-Seq on the Illumina Genome Analyzer. Our analysis of this data re-confirmed the expression of the predicted fusion transcript (data not shown). To determine for recurrence, we additionally screened a cDNA panel of 57 breast cancers and did not find the fusion in any of the cases.

Open in a separate window

Figure 2

RNA-Seq and Sanger sequencing of WWC1-ADRBK2. We identified an in-frame gene fusion transcript in HCC1315: (A) schematic of the predicted gene fusion illustrating RNA-Seq evidence that support the fusion between exon 19 of WWC1 and exon 10 of ADRBK2. Reads are indicated by black solid lines. Paired reads are indicated by the dotted line joining two reads. Reads that span across the fusion junction are highlighted by red solid lines; and (B) the fusion junction was verified using Sanger sequencing of cDNA.

A third in-frame intrachromosomal fusion transcript was identified in a primary tumor T50 involving exon 2 of ADNP and exon 17 of C20orf132, both of which are located on chromosome 20 (Figure ​(Figure3A).3A). Interestingly, while identifying additional reads that mapped across the fusion junction, we found read evidence supporting an additional fusion isoform of ADNP-C20orf132 involving exon 1 of ADNP instead of exon 2 (Figure ​(Figure3C).3C). However, no PE reads were found to be discordantly mapped across this fusion isoform. Notably, exon 2 is known to be alternatively spliced in ADNP, resulting in a long and short wild-type transcript isoforms. Both of these two exons are located in the 5' UTR of ADNP, suggesting that this gene fusion may encode for a truncated version of C20orf132. ADNP encodes for an activity-independent neuroprotector homeobox protein that has been found to be involved in cell survival, due to its proximity to a region that is frequently amplified in cancer, 20q12-13.2 [40]. On the other hand, C20orf132 is an uncharacterized protein and no functional annotation was available. We validated both fusion transcript isoforms using Sanger sequencing (Figure ​(Figure3B3B and ​and3D)3D) and RT-PCR (Additional file 4) of cDNA from the primary tumor. Furthermore, we were not able to find evidence of this gene fusion expressed in the panel of 57 breast cancer samples discussed above, suggesting that it is not recurrent.

Open in a separate window

Figure 3

RNA-Seq and Sanger sequencing of ADNP-C20orf132. We identified another in-frame gene fusion in a primary tumor that was present as two transcript isoforms: (A) schematic of the first predicted gene fusion isoform illustrating RNA-Seq evidence that support the fusion between exon 2 of ADNP and exon 17 of C20orf132; (B) Sanger sequencing of the fusion junction of the first isoform; (C) while searching for fusion junction-spanning reads, we subsequently identified a second isoform in which exon 1 of ADNP was fused with exon 17 of C20orf132; and (D) Sanger sequencing of the fusion junction of the second isoform. Reads are indicated by black solid lines. Paired reads are indicated by the dotted line joining two reads. Reads that span across the fusion junction are highlighted by red solid lines.

Lastly, by sequencing HCC1937 and HCC2337, a tumor cell line and matched lymphocyte cell line derived from the same patient, respectively, we sought to determine whether NFIA-EHF was a germline mutation or an acquired somatic mutation. Our analysis could not detect the expression of this gene fusion in HCC2337, suggesting that it represents the latter. As an additional observation, we confirmed via visual inspection of the mapped reads the presence of three mutations originally reported by Tomlinson et al. [41]: 1) a 5266dupC germline mutation in BRCA1; 2) a somatic deletion of PTEN; and 3) a somatic C->T substitution in TP53 (data not shown).

Expression profiling of candidate gene fusions

We next studied the expression profiles of the candidate gene fusions based on our RNA-Seq data (Methods). The expression levels of every exon of each gene was quantified using a normalized metric described by Mortazavi et al. [31] called reads per kilobase of exon model per million mapped reads (RPKM). For each exon e, we examined the expression fold change (FC) ratio between the sample that the fusion gene is present in (denote as e_S) and the average of our other fusion-negative samples (denote as e_R). Formally, we define FC = RPKM(e_S)/RPKM(e_R).

Our analysis of MTAP-PCHD7 showed that both genes were not highly expressed in SUM149PT compared to the other samples. This was illustrated by plotting the log₂-transformed FCs of all exons in each gene, as shown in Additional file 3B. Interestingly, for WWC1-ADRBK2, we observed a discordance in expression delineated at the predicted fusion breakpoint of both partner genes (Figure 4A-B). Furthermore, the retained 3' part of ADRBK2 showed increased expression over the other samples, which may be indicative of transcriptional activation driven by a foreign promoter or enhancer in WWC1. To experimentally confirm that the wild-type transcripts of these two genes were not expressed, we designed primers using RT-PCR (Figure ​(Figure4C).4C). As expected, only the fusion fragment was detected in HCC3153 (Figure ​(Figure4D).4D). Lastly, we repeated the exon expression analysis for ADNP-C20orf132 and similarly found discordant expression in C20orf132. Here, the exons at the 3' end of the C20orf132 (exons 17-24) were expressed higher than the 5' end (exons 1-16) (Additional file 4).

Open in a separate window

Figure 4

Expression profile of WWC1 and ADRBK2. We studied the gene expression profiles of WWC1 and ADRBK2: (A) Expression plot of WWC1 exons as measured by the log₂ fold change (FC) between the RPKM values of each exon in HCC3153 versus the average of all other WWC1-ADRBK2-negative samples. The predicted gene fusion junction is marked by an arrow. It can be seen that exons downstream of exon 19 in WWC1 is underrepresented in HCC3153. (B) Similarly, the same can be observed in the exons upstream of exon 10 in ADRBK2. Moreover, the 3' end of the gene appears to be overexpressed. (C) Two sets of primers were designed to test for expression. The first set, F₁ and R₁, spanned exons 9 and 10 of ADRBK2 to test for wild-type expression. Similarly, the second set, F₂ and R₂, spanned exons 19 and 20 of WWC1. Both sets were also combined to test for expression of predicted WWC1-ADRBK2 fusion. The expected fragment sizes are shown in brackets. (D) Exon-exon RT-PCR results: a normal control was tested for the wild-type expression of ADRBK2 by F₁/R₁ (lane 1), the wild-type expression of WWC1 by F₂/R₂ (2), and with both sets together (3). In lane 4, both primer sets were applied on HCC3153 and confirmed expression of the expected fusion fragment but no wild-type expression.