Global analysis of somatic structural genomic alterations and their impact on gene expression in diverse human cancers

CNA-Based Evaluation of SV Detection Methodology.

We first noted large discrepancies, for example in terms of SV types, when applying several available tools for SV detection to WGS data from TCGA tumor/normal pairs (Fig. 1A), and therefore sought to carefully evaluate multiple options. A general problem in assessing SV results is the lack of true positive data for comparison (19). To partially overcome this problem, we made use of orthogonal results in the form of CNA profiles from the Affymetrix SNP6 microarray platform (SI Methods). CNAs arise as a consequence of somatic SVs in the genome, and segmented CNA data therefore provide a true positive subset. Although this subset is incomplete, an ideal experimental and computational pipeline should essentially be able to detect all CNA breakpoints.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Evaluating SV detection tools. (A) Four SV detection tools were applied to 11 WGS tumor/normal pairs. The distribution of SV classes is shown for each tool. SVDetect identifies two classes of SVs: inter- and intrachromosomal translocations. The “other” class for Meerkat represents tandem duplications. (B) The fraction of CNA breakpoints overlapping with SV breakpoints was used as a proxy to determine the sensitivity of each tool (red bars). Randomized CNA breakpoints were used to assess the specificity of each tool (blue bars).

We applied four available tools, BreakDancer (20), SVDetect (21), DELLY (22), and Meerkat (11), to high-coverage WGS data from 200 TCGA tumors of diverse origins, and quantified overlaps between SV and CNA breakpoints (<15-kb tolerance). To reduce the computational complexity, only a single chromosome (chr) was considered in this analysis (chr8). To also address the specificity of the methods, we investigated overlaps with a set of randomized CNA breakpoints that maintained a similar positional distribution across the genome (SI Methods). SVDetect had the highest sensitivity score (35%), followed by Meerkat (29%) (Fig. 1B). However, Meerkat performed considerably better than SVDetect in terms of specificity (Fig. 1B). Additionally, Meerkat performed favorably compared with the other methods when applied to 11 complete tumor genomes of diverse origins (Fig. S1), and therefore we chose Meerkat for SV detection.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S1.

Evaluating SV detection tools. SV detection tools were evaluated, shown on the x axis. The fraction of CNA breakpoints overlapping with SV breakpoints was used as a proxy to determine the sensitivity of each tool (red bars). Randomized CNA breakpoints were used to assess the specificity of each tool (blue bars). Meerkat had the highest sensitivity score (38%), followed by SVDetect (28%), DELLY (27%), and BreakDancer (26%).

A Map of Somatic SVs Across 600 Tumor Genomes.

We next used high-coverage WGS data from TCGA to map somatic SVs in 600 tumors across 18 cancer types (Dataset S1), and combined these results with somatic CNA (Affymetrix SNP6 microarrays) and expression (RNA-seq) data from the same samples. In total, 51,446 somatic SVs were identified (on average 87 events per sample). Of these, 36,382 were intrachromosomal (13,745 deletions, 13,335 tandem duplications, and 9,302 inversions), whereas 14,953 were interchromosomal translocations. The number of SVs and the relative distribution between SV classes varied considerably within and across cancer types (Fig. 2A).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Overview of detected somatic SVs and overlaps with CNAs across 18 cancer types. (A) Samples are grouped by cancer type and ordered by the number of detected SVs. (B) The fraction of CNA breakpoints (Affymetrix SNP6 data) detected by the SV detection pipeline is shown for each tumor. Red lines indicate the median in each cancer type. Cancer types are ordered by the degree of copy number instability (average CNA segments). Segments involving breakpoints in centromeres and telomeres (within 2 Mb) were filtered out for this analysis, as these may not be assessable by SV callers. (C) Scatter plot showing correlation between the number of CNA and SV breakpoints in each sample. “Mostly arm level” indicates that more than half of all chromosome arms have arm-level CNAs. “Chromothripsis effected” indicates at least one chromosome arm with more than 30 copy number segments. (D and E) Two examples of complex relationships between CNAs and SVs. Blue represents deletions, and red represents amplifications. The linear copy number amplitude is indicated for the CNA data. (F) Structural basis of copy number amplifications and deletions. Only CNA segments where both breakpoints could be matched to those of a single SV event were considered (15-kb tolerance). Pie chart slices are color coded similar to Fig. 2A.

