"Computer Science, Department of"@en . "Medical Genetics, Department of"@en . "Medicine, Faculty of"@en . "Molecular Medicine and Therapeutics, Centre for"@en . "Pathology and Laboratory Medicine, Department of"@en . "Science, Faculty of"@en . "Non UBC"@en . "DSpace"@en . "Genome Biology. 2015 Apr 23;16(1):84"@en . "Child and Family Research Institute"@en . "Mathelier et al."@en . "Mathelier, Anthony"@en . "Lefebvre, Calvin"@en . "Zhang, Allen W."@en . "Arenillas, David J."@en . "Ding, Jiarui"@en . "Wasserman, Wyeth W."@en . "Shah, Sohrab P."@en . "2016-01-04T18:43:55Z"@* . "2015-04-23"@en . "Background:\r\n With the rapid increase of whole-genome sequencing of human cancers, an important opportunity to analyze and characterize somatic mutations lying within cis-regulatory regions has emerged. A focus on protein-coding regions to identify nonsense or missense mutations disruptive to protein structure and/or function has led to important insights; however, the impact on gene expression of mutations lying within cis-regulatory regions remains under-explored. We analyzed somatic mutations from 84 matched tumor-normal whole genomes from B-cell lymphomas with accompanying gene expression measurements to elucidate the extent to which these cancers are disrupted by cis-regulatory mutations.\r\n \r\n \r\n Results\r\n We characterize mutations overlapping a high quality set of well-annotated transcription factor binding sites (TFBSs), covering a similar portion of the genome as protein-coding exons. Our results indicate that cis-regulatory mutations overlapping predicted TFBSs are enriched in promoter regions of genes involved in apoptosis or growth/proliferation. By integrating gene expression data with mutation data, our computational approach culminates with identification of cis-regulatory mutations most likely to participate in dysregulation of the gene expression program. The impact can be measured along with protein-coding mutations to highlight key mutations disrupting gene expression and pathways in cancer.\r\n \r\n \r\n Conclusions\r\n Our study yields specific genes with disrupted expression triggered by genomic mutations in either the coding or the regulatory space. It implies that mutated regulatory components of the genome contribute substantially to cancer pathways. Our analyses demonstrate that identifying genomically altered cis-regulatory elements coupled with analysis of gene expression data will augment biological interpretation of mutational landscapes of cancers."@en . "https://circle.library.ubc.ca/rest/handle/2429/56170?expand=metadata"@en . "Mathelier et al. Genome Biology (2015) 16:84 DOI 10.1186/s13059-015-0648-7RESEARCH Open AccessCis-regulatory somatic mutations andgene-expression alteration in B-celllymphomasAnthony Mathelier1, Calvin Lefebvre2,3, Allen W Zhang1,3, David J Arenillas1, Jiarui Ding2,4,Wyeth WWasserman1* and Sohrab P Shah2,5*AbstractBackground: With the rapid increase of whole-genome sequencing of human cancers, an important opportunity toanalyze and characterize somatic mutations lying within cis-regulatory regions has emerged. A focus onprotein-coding regions to identify nonsense or missense mutations disruptive to protein structure and/or functionhas led to important insights; however, the impact on gene expression of mutations lying within cis-regulatoryregions remains under-explored. We analyzed somatic mutations from 84 matched tumor-normal whole genomesfrom B-cell lymphomas with accompanying gene expression measurements to elucidate the extent to which thesecancers are disrupted by cis-regulatory mutations.Results: We characterize mutations overlapping a high quality set of well-annotated transcription factor binding sites(TFBSs), covering a similar portion of the genome as protein-coding exons. Our results indicate that cis-regulatorymutations overlapping predicted TFBSs are enriched in promoter regions of genes involved in apoptosis orgrowth/proliferation. By integrating gene expression data with mutation data, our computational approachculminates with identification of cis-regulatory mutations most likely to participate in dysregulation of the geneexpression program. The impact can be measured along with protein-coding mutations to highlight key mutationsdisrupting gene expression and pathways in cancer.Conclusions: Our study yields specific genes with disrupted expression triggered by genomic mutations in either thecoding or the regulatory space. It implies that mutated regulatory components of the genome contribute substantiallyto cancer pathways. Our analyses demonstrate that identifying genomically altered cis-regulatory elements coupledwith analysis of gene expression data will augment biological interpretation of mutational landscapes of cancers.BackgroundTumor genome analyses for mutation and cancer genediscovery have focused primarily on the protein-codingexons, spanning approximately 2% of the genome, as theyare readily interpreted and easy to delineate. Large-scaleconsortia such as The Cancer Genome Atlas have com-pleted interrogation of the protein-coding genome andrevealed the mutation prevalence of previously known*Correspondence: wyeth@cmmt.ubc.ca; sshah@bccrc.ca;1Department of Medical Genetics, Centre for Molecular Medicine andTherapeutics, Child and Family Research Institute, University of BritishColumbia, 950 West 28th Avenue, V5Z 4H4, Vancouver, BC, Canada2Department of Molecular Oncology, British Columbia Cancer Agency, V5Z1L3 Vancouver, BC, CanadaFull list of author information is available at the end of the articlecancer genes across the major tumor types, in additionto discovery of previously unknown biological processesdisrupted by somatic mutations [1]. However, synthesis ofthe vast analyses of The Cancer Genome Atlas projectshas revealed a discovery gap in the search for new cancergenes [2]. We assert this gap can be partially filled throughanalysis of the non-coding genome. In germline geneticdisease studies, evidences for the impact of variationsin the non-coding space of the human genome, includ-ing in cis-regulatory loci, on human phenotypes haveaccumulated over recent decades [3]. Gene-expressionregulation occurs through multiple layers, one of themmediated by DNA-binding transcription factors (TFs).Disruption of sequence-specific TF binding sites (TFBSs)\u00C2\u00A9 2015 Mathelier et al. This is an Open Access article distributed under the terms of the Creative Commons Attribution License(http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium,provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.Mathelier et al. Genome Biology (2015) 16:84 Page 2 of 17has been linked to numerous genetic disorders. For exam-ple, mutation within a HNF4A binding site upstreamof the Factor IX gene is associated with hemophilia BLeyden [4], and alteration of a GATA binding site ina regulatory region upstream of the platelet glycopro-tein gene causes Gilbert\u00E2\u0080\u0099s syndrome [5]. More recently,human melanoma studies have revealed highly recurrentmutations in the TERT promoter, potentially impactingregulatory elements [6-8].The emergence of whole-genome sequencing studies incancer has highlighted the importance of analyzing muta-tions lying within cis-regulatory elements [6,7,9-13]. How-ever, the global relationship between somatic nucleotidevariations and the creation or disruption of TFBSs impact-ing gene expression in cancer is largely unknown [14].Several attempts have been made to predict the degreeto which mutations disrupt TFBSs with TF binding pro-files (that is, position weightmatrices, PWMs), populationgenetics, phylogenetic footprinting and experimental data(DNase-seq, epigenetic, etc.) [12,13,15-19]. Mutations atcritical positions of TF binding profiles (correspondingto high information content) are the most deleterious forTF\u00E2\u0080\u0093DNA binding [17], thus modelling impact by PWMsis an effective strategy for predicting the impact of a muta-tion [20]. However, mutations at the more variable, lowinformation content positions of TFBSs can also be func-tionally constrained [17,21]. Furthermore, relating a TFBSto the gene(s) it regulates presents additional challengesin predicting cis-regulatory mutations impacting geneexpression. A common simplifying assumption is that aTFBS regulates its closest gene. This first approxima-tion does not consider distal regulation; however, recentanalyses of chromatin immunoprecipitation coupled tohigh-throughput sequencing data sets [22] (the so-calledchromatin immunoprecipitation sequencing or ChIP-seqprocedure) showed there was accurate prediction of TFgene targets using this approach [23]. In this study, wepropose that direct measurements of cis-regulatory muta-tions and gene expression in the same tumor samples willoptimally identify mutations in TFBSs impacting the geneexpression program in cancer cells.Ultimately, interpretation of mutations in cis-regulatoryregions of the genome requires accurate annotation ofTFBSs. We have taken the approach of coupling experi-mental data to targeted computational analysis with TFbinding profiles. The ENCODE project [24] and otherindependent analyses provide a rich resource for locat-ing the key regulatory positions by providing genomicregions bound by TFs derived from ChIP-seq data sets.This provides unprecedented means by which to investi-gate altered TFBSs and gene regulation in cancer samples[12]. Together with matched gene expression profiles,analysis of mutations in well-annotated TFBS lying inChIP-seq regions provides a robust set of complementarymeasurements to study the characteristics of dysregula-tion through mutation of cis-regulatory elements.We set out to characterize the impact of cis-regulatorysomatic mutations on gene expression. We focused ontwo cohorts of patients with B-cell lymphomas (BCLs)[25,26], for which 84 trios (the cancer genomes, matchedpatient normal genomes and RNA expression from RNA-sequencing (RNA-seq)) were analyzed. We identifiedcancer-specific somatic mutations across the genome,considering single nucleotide variants (SNVs) and smallinsertions and deletions (indels) and centered our anal-ysis on cis-regulatory elements corresponding to TFBSspredicted within TF-bound regions delineated as ChIP-seq peaks. The regulatory space defined in this analysis bypredicted TFBSs within ChIP\u00E2\u0080\u0099ed regions covered approx-imately 2% of the human genome. We analyzed thelocation of mutations overlapping TFBSs and revealedthat they frequently target promoter regions of apop-totic genes. Integrative analysis of the mutations andgene expression data fromRNA-seq highlighted candidateregulatory-disrupting variations as potentially alteringexpression of genes involved in cancer development.Mutations in cis-regulatory elements were frequent,and high-quality candidates in the regulatory set wereobserved to target genes mutated in the coding spacein other samples. We conclude that analysis and inter-pretation of the cis-regulatory genome of cancers willmeaningfully augment biological discovery in future stud-ies, resulting in novel mechanistic insight into the genesismalignant phenotypes.ResultsWe analyzed somatic mutations extracted from whole-genome sequencing of 84 BCL samples along with match-ing normal samples from the same individuals. The fullset of samples is composed of 40 diffuse large B-cell lym-phomas (DLBCLs) (cohort 1) and 44 patients of mixedhistology (cohort 2: 14 Burkitt lymphomas, 15 DLBCLs,1 primary mediastinal large B-cell lymphoma (PMBCL)and 14 follicular lymphomas). RNA expression profilingdata (from RNA-seq) were also available for each of the84 cancer samples plus 62 additional lymphoma samples(52 associated with cohort 1 and 10 with cohort 2). SNVand indel analyses of the data sets from the two cohortswere performed independently as data were derivedfrom different sequencing methods. Somatic mutations incohort 1 were identified with MutationSeq [27], whereasmutations from cohort 2 were retrieved from the origi-nal publication [26]. In aggregate, we observed 406,611SNVs (from 146 to 31,874 per sample; mean = 10, 165,median = 7, 821 and standard deviation (sd) = 6, 995)and 15,739 indels (from 65 to 4,810 per sample; mean =393, median = 222 and sd = 735) in samples fromMathelier et al. Genome Biology (2015) 16:84 Page 3 of 17cohort 1 and 282,636 SNVs (from 1,242 to 37,987 persample; mean = 6, 424, median = 3, 577 and sd =7, 165), and 8,080 indels (from 67 to 871 per sample;mean = 184, median = 136 and sd = 142) in samplesfrom cohort 2 (Figure 1). The distribution of mutationsand mutation types over the samples followed a similarpattern in the data sets from the two cohorts includingthe maximum number of mutations, >30,000 (Figure 1and Additional file 1: Figure S1). Histological types fromcohort 2 clustered by the number of mutations. Namely,Burkitt lymphomas harbored fewer mutations than fol-licular lymphomas, while DLBCLs harbored the highestnumber of mutations, consistent with the number ofmutations observed within cohort 1 (Figure 1).Defined cis-regulatory elements showed a higher mutationrate than protein-coding exons but were less mutated thantheir flanking regionsWe began by first identifying mutations lying withincis-regulatory elements. We considered TFBSs to be cis-regulatory elements and mutations overlapping TFBSswere assumed to be cis-regulatory mutations. TFBSs werepredicted within ChIP-seq peak regions, collected frommultiple cell types and tissues, at whole-genome scaleusing TFBS profiles from the JASPAR database [28] (see\u00E2\u0080\u0098Materials and methods\u00E2\u0080\u0099). We used 477 ChIP-seq datasets (collected for the last update of the JASPAR database[28]) to predict TFBSs associated with 103 TFs (107 JAS-PAR profiles). Predicted TFBSs covered 76,160,599 bpFigure 1 Distribution of the number of mutations per sample in cohorts 1 and 2. (A) Number of SNVs (blue) and indels (red) on the y-axis are givenfor all the samples in cohort 1 on the x-axis. The samples are ordered from the least number of mutations (left) to the most (right). (B) The same typeof distribution for the samples in cohort 2. Sample names on the x-axis are color-coded by tumor subtype: Burkitt lymphomas (green), diffuse largeB-cell lymphomas (DLBCLs, black), primary mediastinal large B-cell lymphomas (PMBCLs, gray) and follicular lymphomas (FLs, red). The same y-axisscale has been used for (A) and (B) for comparison. DLBCL, diffuse large B-cell lymphoma; FL, follicular lymphoma; PMBCL, primary mediastinal largeB-cell lymphoma; SNV, single nucleotide variant.Mathelier et al. Genome Biology (2015) 16:84 Page 4 of 17of the human genome in our analysis. As expected, weobserved a very strong enrichment for predicted TFBSs inpromoter regions of protein-coding genes with >5 timesmore nucleotides covered by TFBSs than expected bychance (11,073,418 bp from TFBSs overlapping promot-ers covering 85,296,239 bp). The portion of the genomecovered by predicted TFBSs represents approximately 2%of the chromosomes. We noted the cis-regulatory space,which only overlapped the protein-coding space by 6%,covered a proportion of the human genome similar toprotein-coding exons.From the 422,350 mutations predicted from cohort 1,8,184 (approximately 2%) overlapped TFBSs. Likewise,6,608 of 290,716 mutations (approximately 2%) over-lapped TFBSs in cohort 2. By comparison, 4,990mutations (approximately 1%) and 5,098 mutations(approximately 2%) overlapped protein-coding exons incohorts 1 and 2, respectively. SNV mutation rates werehigher in TFBSs than protein-coding exons for 38 (95%)cohort 1 and 25 (57%) cohort 2 samples (Figure 2). Anal-ysis of 1,000 simulated genomes with a random mutationdistribution shows TFBSs with a higher mutation ratethan exons is expected by chance for 540 (respectively,616) samples in cohort 1 (respectively, cohort 2). Themajority of DLBCL (10 of 15) and follicular lymphoma (9of 14) cohort 2 samples showed a higher mutation rate inTFBSs than in exons; however, the reverse was observedfor Burkitt lymphomas (5 of 14) (Figure 2B). Indel muta-tion rates were similar in TFBSs and exons (Additionalfile 1: Figure S2).TFBSs were less mutated than their flanking regions inboth cohorts (Figure 3). Namely, 38/40 cohort 1 and 37/44cohort 2 samples exhibited lower SNV mutation rates inTFBSs compared to flanking regions. In contrast, only30/40 cohort 1 and 3/44 cohort 2 samples showed lowerSNV mutation rates in protein-coding exons compared toflanking regions (Figure 3). The difference observed forcohort 2 is consistent with the above stated comparison ofTFBS and exon SNV mutation rates. Local mutation ratesin TFBSs and exons were found to be similar to their flank-ing regions for indels (Additional file 1: Figure S3). Takentogether, these results indicate that predicted TFBSs havelower SNV mutation rates than their flanking regions inboth cohorts, while indels are more randomly distributed.Promoters of apoptotic genes are frequently targetedregions for cis-regulatory mutationsWe further explored the impact of cis-regulatory muta-tions by investigating their distribution along the humangenome. We sought to characterize the accumulationof mutations in TFBSs lying within the promoters ofgenes implicated in pathways known to be disrupted incancer development. We quantified mutation rates in 1-kb-long sliding windows across the genome, identifyingwindows where at least three mutations were found (thetwo cohorts were analyzed independently and combined).Frequently mutated regions are significantly enrichedfor promoters of protein-coding genes (Figure 4). Namely,135 mutations in frequently targeted regions were \u00E2\u0089\u00A42 kbaway from a protein-coding gene\u00E2\u0080\u0099s transcription startsite (TSS) using the samples from cohort 1 (represent-ing approximately 49% of all 273 mutations found infrequently targeted regions, P = 1.16 \u00C3\u0097 10\u00E2\u0088\u009275, hyper-geometric test with 680 mutations overlapping TFBSs inpromoters out of 8,185 in TFBSs). In cohort 2 samples, 348mutations in frequently targeted regions were within pro-moters (approximately 65% of the 534 found in frequentlytargeted regions, P = 3.28 \u00C3\u0097 10\u00E2\u0088\u0092156, hypergeometric testwith 1,102 mutations overlapping TFBSs in promoters outof 6,608 in TFBSs). We compiled the set of mutationsfound in the frequently targeted regions within promot-ers and extracted the closest protein-coding gene to eachmutation. Twelve genes were frequently targeted in bothcohorts independently (Figure 4A,B), including BCL2,BCL6, BCL7A, CD74 and CIITA, all listed as oncogenesin the Cancer Gene Census [29] and known to be involvedin lymphomagenesis. An additional 13 genes (ARID2,BCL2L11, BZRAP1, EPS15, HIST1H2BG, ID3, IGLL5,IL2R1, IRF1, KIAA0226L, NEDD9, RARS and ZNF860)from combined cohort analysis (Figure 4C) had not beenpreviously described as aberrant somatic hypermutatedregions [30]. Six of these genes (ARID2, BCL2L11, EPS15,IL2R1, NEDD9 and ZNF860) were exclusively mutatedin their promoters (that is, no mutations in exons wereobserved). Our data indicated for the first time that ID3(recurrently mutated in Burkitt lymphomas [26]) can betargeted through TFBS mutations in its promoter region.Thus, both exonic and promoter portions of the gene arerecurrently mutated, suggesting complementary geneticmechanisms for gene disruption.Mutations in frequently targeted regions overlapped fiveenhancers in cohort 1 and five enhancers in cohort 2(there were two enhancers in common: intronic enhancersof BIRC3 and ST6GAL1). Four of the enhancers targetedin cohort 1 are intronic enhancers for the genes BCL2,BCL7A, BIRC3 and ST6GAL1 while the fifth is locatedin the intergenic region between BCL6 and LPP. All fiveenhancers found in cohort 2 are intronic enhancers ingenes BCL2, BIRC3, CIITA, IGLL5 and ST6GAL1. Allthese genes have already been associated with hyper-mutated regions in BCLs (BCL2, BCL6, BCL7A, BIRC3,CIITA and ST6GAL1) [30] or listed in the Cancer GeneCensus (BCL2, BCL6, BCL7A, BIRC3, CIITA, IGLL5 \u00E2\u0080\u0093 inthe IGL@ locus \u00E2\u0080\u0093 and LPP).To synthesize our observations from the gene level,we next analyzed genes with frequently targeted pro-moters through pathway enrichment analysis. All thegenes highlighted in Figure 4 were submitted to EnrichrMathelier et al. Genome Biology (2015) 16:84 Page 5 of 17Figure 2 Comparison of the mutation rates in the cis-regulatory and protein-coding spaces. Only SNVs from cohort 1 (A) and cohort 2 (B) havebeen considered (see Additional file 1: Figure S2 for indels). TFBS mutation rates (y-axis) and protein-coding mutation rates (x-axis) are plotted for allthe samples in cohort 1 (A) and cohort 2 (B). Each triangle represents a sample and is color-coded depending on the tumor subtype as in Figure 1.Dashed gray lines represent the identity function (x = y). Blue lines represent the linear regressions computed from the samples in the two datasets. The equations corresponding to the linear regressions (y \u00E2\u0088\u00BC x) are written on top of the plots along with the computed r2 statistical measures.Dark gray areas surrounding the blue lines provide the 95% confidence region. The same x- and y-axis scales have been used for both cohort 1 (A)and cohort 2 (B). DLBCL, diffuse large B-cell lymphoma; FL, follicular lymphoma; PMBCL, primary mediastinal large B-cell lymphoma; SNV, singlenucleotide variant; TFBS, transcription factor binding site.[31]. Both cohorts were analyzed separately (genes fromFigure 4A,B) and combined (genes from Figure 4C). Weidentified enrichment (adjusted P < 0.05) for apoptoticprocesses (Figure 5 and Additional file 2) includingapoptosis, regulation of the B-cell apoptotic process andcell-type-specific apoptotic processes. The genes asso-ciated with the apoptotic terms are BCL2, BCL2L11,BIRC3, BTG1, CD74, IRF1, IRF4 and MYC. Moreover,Mathelier et al. Genome Biology (2015) 16:84 Page 6 of 17Figure 3 Local SNV mutation rates for TFBSs and protein-coding exons. Mutation rates are plotted for TFBSs and exons (x-axis) versus their 1-kbflanking regions on both sides (y-axis). Each sample from cohort 1 (A) and cohort 2 (B) is represented by a square for mutation rates in TFBSs (andtheir flanking regions) and a triangle for mutation rates in exons (and their flanking regions). Tumor subtypes are color-coded (see legend) like inFigure 1. Results are for SNVs. Figures corresponding to local indel mutation rates are provided in Additional file 1: Figure S3. DLBCL, diffuse largeB-cell lymphoma; FL, follicular lymphoma; PMBCL, primary mediastinal large B-cell lymphoma; SNV, single nucleotide variant; TFBS, transcriptionfactor binding site.we observed enrichment for B-cell and oncogenic relatedpathways (for example, the B-cell receptor signaling path-way, small cell lung cancer, regulation of B-cell prolifera-tion, lymphoma and leukemia) as shown in Figure 5 andAdditional file 2. Taken together, these results highlightthat apoptotic genes, and oncogenic processes in general,are frequently targeted by mutations within TFBSs foundat their promoter regions.Landscape of cis-regulatory mutations impacting geneexpression in B-cell lymphomasWe next assessed the impact of cis-regulatory mutationson gene expression. We used a novel probabilistic model,called xseq [32], to relate specific mutations to expres-sion disruption in pathways (see \u00E2\u0080\u0098Materials and methods\u00E2\u0080\u0099).The approach assesses the likely association of the pres-ence of mutations with observed deviations from neutralexpression measurements taken from the same tumor.The method takes as input a patient-gene expressionmatrix and a binary patient-gene mutation matrix andoutputs the probabilities that: (a) a mutated gene (overthe whole patient population) impacts gene expressionand (b) a patient-specific mutated gene impacts expres-sion in the patient. xseq was originally developed forgenes harboring mutations in their protein-coding exonsonly. Here, we extended xseq to highlight cis-regulatorymutations potentially deregulating transcriptional regula-tion. We encoded a gene as mutated in the patient-genemutation matrix when it was the closest gene to a cis-regulatory mutation and it was up- or down-regulated inthe mutated sample compared to other samples. With theapplied criteria, a TFBS was associated with a single genebut a gene might be associated with several TFBSs. Toprovide xseq with a complete view of mutated genes, weincorporated both genes mutated in their protein-codingregions and genes showing altered expression associatedwith mutations in their regulatory regions (Additionalfile 1: Figure S4 and \u00E2\u0080\u0098Materials and methods\u00E2\u0080\u0099). By combin-ing these tagged genes with expression data from RNA-seq, we used xseq to predict candidate mutated genesassociated with altered expression and linked to genes inbiological networks harboring altered expression in thesame samples.A total of 42 genes were predicted in cohort 1 sam-ples along with 5,412 biological network neighbors withaltered expression (Figure 6A). The same analysis appliedto cohort 2 samples led to 1,533 deregulated biologi-cal network genes connected to the 52 xseq-predictedgenes with altered expression associated with muta-tions (Figure 6B). The sets of genes captured by xseqalong with their deregulated neighbors were enrichedfor pathways related to cancer and cancer development(Figure 7A,B,C,D and Additional file 2). The sets of 5,554genes from cohort 1 and 1,585 genes from cohort 2had an intersection of 829 genes (Additional file 2).Functional enrichment analysis highlighted strong over-representation of cancer-related genes (Figure 7E,F andAdditional file 2), reinforcing the predictions from xseq asbeing involved in cancer development a posteriori. Notethat four genes in cohort 1 (HIST1H1B, RHOH, SGK1and ZFP36L1) and seven in cohort 2 (BCL6, DUSP2, ID3,FOXO1, MYC, PIM1 and SGK1) were associated with fre-quently targeted promoters (Figure 4) and predicted byxseq (Figure 6).Ranking the predicted genes by the number of sam-ples in which they were dysregulated highlighted knownMathelier et al. Genome Biology (2015) 16:84 Page 7 of 17Figure 4 Regions frequently targeted by somatic mutations overlapping cis-regulatory elements are enriched in promoters. A 1-kb-long windowwas slid with a 500-bp step along the human chromosomes and we recorded the number of overlapping mutations at each position. Thecorresponding histogram is given in the inner gray circles of the Circos plots where only positions containing at least three mutations have beenretained. The y-axis range for the histograms is [0, 40] for (A) and (B) and [0, 80] for (C). The outer circles contain the names of the genes closest tothe mutations in the considered windows if the mutation is at most 2 kb away from the TSS of the gene. Names highlighted in red correspond togenes shared between the analyses for the cohort 1 (A) and the cohort 2 (B) data sets. Genes highlighted in blue are specific to the analysis of themutations when combining somatic mutations from the two cohorts. Analyses have been applied to the set of mutations from cohort 1 (A),cohort 2 (B) and both cohorts combined (C). TSS, transcription start site.cancer driver genes such as MYC, TP53, ID3 and BCL6(Figure 6). Burkitt lymphomas tended to be segregatedfrom the other types of BCLs where MYC was predictedas a mutated gene with altered expression. MYC waspredicted by xseq for 11 samples, 10 of which areBurkitt lymphomas. In all of the 11 samples, MYCMathelier et al. Genome Biology (2015) 16:84 Page 8 of 17Figure 5 Functional enrichment analyses of genes associated with frequently mutated regions. Enrichr [31] functional enrichment analyses wererealized for the sets of genes listed in Figure 4 for cohort 1 (A), cohort 2 (B) and the two data sets combined (C) (see Additional file 2). The enrichedpathways with the 20 lowest Bonferroni corrected P values are shown for each category only where Bonferroni corrected P < 0.05 (see Additionalfile 2 for the complete Enrichr results). Each node of the graphs represents an enriched pathway where the color of a node represents its Bonferronicorrected P value. An edge between two nodes indicates that the pathways share genes. The larger the width of the edge, the larger the number ofshared genes. FDR, false discovery rate; GO, gene ontology; OMIM, online mendelian inheritance in man.up-regulation was observed (Additional file 1: Figure S5),in agreement with the oncogene function of MYC incancers [33].We next characterized the distribution of mutationsin genes impacting gene expression in protein-codingand cis-regulatory regions. We categorized each gene asassociated with: (1) a protein-coding mutation, (2) a cis-regulatory mutation or (3) both (Figure 6). Some geneswere predicted with altered expression and associatedwith mutations in their exons only (for example, TP53,Mathelier et al. Genome Biology (2015) 16:84 Page 9 of 17Figure 6 xseq results. Cancer genes predicted by the xseq tool from the cohort 1 (A) and cohort 2 (B) data sets. Each row corresponds to apredicted gene and each column to a cancer sample. When a gene is predicted in a specific sample, a colored box is drawn. Box colors indicate thetype of mutation associated with the gene (in protein-coding exons only, brown; in TFBS only, green; in protein-coding exon and TFBS, orange; in aTFBS only and predicted to disrupt the TFBS, pink; and in protein-coding exon and TFBS and predicted to disrupt the TFBS, purple). The histogramson the right sum the number of samples where the gene is predicted (using the same box colors). The histograms at the top sum the number ofgenes predicted by xseq in samples (using the same box colors). Cohort 2 sample names (B) are color-coded as defined in Figure 1. PC, proteincoding; TFBS, transcription factor binding site.RYR2 and SIN3A in cohort 2 and COL3A1, IRF8 andNRIP1 in cohort 1). We also observed multiple genes pre-dicted inmultiple samples consistent with alternatemech-anisms of alteration. For instance, HAS2 and ZFP36L1 inthe cohort 1 data set and MYC and BCL6 in the cohort 2data set were associated with mutations either in the cis-regulatory or the protein-coding spaces. For ID3, geneexpression alteration was associated with a mutation ina TFBS in the SA320932 sample whereas it was associ-ated with mutations in the exons in the two other samples(SA321012 and SA320818) (Figure 6B).Examples of genes with cis-regulatory mutationsassociated with expression dysregulationxseq analyses highlighted the specific mutations associ-ated with gene-expression dysregulation along with cas-cading effects on interactors through functional proteinassociation networks. For instance, specific SNVs werepredicted as deleterious for TFBSs and associated withexpression dysregulation of the genes HAS2 and GNA13(Figure 8 and Additional file 3). Recurrent dysregula-tion associated with cis-regulatory mutations was alsoobserved as exemplified in the promoter of BCL6 alongwith a potential cascading effect on interactors knownto be involved in cancer development (Additional file 1:Figures S7 and S8 and Additional file 3). As a last exam-ple, our approach highlighted SNVs in TFBSs associatedwith the promoters of ROBO1 for five DLBCL samples(Figure 6 and Additional file 1: Figure S9). ROBO1 wasdown-regulated in all of these five samples (Additionalfile 1: Figure S10). We hypothesize that ROBO1 is a tumorsuppressor (as suggested in [34-36]), whose dysregula-tion shows recurrent altered expression of its interactorsSOS1, SOS2 and RAC1, which are associated with car-cinogenesis [36,37] (Additional file 1: Figure S10 andAdditional file 3). These observations shed light on thesupposed tumor suppressor role of the ROBO1 gene. Wehighlight that ROBO1 might be down-regulated in someDLBCLs at the transcriptional level by cis-regulatorymutations since no mutations were found in the protein-coding space in these samples.DiscussionOur results reveal the importance of fully characterizingsomatic mutations in cis-regulatory regions of can-cer genomes. Whole-genome sequencing data fromMathelier et al. Genome Biology (2015) 16:84 Page 10 of 17Figure 7 Functional enrichment analyses of disrupted pathways. xseq-predicted genes along with their neighbors in biological pathways showingaltered expression were derived from the xseq analyses (see \u00E2\u0080\u0098Materials and methods\u00E2\u0080\u0099 and Additional file 2). Functional enrichment was performedwith Enrichr [31] using the genes obtained from the cohort 1 (A,B) and cohort 2 (C,D) data sets and their intersection (E,F) (see Additional file 2).The enriched terms from KEGG (A,C,E) and WikiPathways (B,D,F) with the 20 lowest Bonferroni adjusted P values are shown. Only terms with aBonferroni corrected P < 0.05 are conserved (see Additional file 2 for the complete Enrichr results). Each node of the graphs represents an enrichedpathway where the color of a node represents its Bonferroni corrected P value. An edge between two nodes indicates that the pathways sharegenes. The larger the width of the edge, the larger the number of shared genes. FDR, false discovery rate.Mathelier et al. Genome Biology (2015) 16:84 Page 11 of 17Figure 8 Examples of predicted cis-regulatory mutations potentially impacting gene expression. HAS2 (A,B,C,D,E) and GNA13 (F,G,H,I,J) havebeen predicted by xseq for samples RG116 and SA320848, respectively. In RG116, a CEBPA TFBS (TF binding profile in (A)) is predicted to bedisrupted (see reference and alternative sequences in (B) where the SNV is highlighted with the reference nucleotide in green and the alternative inpurple). Score differences between the reference TFBS and all possible alternative TFBSs are plotted in (C). The distribution of HAS2 expression fromRNA-seq data is plotted in (D) with an arrow pointing to the expression value in sample RG116. (E) represents the network of genes associated withHAS2, which are predicted to be either down- (blue) or up-regulated (red) in RG116. The higher the opacity, the stronger the down- orup-regulation. Similar plots are given in (F,G,H,I,J) for GNA13 in SA320848 with a potentially disrupted GATA3 TFBS. alt., alternative; ref., reference;SNV, single nucleotide variant.Mathelier et al. Genome Biology (2015) 16:84 Page 12 of 17lymphoma samples indicated somatic mutations impact-ing TFBSs and events associated with alteration in tran-scription. These results demonstrated that interpretationof mutations in cancer genomes will be substan-tially enhanced by consideration of mutations impact-ing cis-regulatory regions and with joint analysis ofgene expression data acquired from the same tumortissue.We expect our results to be an underestimate ofthe functional non-coding mutational landscape. Ourapproach relied on high-quality annotations of the cis-regulatory space, capitalizing on the availability of a largevolume of experimentally derived TF-DNA interactionsfrom ChIP-sequencing. We used manually curated TFbinding profiles from the JASPAR database to predictTFBSs within regions bound by ChIP\u00E2\u0080\u0099ed TFs. Although wecombine ChIP-seq experiments from multiple cell typesand conditions, the experimental ChIP-seq informationprovides the best current opportunity to focus on the non-coding space. The set of predicted TFBSs covered approx-imately 2% of the human genome, a similar proportionto the coding regions, and harbored lower SNV muta-tion rates than their surrounding regions. However, wecan expect that the robustly annotated regulatory spaceof the human genome will grow over the coming yearswith the availability of more antibodies, decreasing costsand broader coverage of cell types. As such, it is likelythat additional cis-regulatory regions will be found aber-rant in tumor genomes, allowing for more comprehensiveinterpretation of genome-wide somatic mutations drivingmalignant phenotypes.Detection of genes with altered expression due to thedisruption of regulatory TFBSs is the subject of ongo-ing research. There is a need for better prediction ofthe impact of mutations on TF\u00E2\u0080\u0093DNA binding affinity.Multiple approaches have been explored over the yearsto tackle this problem by considering the score dif-ference between reference and alternative sites [15] orthe decrease of the reference binding compared to thealternative binding score [12]. Here, we considered allmutations lying within TFBSs of potential interest, high-lighting the ones that are the most likely to disrupt TFBSswhere the alternative score was below a defined thresh-old. We suggest this approach is simple and conservative.Therefore, future improvements are likely to increase sen-sitivity when coupling ChIP-seq data to TFBS variantprediction.The prediction of TFBSs within ChIP-seq peaks is per-formed without considering the competing environmentbetween TFs with different specificities. For instance, thedown-regulation of BCL6 in SA320962 is associated witha mutation not predicted to disrupt the STAT3 TFBSwhile STAT3 is known as an activator. A hypothesis isthat competition between STAT3 and STAT5 occurs atthis binding site, since they recognize similar motifs anda previous study highlighted that STAT5 outcompetesSTAT3 for repressing the expression of BCL6 [38]. Themutation might then provide an advantage for STAT5binding at this location. Functional studies based on ouranalyses would be required to decipher the mechanismsunderlying the competition between TFs at TFBS locito understand further the impact of mutations on generegulation.We aimed to predict the mutations most likely to havean impact on gene expression through the disruption ofregulatory elements (TFBSs). The approach focused ongene expression alteration was built on top of the protein-coding changes to provide new insights into gene expres-sion dysregulation of cancer driver genes. We showedthat some genes have been predicted in multiple sam-ples using different mechanisms of dysregulation througheither protein-coding or TFBS alterations. We took anaive approach in this study to associate mutations/cis-regulatory elements with a gene. This approach is relevantwhen looking at promoter regions but we will ultimatelyrequire more information about the association of dis-tal regulatory elements with promoters. Cell-specific (orcell-type-specific) experimental profile comparisons andexpansion of chromatin conformation capture data setswill empower analysis linking distal regulatory elementsto their targets.With the forthcoming availability of cancer whole-genome sequence data coupled with gene expression dataat large scale, the analysis of non-coding cis-regulatoryelements will be critical for understanding cancer. Ourresults indicated this will be fruitful, yielding addi-tional cancer biology to aid in closing the discoverygap in large-scale studies that have focused exclusivelyon the protein-coding component of the genome. Wesuggest that combining the impact of mutations ontranscriptional regulation, protein products and post-transcriptional regulation at genome scales will empowercomprehensive biological interpretation of humanmalignancy.ConclusionsIn this report, we analyzed a set of approximately 700,000somatic SNVs and indels in 84 BCL samples to pro-vide an initial genome-scale foray into the analysis ofcis-regulatory mutations impacting gene expression incancer. By overlapping the somatic mutations with pre-dicted TFBSs within ChIP-seq regions, we looked atthe distribution of mutations overlapping TFBSs in thehuman genome. We highlighted that cis-regulatory muta-tions are frequently situated in promoters. The set ofgenes with promoters targeted bymutations within TFBSsare enriched for apoptosis-related and carcinogenesispathways. Finally, by combining mutation informationMathelier et al. Genome Biology (2015) 16:84 Page 13 of 17with gene expression from RNA-seq data, we predictedcancer genes with altered expression associated withmutations found either in exons or in TFBSs associatedwith the genes. The approach revealed samples wheregenes were potentially dysregulated through the disrup-tion of cis-regulatory elements and highlighted the impor-tance of interrogating the cis-regulatory genomic space forsomatic mutations in cancer.Materials andmethodsTranscription start site, exon and enhancer coordinatesTSS and exonic positions of protein-coding genes (tran-script accessors starting with NM_) have been retrievedfrom the UCSC hg19 Table Browser [39] by selecting theknownGene table from the RefSeq gene track. Enhancercoordinates were retrieved from [40].Cancer genome dataThe raw sequencing data for the DLBCL data set ofcohort 1 were retrieved from [25]. We obtained the cor-responding RNA-seq data from the same publication. Thealready processed sets of mutations and RNA-seq expres-sion level data for cohort 2 [26] were retrieved from theICGC data portal [41].Expression data computationExpression data for samples in cohort 1 were processedusing the Rsamtools andGenomicFeature [42] Bioconduc-tor [43] packages to generate gene expression levels fromthe RNA-seq raw data. We only considered genes withan official HGNC symbol [44]. Finally, genes with nullexpression over all the samples were filtered out. The finalset of official HGNC symbols for the considered genes canbe found in Additional file 4. We did not consider copynumber alteration information for both cohorts since thedata were not available for cohort 2.Single nucleotide variant predictionsSNVs were identified for cohort 1 samples using a modi-fied version of MutationSeq [27,45]. We filtered out SNVswith probability <0.9.Indel predictionsWe used Dindel [46] to call indels in the samples fromcohort 1. Dindel identified indels of length 1 to 50 bp. Fol-lowing Dindel\u00E2\u0080\u0099s manual recommendations, we providedDindel with the BAM file and the set of candidate indelsobtained from Pindel [47] for each sample. Default param-eters were used for Dindel. Pindel indel candidates wereobtained using the default parameters except the insertsize, which was provided by the CollectInsertSizeMetricsPicard subtool [48]. We only considered indels with Din-del quality scores greater than or equal to ten, which isequivalent to 90% confidence. All variants reported indbSNP (version 132) [49] and the 1,000 genomes project[21] were filtered out. Identifying the specific locationof an indel within a homopolymer or tandem repeat ischallenging, which effects the ability to label an indelproperly as somatic or germline. Therefore, we labeledindels as germline mutations if the distance of the repeat-ing region from the start position of the indel was longerthan the distance to the closest indel in the normal sam-ple. The repetitive sequences that were considered canbe any combination of base pairs between 1 and 6 bp inlength.Mutation ratesMutation rates within TFBSs were computed by divid-ing the number of SNVs or indels lying within TFBSs bythe total number of nucleotides within TFBSs (that is,76,160,599). The included TFBSs were predicted withinChIP-seq peak regions as described below. A similar com-putation was performed for protein-coding exons usingthe exonic start and end positions from RefSeq. Thetotal number of nucleotides covered by the exons is65,469,364.When computing local mutation rates, we consideredregions directly flanking TFBSs and exons. The flankingregions were obtained using the flank and subtract sub-commands of BEDTools [50] by extracting 1 kb upstreamand 1 kb downstream of the TFBSs (respectively, exons)and filtering out sequences overlapping TFBSs (respec-tively, exons).Genomes with randomly distributed mutations wereconstructed by shuffling all the mutations from cohort 1or cohort 2 in the human genome using the shuffle sub-command of BEDTools [50]. Then 1,000 genomes werecomputed for each cohort. For each genome, we calcu-lated the mutation rates in TFBSs and exons using therandomly positioned mutations.ChIP-seq data and transcription factor binding sitepredictionsWe collected 477 human TF ChIP-seq data sets fromboth ENCODE [24] and publications collected in PAZAR[51] with an associated TF binding profile described inthe JASPAR database [28] (Additional file 5). ChIP-seqpeak regions called in the studies were retrieved from thecorresponding analyses. TF binding profiles for the corre-sponding 103 TFs were retrieved from the 2014 release ofthe JASPAR database [28].TFBS predictions were obtained by scanning PWMsderived from the TF binding profiles (see [52] for theTF binding profile to PWM conversion) using the TFBSPerl module [53]. The PWMs were applied to the wholelength ChIP-seq peaks and we further considered in ouranalysis the hits for which the relative PWM score wasover 85% (see [52,53]). The default threshold of 85%Mathelier et al. Genome Biology (2015) 16:84 Page 14 of 17was used to call TFBSs as in previous studies [54,55].We predicted TFBSs covering 76,160,599 bp. Note thatboth strands on the reference and alternative genomeswere scanned with the PWMs to search for the optimalhits.Frequently targeted regionsFrequently targeted regions were obtained by sliding a 1-kb window along the human genome with steps of 500 bp.The makewindows subcommand of BEDTools [50] wasused to construct the set of window coordinates. For eachposition, we recorded the number of cis-regulatory muta-tions overlapping the window using the BEDTools [50]intersect subcommand. Only windows containing at leastthree mutations were considered and plotted in Figure 4.Mutations from the two cohorts were analyzed separately(Figure 4A,B) and combined (Figure 4C). For each muta-tion lying within a frequently targeted region, the closestgene was extracted using the TSS positions of the RefSeqgenes by applying the closest subcommand of BEDTools[50]. Genes with a TSS at a distance of at most 2 kb wereused for Figure 4.Prediction of alternative transcription factor binding sitesFor each mutation, we computed the best PWM scoreusing the alternative sequence containing the mutation.To compute the PWM score, we extracted sequenceswith a length of 2n \u00E2\u0088\u0092 1 bp (with n being the lengthof the considered PWM, Additional file 1: Figure S11A)centered around the SNV to identify regions that couldcontain a better TFBS at any overlapping position onthe alternative sequence (Additional file 1: Figure S11B).Alternative TFBS scores resulting from an insertionwere computed for sequences of length 2n \u00E2\u0088\u0092 2 + i bpwhere i represents the length of the insertion (Additionalfile 1: Figure S11C). Similarly, alternative TFBS scoresresulting from a deletion were obtained by scanning the2(n \u00E2\u0088\u0092 1) bp sequence centered at the deletion region(Additional file 1: Figure S11D). When scanning alter-native sequences with the PWMs, only the best hitper sequence was recorded. We considered a mutation(SNV or indel) to be deleterious for a TFBS if thePWM relative score was below 80% for the alternativesequence.MANTAAll the predicted TFBS positions can be scanned foroverlap with SNVs using our dedicated Mongo databasefor the analysis of TFBS alterations (MANTA). MANTAstores the positions of all predicted TFBSs as well as allthe potential SNVs overlapping these positions. For eachpotential SNV, you can retrieve information about the ref-erence and alternative best TFBSs along with their scores(see \u00E2\u0080\u0098Prediction of alternative transcription factor bindingsites\u00E2\u0080\u0099 for the computation of the alternative scores). TheMANTA source code is available at [56] and the systemcan be interrogated at [57].xseqxseq analyses were performed as follows.Preprocessing: Identifymutated genesMutations lying within TFBSs impacting gene expres-sion Mutations lying within TFBSs were obtained usingMANTA. The closest gene to each mutation was obtainedusing the set of TSSs of known refSeq genes from UCSC.When finding the closest gene, we consider the start andend positions for the mutation and the start positions forall protein-coding TSSs from refSeq. Only mutations lyingwithin TFBSs with a potential impact on gene expres-sion were considered. Namely, we require that the closestgene to the corresponding mutations to be either up- ordown-regulated in the sample of interest. To determine ifa gene is deregulated, we considered the distribution ofexpression of the gene in the cancer samples and requiredthat the expression in the corresponding sample was >\u00CE\u00BC + 1\u00CF\u0083 or < \u00CE\u00BC \u00E2\u0088\u0092 1\u00CF\u0083 where \u00CE\u00BC and \u00CF\u0083 represent the meanand standard deviation of the distribution of expressionvalues.Mutations lying within protein-coding exons Allmutations lying within a protein-coding exon were con-sidered. Namely, SNPeff [58] was used to extract themutations overlapping protein-coding regions and theirpredicted impact on the protein. Additional file 6 lists themutation impacts that were considered in the analysis.xseq analysisAll genes obtained from the previous step were used in theinput to the xseq tool, which is a probabilistic model thataims to encode the impact of somatic mutations on geneexpression profiles. The model uses a generative hierar-chical Bayes approach, which has as input three observedquantities: a patient-gene expression matrix, a patient-gene mutation matrix and a graph containing knowninteractions between genes (for example, from pathwaydatabases). The model has two key unobserved randomvariables, which constitute the output: Dg is a Bernoullirandom variable where Dg = 1 indicates that gene ginfluences expression when mutated; Fpg |Dg is a Bernoullirandom variable where Fpg indicates that mutated geneg influences expression in patient p. As such, we modelexpression influence at two levels: over the patient popu-lation and at the level of individual mutations in individualpatients. Random variables are estimated using the beliefpropagation algorithm, with outputs consisting of tworelevant probabilities: Pr(Dg) and Pr(Fpg ). Software imple-mented in an R package encoding xseq is available atMathelier et al. Genome Biology (2015) 16:84 Page 15 of 17[59,60]. By considering the disruption likelihoods of eachmutated gene and its neighbors in biological networks,xseq computes the probability of a mutated gene beingderegulated and causing cascading dysregulation effectson its neighbors.Post-processingxseq provides the probability of each input gene being adriver gene in the specific samples where it is mutated(single-sample probability, Pr(Fpg )) as well as its driverpotential considering all samples (all-samples probabil-ity, Pr(Dg)). Potential false positives from xseq are pro-duced when a gene is only mutated in a single samplebecause there is minimal information for calculating theall-samples probability properly. By plotting a histogramof all the calculated probabilities when considering allsamples (Additional file 1: Figure S12), a distinct peakwas observed, which is formed by a large number ofthese false positives. These distinct peaks were used asa threshold; genes must have an all-samples probabil-ity greater than or equal to 0.5 in cohort 1 and 0.8 incohort 2 to be considered in our analyses. Furthermore,we require that the single-sample probability of a geneis greater than or equal to 0.5 and the gene is predictedin at least two samples to be considered in the analysis(Figure 6).Circos plotsThe Circos plots in Figure 4 were drawn using the Circostool version 0.64 [61].Functional enrichment analysesFunctional enrichment analyses were performed using theEnrichr tool [31] (as of 20 January 2015) through its APIusing the poster library of Python2.7. A term is consid-ered to be enriched if the associated adjusted P \u00E2\u0089\u00A4 0.05.Visualization for the enrichment plots were constructedmanually using Cytoscape 3.1.0 [62]. Enrichment resultsassociated with Mus musculus in WikiPathways were fil-tered out and only Homo sapiens associated terms wereconserved.The functional enrichment analysis illustrated inFigure 5 was obtained from the list of genes provided inFigure 4 (Additional file 2). The functional enrichmentanalysis in Figure 7 was computed using the genes pre-dicted by xseq (Figure 6 and Section \u00E2\u0080\u0098xseq\u00E2\u0080\u0099) along withtheir biological network neighbors with altered expression(that is, predicted by xseq to have a higher probability ofbeing up- or down-regulated than being neutral) in cohort1, cohort 2 and their intersection (Additional file 2).Statistical analysesHypergeometrical P values were computed using the phy-per function of the R environment [63].Additional filesAdditional file 1: Additional figures. Supplementary figures referencedin the manuscript along with their descriptions.Additional file 2: Functional enrichments. Lists of genes used for thefunctional enrichment analyses along with the complete Enrichr resultswith a corrected P < 0.05.Additional file 3: Case examples. Details of the case examples predictedby xseq (HAS2, GNA13, BCL6 and ROBO1).Additional file 4: HGNC symbols. List of HGNC symbols associated withthe genes that were analyzed in this study.Additional file 5: ChIP-seq experiments. List of ChIP-seq experimentsused in this analysis along with the associated JASPAR TF binding profiles.Additional file 6: Mutation impacts. List of mutation impacts fromSNPeff that have been considered in this study.AbbreviationsBCL: B-cell lymphoma; bp: base pair; ChIP-seq: chromatin immunoprecipitationsequencing; kb: kilobase; PMBCL: primary mediastinal large B-cell lymphoma;PWM: position weight matrix; RNA-seq: RNA-sequencing; sd: standarddeviation; SNV: single nucleotide variant; TF: transcription factor; TFBS:transcription factor binding site; TSS: transcription start site.Competing interestsThe authors declare that they have no competing interests.Authors\u00E2\u0080\u0099 contributionsWWW and SPS were responsible for project conception and oversight. AMwas responsible for the analysis design and execution. AM, CL and AWZundertook the bioinformatics analysis. JD developed the xseq tool. AM, AWZand DJA developed the MANTA database. AM, WWW and SPS wrote themanuscript. All authors read and approved the final manuscript.AcknowledgementsWe thank Jonathan Lim for his help with the PAZAR database, Miroslav Hatasfor systems support and Dora Pak for management support. AM, DJA andWWW are supported by the Genome Canada Large Scale Applied ResearchGrant 174CDE. CL is supported by a Genome Canada/Genome BC ResearchGrant for Bioinformatics/Computational Biology. SPS is supported by a CanadaResearch Chair and is a Michael Smith Foundation for Health Research Scholar.Funding has been provided by the Child and Family Research Institute,Vancouver, to AM and AWZ and the BC Cancer Foundation to SPS.Author details1Department of Medical Genetics, Centre for Molecular Medicine andTherapeutics, Child and Family Research Institute, University of BritishColumbia, 950 West 28th Avenue, V5Z 4H4, Vancouver, BC, Canada.2Department of Molecular Oncology, British Columbia Cancer Agency, V5Z1L3 Vancouver, BC, Canada. 3Bioinformatics Graduate Program, University ofBritish Columbia, V5Z 1L3, Vancouver, BC, Canada. 4Department of ComputerScience, University of British Columbia, V6T 1Z4 Vancouver, BC, Canada.5Department of Pathology and Laboratory Medicine, University of BritishColumbia, G227-2211 Vancouver, BC, Canada.Received: 1 December 2014 Accepted: 7 April 2015References1. Kandoth C, McLellan MD, Vandin F, Ye K, Niu B, Lu C, et al. Mutationallandscape and significance across 12 major cancer types. Nature.2013;502:333\u00E2\u0080\u00939. doi:10.1038/nature12634.2. Lawrence MS, Stojanov P, Mermel CH, Robinson JT, Garraway LA, GolubTR, et al. Discovery and saturation analysis of cancer genes across 21tumour types. Nature. 2014;505:495\u00E2\u0080\u0093501.3. Mathelier A, Shi W, Wasserman WW. Identification of alteredcis-regulatory elements in human disease. Trends Genet. 2015;31:67\u00E2\u0080\u009376.Mathelier et al. Genome Biology (2015) 16:84 Page 16 of 174. Reijnen MJ, Sladek FM, Bertina RM, Reitsma PH. Disruption of a bindingsite for hepatocyte nuclear factor 4 results in hemophilia B Leyden. ProcNatl Acad Sci. 1992;89:6300\u00E2\u0080\u00933.5. Ludlow LB, Schick BP, Budarf ML, Driscoll DA, Zackai EH, Cohen A, et al.Identification of a mutation in a GATA binding site of the plateletglycoprotein ib\u00CE\u00B2 promoter resulting in the Bernard\u00E2\u0080\u0093Soulier syndrome.J Biol Chem. 1996;271:22076\u00E2\u0080\u009380. doi:10.1074/jbc.271.36.22076.6. Huang FW, Hodis E, Xu MJ, Kryukov GV, Chin L, Garraway LA, et al.Highly recurrent TERT promoter mutations in human melanoma. Science.2013;339:957\u00E2\u0080\u00939. doi:10.1126/science.1229259.7. Horn S, Figl A, Rachakonda PS, Fischer C, Sucker A, Gast A, et al. TERTpromoter mutations in familial and sporadic melanoma. Science.2013;339:959\u00E2\u0080\u009361. doi:10.1126/science.1230062.8. Fredriksson NJ, Ny L, Nilsson JA, Larsson E. Systematic analysis ofnoncoding somatic mutations and gene expression alterations across 14tumor types. Nat Genet. 2014;46:1258\u00E2\u0080\u009363.9. Kapeller J, Houghton LA, Monnikes H, Walstab J, Moller D, Bonisch H,et al. First evidence for an association of a functional variant in themicroRNA-510 target site of the serotonin receptor-type 3e gene withdiarrhea predominant irritable bowel syndrome. Hum Mol Genet.2008;17:2967\u00E2\u0080\u009377. doi:10.1093/hmg/ddn195.10. VanderMeer JE, Ahituv N. Cis-regulatory mutations are a genetic cause ofhuman limb malformations. Dev Dyn. 2011;240:920\u00E2\u0080\u009330.doi:10.1002/dvdy.22535.11. Spielmann M, Brancati F, Krawitz PM, Robinson P, Ibrahim DM, FrankeM, et al. Homeotic arm-to-leg transformation associated with genomicrearrangements at the PITX1 locus. Am J Hum Genet. 2012;91:629\u00E2\u0080\u009335.doi:10.1016/j.ajhg.2012.08.014.12. Khurana E, Fu Y, Colonna V, Mu XJ, Kang HM, Lappalainen T, et al.Integrative annotation of variants from 1092 humans: application tocancer genomics. Science. 2013;342:1235587.doi:10.1126/science.1235587.13. Weinhold N, Jacobsen A, Schultz N, Sander C, Lee W. Genome-wideanalysis of noncoding regulatory mutations in cancer. Nat Genet. 2014;11:1160\u00E2\u0080\u00935. doi:10.1038/ng.3101.14. Borneman AR, Gianoulis TA, Zhang ZD, Yu H, Rozowsky J, SeringhausMR, et al. Divergence of transcription factor binding sites across relatedyeast species. Science. 2007;317:815\u00E2\u0080\u00939. doi:10.1126/science.1140748.15. AndersenMC, Engstr\u00C3\u0083u\u00CB\u009Dm PG, Lithwick S, Arenillas D, Eriksson P, LenhardB, et al. In silico detection of sequence variations modifying transcriptionalregulation. PLoS Comput Biol. 2008;4:5. doi:10.1371/journal.pcbi.0040005.16. Kasowski M, Grubert F, Heffelfinger C, Hariharan M, Asabere A, WaszakSM, et al. Variation in transcription factor binding among humans.Science. 2010;328:232\u00E2\u0080\u00935. doi:10.1126/science.1183621.17. Spivakov M, Akhtar J, Kheradpour P, Beal K, Girardot C, Koscielny G,et al. Analysis of variation at transcription factor binding sites inDrosophila and humans. Genome Biol. 2012;13:49.18. Chen C-C, Xiao S, Xie D, Cao X, Song C-X, Wang T, et al. Understandingvariation in transcription factor binding by modeling transcription factorgenome-epigenome interactions. PLoS Comput Biol. 2013;9:1003367.doi:10.1371/journal.pcbi.1003367.19. Cusanovich DA, Pavlovic B, Pritchard JK, Gilad Y. The functionalconsequences of variation in transcription factor binding. PLoS Genet.2014;10:1004226. doi:10.1371/journal.pgen.1004226.20. Kheradpour P, Ernst J, Melnikov A, Rogov P, Wang L, Zhang X, et al.Systematic dissection of regulatory motifs in 2000 predicted humanenhancers using a massively parallel reporter assay. Genome Res. 2013;23:800\u00E2\u0080\u009311. doi:10.1101/gr.144899.112.21. McVean GA, Altshuler DM, Durbin RM, Abecasis GR, Bentley DR,Chakravarti A, et al. An integrated map of genetic variation from 1,092human genomes. Nature. 2012;491:56\u00E2\u0080\u009365. doi:10.1038/nature11632.22. Johnson DS, Mortazavi A, Myers RM, Wold B. Genome-wide mapping ofin vivo protein-DNA interactions. Science. 2007;316:1497\u00E2\u0080\u0093502.doi:10.1126/science.1141319.23. Sikora-Wohlfeld W, Ackermann M, Christodoulou EG, Singaravelu K,Beyer A. Assessing computational methods for transcription factor targetgene identification based on ChIP-seq data. PLoS Comput Biol. 2013;9:1003342. doi:10.1371/journal.pcbi.1003342.24. The ENCODE Project Consortium. The ENCODE (ENCyclopedia of DNAelements) project. Science. 2004;306:636\u00E2\u0080\u009340.doi:10.1126/science.1105136.25. Morin RD, Mungall K, Pleasance E, Mungall AJ, Goya R, Huff RD, et al.Mutational and structural analysis of diffuse large B-cell lymphoma usingwhole-genome sequencing. Blood. 2013;122:1256.doi:10.1182/blood-2013-02-483727.26. Richter J, Schlesner M, Hoffmann S, Kreuz M, Leich E, Burkhardt B, et al.Recurrent mutation of the ID3 gene in Burkitt lymphoma identified byintegrated genome, exome and transcriptome sequencing. Nat Genet.2012;44:1316\u00E2\u0080\u009320. doi:10.1038/ng.2469.27. Ding J, Bashashati A, Roth A, Oloumi A, Tse K, Zeng T, et al.Feature-based classifiers for somatic mutation detection intumour-normal paired sequencing data. Bioinformatics. 2012;28:167\u00E2\u0080\u009375.doi:10.1093/bioinformatics/btr629.28. Mathelier A, Zhao X, Zhang AW, Parcy F, Worsley-Hunt R, Arenillas DJ,et al. JASPAR 2014: an extensively expanded and updated open-accessdatabase of transcription factor binding profiles. Nucleic Acids Res.2013;42:142\u00E2\u0080\u00937. doi:10.1093/nar/gkt997.29. Forbes SA, Tang G, Bindal N, Bamford S, Dawson E, Cole C, et al.COSMIC (the catalogue of somatic mutations in cancer): a resource toinvestigate acquired mutations in human cancer. Nucleic Acids Res.2010;38:652\u00E2\u0080\u00937. doi:10.1093/nar/gkp995.30. Khodabakhshi AH, Morin RD, Fejes AP, Mungall AJ, Mungall KL,Bolger-Munro M, et al. Recurrent targets of aberrant somatichypermutation in lymphoma. Oncotarget. 2012;3:1308\u00E2\u0080\u009319.31. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al. Enrichr:interactive and collaborative HTML5 gene list enrichment analysis tool.BMC Bioinformatics. 2013;14:128.32. Ding J, McConechy MK, Horlings HM, Ha G, Fong CC, Funnell T, et al.Systematic analysis of somatic mutations impacting gene expression intwelve tumor types. Submitted.33. Pelengaris S, Khan M, Evan G. c-MYC: more than just a matter of life anddeath. Nat Rev Cancer. 2002;2:764\u00E2\u0080\u009376. doi:10.1038/nrc904.34. Bonner AE, Lemon WJ, You M. Gene expression signatures identify novelregulatory pathways during murine lung development: implications forlung tumorigenesis. J Med Genet. 2003;40:408\u00E2\u0080\u009317.35. Xian J, Aitchison A, Bobrow L, Corbett G, Pannell R, Rabbitts T, et al.Targeted disruption of the 3p12 gene, DUTT1/ROBO1, predisposes miceto lung adenocarcinomas and lymphomas with methylation of the genepromoter. Cancer Res. 2004;64:6432\u00E2\u0080\u00937.doi:10.1158/0008-5472.CAN-04-2561.36. Parray A, Siddique HR, Kuriger JK, Mishra SK, Rhim JS, Nelson HH, et al.ROBO1, a tumor suppressor and critical molecular barrier for localizedtumor cells to acquire invasive phenotype: study in African\u00E2\u0080\u0093Americanand Caucasian prostate cancer models. Int J Cancer. J Int Du Cancer.2014;135:2493\u00E2\u0080\u0093506. doi:10.1002/ijc.28919.37. Orton RJ, Sturm OE, Gormand A, Wolch W, Gilbert DR. Computationalmodelling reveals feedback redundancy within the epidermal growthfactor receptor/extracellular-signal regulated kinase signalling pathway.IET Syst Biol. 2008;2:173\u00E2\u0080\u009383. doi:10.1049/iet-syb:20070066.38. Walker SR, Nelson EA, Yeh JE, Pinello L, Yuan G-C, Frank DA, et al. STAT5outcompetes STAT3 to regulate the expression of the oncogenictranscriptional modulator BCL6. Mol Cell Biol. 2013;33:2879\u00E2\u0080\u009390.doi:10.1128/MCB.01620-12.39. UCSC hgTables. http://genome.ucsc.edu/cgi-bin/hgTables.40. Andersson R, Gebhard C, Miguel-Escalada I, Hoof I, Bornholdt J, Boyd M,et al. An atlas of active enhancers across human cell types and tissues.Nature. 2014;507:455\u00E2\u0080\u009361. doi:10.1038/nature12787.41. ICGC Malignant Lymphoma webpage. https://dcc.icgc.org/projects/MALY-DE.42. Lawrence M, Huber W, Pag\u00C3\u0083l\u00C2\u00B4s H, Aboyoun P, Carlson M, Gentleman R,et al. Software for computing and annotating genomic ranges. PLoSComput Biol. 2013;9:1003118. doi:10.1371/journal.pcbi.1003118.43. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S,et al. Bioconductor: open software development for computationalbiology and bioinformatics. Genome Biol. 2004;5:80.doi:10.1186/gb-2004-5-10-r80.44. Gray KA, Daugherty LC, Gordon SM, Seal RL, Wright MW, Bruford EA,et al. Genenames.org: the HGNC resources in 2013. Nucleic Acids Res.2013;41:545\u00E2\u0080\u009352. doi:10.1093/nar/gks1066.45. MutationSeq webpage. http://compbio.bccrc.ca/software/mutationseq/.46. Albers CA, Lunter G, MacArthur DG, McVean G, Ouwehand WH, DurbinR, et al. Dindel: Accurate indel calls from short-read data. Genome Res.2011;21:961\u00E2\u0080\u009373. doi:10.1101/gr.112326.110.Mathelier et al. Genome Biology (2015) 16:84 Page 17 of 1747. Ye K, Schulz MH, Long Q, Apweiler R, Ning Z. Pindel: a pattern growthapproach to detect break points of large deletions and medium sizedinsertions from paired-end short reads. Bioinformatics. 2009;25:2865\u00E2\u0080\u009371.doi:10.1093/bioinformatics/btp394.48. Picard webpage. http://broadinstitute.github.io/picard.49. Sherry ST, Ward M-H, Kholodov M, Baker J, Phan L, Smigielski EM,Sirotkin K, et al. dbSNP: the NCBI database of genetic variation. NucleicAcids Res. 2001;29:308\u00E2\u0080\u009311. doi:10.1093/nar/29.1.308.50. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparinggenomic features. Bioinformatics. 2010;26:841\u00E2\u0080\u009342.doi:10.1093/bioinformatics/btq033.51. Portales-Casamar E, Kirov S, Lim J, Lithwick S, Swanson MI, Ticoll A, et al.PAZAR: a framework for collection and dissemination of cis-regulatorysequence annotation. Genome Biol. 2007;8:207.52. Wasserman WW, Sandelin A. Applied bioinformatics for the identificationof regulatory elements. Nat Rev Genet. 2004;5:276\u00E2\u0080\u009387.doi:10.1038/nrg1315.53. Lenhard B, Wasserman WW. TFBS: Computational framework fortranscription factor binding site analysis. Bioinformatics. 2002;18:1135\u00E2\u0080\u00936.doi:10.1093/bioinformatics/18.8.1135.54. Kwon AT, Arenillas DJ, Hunt RW, Wasserman WW. oPOSSUM-3:Advanced analysis of regulatory motif over-representation across genesor ChIP-seq datasets. G3; Genes|Genomes|Genet. 2012;2:987\u00E2\u0080\u00931002.doi:10.1534/g3.112.003202.55. Lenhard B, Sandelin A, Mendoza L, Engstr\u00C3\u0083u\u00CB\u009Dm P, Jareborg N,Wasserman WW, et al. Identification of conserved regulatory elements bycomparative genome analysis. J Biol. 2003;2:13.56. MANTA source code. https://github.com/wassermanlab/manta.57. MANTA webpage. http://manta.cmmt.ubc.ca/manta/.58. Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, et al. Aprogram for annotating and predicting the effects of single nucleotidepolymorphisms, SnpEff: SNPs in the genome of Drosophila melanogasterstrain w1118; iso-2; iso-3. Fly. 2012;6:80\u00E2\u0080\u009392. doi:10.4161/fly.19695.59. xseq webpage. http://compbio.bccrc.ca/software/xseq/.60. xseq source code. https://bitbucket.org/MO_BCCRC/xseq.61. Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, et al.Circos: an information aesthetic for comparative genomics. Genome Res.2009;19:1639\u00E2\u0080\u009345. doi:10.1101/gr.092759.109.62. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al.Cytoscape: a software environment for integrated models ofbiomolecular interaction networks. Genome Res. 2003;13:2498\u00E2\u0080\u0093504.doi:10.1101/gr.1239303.63. R webpage. http://www.R-project.org.Submit your next manuscript to BioMed Centraland take full advantage of: \u00E2\u0080\u00A2 Convenient online submission\u00E2\u0080\u00A2 Thorough peer review\u00E2\u0080\u00A2 No space constraints or color figure charges\u00E2\u0080\u00A2 Immediate publication on acceptance\u00E2\u0080\u00A2 Inclusion in PubMed, CAS, Scopus and Google Scholar\u00E2\u0080\u00A2 Research which is freely available for redistributionSubmit your manuscript at www.biomedcentral.com/submit"@en . "Article"@en . "10.14288/1.0228385"@en . "eng"@en . "Reviewed"@en . "Vancouver : University of British Columbia Library"@en . "BioMed Central"@en . "10.1186/s13059-015-0648-7"@en . "Attribution 4.0 International (CC BY 4.0)"@en . "http://creativecommons.org/licenses/by/4.0/"@en . "Faculty"@en . "Cis-regulatory somatic mutations and gene-expression alteration in B-cell lymphomas"@en . "Text"@en . "http://hdl.handle.net/2429/56170"@en .