Confirming previous data (10, 23⇓–25), breast (BRCA), ovarian (OV), and lung squamous (LUSC) tumors had a large number of SVs (>100 SVs per sample on average). Likewise, this number was high in uterine (UCEC) and stomach (STAD) tumors (107 and 125 per sample, respectively). Few SVs were found in copy number stable cancers such as thyroid carcinoma (THCA) (26) (∼8 per sample on average). Similarly, clear cell and chromophobe kidney carcinomas (KIRC and KICH) and low-grade glioma (LGG) tumors exhibited a small number of SVs (Fig. 2A).

We next investigated the relationship between somatic CNA and SV breakpoints. On average, 34% of CNA breakpoints (segments with absolute log₂ amplitude ≥0.3, SNP6 platform) detected across all 600 tumors were confirmed in the SV data. Overlaps were reduced in CNA stable tumors such as KICH, which are dominated by arm-level events and therefore may lack breakpoints detectable by WGS (Fig. 2B). Absolute numbers of CNA and SV breakpoints were strongly correlated across tumors (Pearson’s r = 0.81; Fig. 2C). However, a small number of samples had elevated CNAs not reflected in the SV data (Fig. 2C), again explainable by the predominance of arm-level changes in these tumors (Fig. S2). These results further support the validity of the SV data and show that both approaches are effective for quantifying overall chromosomal instability.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S2.

Copy number profile of arm-level dominant tumors. Screenshot from Integrative Genomics Viewer (IGV) showing the copy number changes in cases with mostly arm-level events. The number of CNV breakpoints in these cases is much higher than SV breakpoints detected using WGS. Chromosome-level or arm-level CN change is hard to detect using WGS. GBM, glioblastoma.

The overlap between CNA and SV data at the level of complete events (for example, both breakpoints for an amplification segment) was considerably lower (∼10%) compared with breakpoint-level overlaps. In many cases, this was explained in the SV data by a more complex structural basis than immediately evident from CNA profiles. Examples include apparent closely spaced large copy number deletion segments, revealed by the SV data to be due to a large (for example, whole-chromosome) deletion event superimposed on several small-amplitude–neutralizing tandem duplications (Fig. 2D). Similarly, apparent clustered copy number amplification segments could sometimes be explained by a larger duplication in combination with several small deletions (Fig. 2E). These results emphasize that care needs to be taken when drawing conclusions about structural chromosomal changes based on CNA data alone.

Of all CNA deletion segments (log₂ amplitude ≤−0.3) with a clear correspondence in the SV data, we found that 97% were categorized as deletions in the SV analysis (Fig. 2F). Likewise, 90% of amplified regions were due to tandem duplication, whereas only 5% each overlapped with deletions or inversions.

Global Impact of Somatic SVs on Tumor mRNA Levels.

For somatic SVs to be functional in cancer, they need to have an influence on transcription in some way, by altering mRNA structures or levels. We found that a large fraction (56%) of SV breakpoints were in genic regions, and that on average 112 genes per sample overlapped with at least one SV breakpoint, showing that SVs have considerable potential to alter tumor transcriptional output. A known mechanism for gene activation by SVs is the swapping of strong and weak promoters in the context of gene fusions. This is described for many individual genes, both with and without alteration of the protein-coding region (6, 27, 28), but it remains unclear to what extent such events have a measurable global influence on tumor expression profiles.

Based on the SV data, we identified local rearrangements predicted to result in a fusion gene involving a strong promoter juxtaposed to a normally weakly expressed gene. We thus considered only gene pairs where the 5′ partner was expressed at a considerably higher average level within a given cancer type (greater than fourfold difference; SI Methods). We found that these criteria were often predictive of strong induction of the 3′ fusion partner: In 62/313 cases (20%), we observed induction of the 3′ partner above threefold in the structurally affected tumor compared with other tumors of the same cancer type (Fig. 3A). This should be compared with only 8% and 3% for fusions involving incorrect gene orientation or a weakly expressed promoter, respectively (Fig. 3A; P = 1.1e-5 compared with incorrect orientation, P = 2.7e-7 compared with weak promoters, and P = 1.3e-6 compared with randomly picked cases; two-sided Kolmogorov–Smirnov test).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Impact of SVs on gene expression levels in tumors. (A) Gene fusion events leading to swapping of strong and weak promoters have a global influence on mRNA expression levels in tumors. Cumulative distribution function (CDF) plot showing log₂-transformed expression differences for genes predicted in the SV data to be activated due to fusion with a stronger promoter, comparing the SV-affected sample with the mean of other samples in the same cancer type. Control analyses involving fusions not predicted to cause activation of the 3′ partner (incorrect strand orientation or weakly expressed fusion partner) are shown by dashed lines, in addition to a randomized set (randomly picked sample within the relevant cancer type, median of 100 randomizations). Inversions and deletions were used for this analysis. (B) Global influence of promoter-proximal SVs on mRNA levels. CDF plot showing expression differences for genes harboring an upstream SV event (−10 to 0 kb), comparing the affected sample with the mean of other samples in the same cancer type. Genes with SVs in a region farther upstream (−110 to 100 kb) and randomized controls are shown. (C) Expression vs. copy number change for TERT in KICH. Samples with promoter-proximal SV breakpoints are marked in orange. P values were calculated using the Wilcoxon rank-sum test comparing the expression of the altered tumors with other samples. (D–F) Same as C for KIRC, READ, and SKCM, respectively. Samples with TERT promoter mutations are indicated in purple.

Among the 62 cases showing induced expression were several confirmatory observations (Fig. S3), including prominent activation of ERG in prostate cancer (29) as well as RET (30) and ALK (31) in THCA, by fusion with strong promoters from TMPRSS2, CCDC6, and STRN, respectively. Notable was strong (8.1-fold relative mean) mRNA induction of protein kinase C beta (PRKCB) through fusion with USP7, a highly expressed deubiquitinase, in colorectal carcinoma (COAD). PRKCB is a protein kinase linked to tumor growth and survival (32) that is druggable by the small-molecule inhibitor enzastaurin (33). Although fusions involving PRKCB were recently detected in fibrous histiocytoma (34), they have not been described in COAD. Additionally, among the induced transcripts were several long noncoding RNAs (lncRNAs), including the intergenic RP11-521D12.1, strongly (39.9-fold) induced through fusion with ASAP2 in BRCA. Although the functional relevance remains to be established, this establishes promoter swapping as a mechanism for activation of lncRNAs in cancer.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S3.

Promoter hijacking due to gene fusions resulting in mRNA induction. (A) Cases caused by deletions. The box plots show mRNA levels of fusion partners in relevant cancers. Red dots represent the affected samples. 3′ partners of fusions are shown in green. Only cases where the 3′ partner in the affected sample was expressed above threefold compared with other samples of the same cancer type are shown here. (B) Same as A but for inversions. Here, weakly expressed partners of the fusions are shown in green. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers.

A few recent studies show that structural rearrangements can contribute to transcriptional activation through shuffling of cis-regulatory elements near the TSS of genes, without altering mRNA structure (35, 36). This has been described in individual cases, but the impact on a more global scale remains unknown. To investigate this, we identified cases of genes being affected by upstream SVs (−10 to 0 kb from the TSS). We next considered the relative expression levels for these genes in affected samples compared with unaffected samples of the same cancer type while avoiding confounding effects from copy number changes and chromothripsis (SI Methods and Fig. S4). Cases with SVs in a region −110 to 100 kb upstream, or randomly picked cases, were considered as controls. We found that SV-affected genes/cases showed a weak trend toward overall increased mRNA levels (Fig. 3B; P = 0.069 and P = 0.0099, respectively, compared with the control sets). Although it remains to be determined to what extent the observed events are under positive selection, the results suggest upstream SVs may often contribute to transcriptional activation of genes in human tumors.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S4.

Relationship between chromotripsis and gene expression levels. The cumulative distribution function plot shows relative expression levels for genes/samples affected by chromthripsis. For a given gene, these relative levels were determined by comparing the expression level in the affected sample to the mean level of all samples with the same cancer type. Similarly, relative levels for genes/cases affected by overlapping simple SVs are shown. Randomly picked samples where it was affected by neither chromothripsis nor simple SVs were considered as controls. One hundred randomizations were performed, and the plot shows the mean. The plots show log₂-transformed relative expression levels.

Notably, among top cases showing increased expression in relation to upstream SVs (40 cases above fourfold) was telomere reverse transcriptase (TERT) in two different cancers, KICH and KIRC (Fig. S5). The data from KICH, where TERT upstream SVs were found in two samples with high TERT mRNA levels, confirm recently published results (7) (Fig. 3C). In KIRC, the sample with the highest TERT expression also harbored an upstream SV, whereas an additional affected sample lacked mRNA elevation in this cancer (Fig. 3D). Although not considered in the global screen due to simultaneous copy number gain, TERT upstream SVs were found also in one rectal adenocarcinoma (READ) and one melanoma (SKCM) sample (Fig. 3 E and F), the latter notably lacking TERT promoter mutations. These results further support a role for upstream SVs in activation of TERT in cancer.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S5.

Genes having an SV breakpoint within 10 kb upstream of the TSS and showing mRNA induction in the affected sample. mRNA level and copy number amplitude are shown for each sample in a given cancer type, with SV-affected samples indicated by orange dots. Only genes with at least a fourfold increase in the SV-affected sample are shown here.

Improved Gene Fusion Detection by Combined DNA/RNA Analysis.

Recent studies have relied on RNA-seq data to detect fusion genes in comprehensive tumor materials (37, 38). However, similar to SV analysis, RNA-seq–based fusion detection is prone to both false positives and negatives (39), for example due to read-mapping issues or misinterpreted germ-line events (40). Therefore, we explored whether fusion gene detection can be improved by combining RNA and DNA data.

We applied FusionCatcher (41) to identify fusion transcripts in 431/600 samples with available paired-end RNA-seq data. In total, 9,657 fusions (on average 22.4 per sample) were reported, of which 1,467 were in-frame (Dataset S2). In parallel, we annotated all genomic SV breakpoints and identified 9,854 predicted gene fusions, of which 740 were inversions, 1,240 were deletions, 2,511 were duplications, and 5,363 were interchromosomal translocations (Fig. 4A). Of these, 5,209 (52%) involved genes fused in the correct orientation, thus having the potential to produce a fusion transcript. We found that only 299 unique pairs of genes were detected in both DNA and RNA data, of which 97 were in-frame according to FusionCatcher (Fig. S6). In this high-confidence set with dual DNA/RNA support, there were 78 genes recurrently involved in fusions in more than one sample (Table S1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Combined DNA/RNA approach improves detection of gene fusions. (A) Venn diagram (not proportional) showing the number of fusions detected by analysis of RNA-seq data and in the genomic SV data, as well as overlaps between the two. (B) Enrichment of known cancer genes (from the Cancer Gene Census) among the top recurrently fused genes identified by analysis of RNA-seq data, genomic SVs, or a combined DNA/RNA approach. The graphs show the cumulative fraction of known cancer genes as a function of gene rank. Rank was determined for individual genes based on the frequency of recurrence, taking into account all fusions in which a given gene was predicted to participate. (C) A fusion between EIF2AK4 and BUB1B in one PRAD tumor, due to a deletion in chromosome 15. The deleted region is shown in blue, whereas nonaltered regions are shown in gray. mRNA levels [represented as reads per kilobase of transcript per million mapped reads (RPKM)] for both fusion partners across all PRADs are shown, with the affected sample indicated in red.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S6.

List of fusions detected by both SV-based and RNA-seq–based analyses. For each fusion, the underlying structural basis is shown. Bars on the right show the number of recurrences across all tumors, color-coded by cancer type similar to Fig. 2A. Genes listed in cancer gene census are shown in bold.

Next, we evaluated the two approaches, either alone or in combination, by investigating the enrichment of known cancer drivers in the Cancer Gene Census (CGC) database among the top recurrently fused genes (42). We found that neither method alone showed notable CGC enrichment, although the SV-based analysis performed slightly better (Fig. 4B). However, combining both approaches resulted in a notable CGC enrichment among top-scoring genes, with 14 of the 78 recurrently fused genes (18%) being listed in the CGC (random expectation ∼2.5%). This suggests that DNA and RNA approaches to fusion gene detection complement each other, such that a combination of both provides more relevant results.

In the set of high-confidence detections with dual DNA/RNA support (Fig. S6) were several known functional fusions, including TMPRSS2–ERG in 6/20 prostate adenocarcinoma tumors (PRAD) and RET–CCDC6 in THCA (2/34 tumors). Results also included EIF2AK4 fused to BUB1B, a mitotic checkpoint kinase, by a deletion in one PRAD tumor, with strong transcriptional induction of BUB1B in this sample (Fig. 4C).

Novel Mechanism for Nuclear Factor, Erythroid 2 Like 2 Activation in Cancer.

Notable among the high-confidence fusions was paired box 8 (PAX8) fusing with nuclear factor, erythroid 2 like 2 (NFE2L2) in one THCA tumor, revealed by the SV data to be due to tandem duplication on chromosome 2 (Fig. 5A). This event was associated with potent NFE2L2 transcriptional activation (Fig. 5B) as well as loss of the KEAP1 interaction domain in NFE2L2 (Fig. 5A). NFE2L2 (NRF2) is a transcription factor that coordinates the cellular response to oxidative stress, which is activated in some tumors by mutations that disrupt interaction with its inhibitor KEAP1 (43). However, NFE2L2 fusions involving loss of the KEAP1 interaction domain have not been described, and the result is thus suggestive of a novel mechanism for NFE2L2 activation in cancer.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Activation of the NFE2L2 pathway by a PAX8–NFE2L2 fusion gene in thyroid carcinoma. (A) PAX8 fused to the transcription factor NFE2L2 as a result of tandem duplication on chromosome 2. Nonaltered regions are shown in gray, and the tandem duplication is indicated in red. Major domains are indicated. The 5′ end of PAX8 replaces the inactivating KEAP1 interaction domain in NFE2L2. (B) Fusion with PAX8 is associated with strong transcriptional induction. The plot shows mRNA levels for both fusion partners in all thyroid carcinoma samples, with the affected tumor indicated in red. (C) Strong activation of NFE2L2 targets in the tumor expressing the PAX8–NFE2L2 fusion transcript. mRNA levels of NFE2L2 and four canonical target genes (GCLC, GCLM, TXNRD1, and NQO1) are shown across 572 samples with available exome data from TCGA. The PAX8–NFE2L2–expressing tumor, as well as one additional tumor identified as having a predicted NFE2L2-activating mutation in KEAP1, is indicated in green. P values were calculated using Student’s t test, comparing the expression of the two altered tumors with the remaining samples. (D) Transcriptome profiling of TRAMP-C1 murine prostate carcinoma cells transduced with a lentivirus expressing the PAX8–NFE2L2 fusion transcript or a control lentivirus. Volcano plot showing 534 genes significantly induced and 215 genes repressed at q < 0.01. (E) Enriched gene sets (P_nominal < 0.01) calculated by GSEA. The familywise error rate (FWER) calculated by GSEA is shown on the x axis, and the dashed line represents P_FWER = 0.05. WikiPathway (57) gene sets were used for this analysis. (F) Enrichment plot for the transcriptional activation by NRF2 gene set. (G) Altered expression of NFE2L2 target genes Gclc and Gclm was confirmed using RT-qPCR. Expression levels are shown relative to ubiquitin (mUb). P values were determined using Student’s t test. Bars indicate the mean, and error bars indicate SD.

To further evaluate the PAX8–NFE2L2 fusion, we investigated the expression of known target genes downstream of NFE2L2. We found that known NFE2L2 targets were strongly induced in the tumor that carried the PAX8–NFE2L2 fusion, supporting its function as a potent activator of the NFE2L2 transcriptional program (Fig. 5C). By considering an expanded set of THCA tumors with available exome sequencing data, we identified one tumor with a KEAP1 mutation, likewise associated with strong induction of NFE2L2 target mRNAs (Fig. 5C). This further supports that NFE2L2 activation may contribute to tumorigenesis in some thyroid carcinomas.

We next synthesized the PAX8–NFE2L2 fusion transcript and cloned it into a lentiviral vector (Fig. S7), and used this to overexpress it in the murine prostate adenocarcinoma cell line TRAMP-C1 (SI Methods). We used high-coverage RNA sequencing (>30 million reads per sample) to compare the expression profiles of cells expressing the PAX8–NFE2L2 fusion transcript with cells transduced with an empty vector. We observed a prominent transcriptional response dominated by inductions (534 induced and 215 repressed genes at q < 0.01; Fig. 5D).

• Download figure
• Open in new tab
• Download powerpoint

Fig. S7.

PAX8–NRF2 fusion gene constructs.

To identify affected pathways, we next performed gene set enrichment analysis using GSEA (44). Whereas no gene sets were significantly down-regulated, four sets were enriched among induced genes (adjusted P < 0.05; Fig. 5E). Notably, these included experimentally determined targets of NFE2L2 (“transcriptional activation by NRF2”; Fig. 5F), as well as “glutathione metabolism,” “oxidative stress,” and “arylhydrocarbon receptor” signaling, all of which are known to be regulated by NFE2L2 (43, 45). Induction of the canonical NFE2L2 targets Gclc and Gclm was confirmed by quantitative RT-PCR (2.4- and 3.1-fold, respectively; P < 0.022), further showing that the PAX8–NFE2L2 fusion protein is a capable activator of NFE2L2 target mRNAs (Fig. 5G).