Variants in potential regulatory elements

Regulatory elements, especially promoters and enhancers, are capable of regulating the expression of specific genes. We collected functional annotations from ENCODE [15] (transcription factor binding sites and the high-resolution motifs within them, enhancers from genome segmentations and DNase I hypersensitive sites) and regions that are highly occupied by transcription factors (HOT) from Yip et al. [27] to annotate variants in potential regulatory sequences. We evaluated their functional effect from sequence and network levels in the following sections.

Nucleotide-level impact of regulatory variants

Regulatory mutations can cause transcriptional alterations by either loss-of or gain-of- function effects. Loss-of-function events in transcription factor binding motifs are likely to cause deleterious impact [26,28,29]. Variants decreasing the position weight matrix (PWM) scores could potentially alter the binding strength of transcription factors, or even eliminate the binding. Our framework consists of a module to detect motif-breaking events - defined as variants decreasing PWMs (Material and methods). Meanwhile, gain of new binding sites caused by somatic mutations can constitute driver events [7–10]. To the best of our knowledge, there is no automated tool to detect such events in whole tumor genomes. We incorporated a gain-of-motif scheme to scan and statistically evaluate [30] all possible motifs created by mutated alleles in promoter or enhancer regions. For each variant (SNV or indel), we concatenated it with +/- 29 bp nucleotides around it and calculated sequence scores for each possible motif against the PWMs. Gain-of-motif events are identified when the sequence score with mutated allele is significantly higher than the background (P <4e-8), whereas that with germline allele is not. As discussed later, our scheme is validated by the detection of motif-gaining events caused by the two driver mutations in TERT promoter (Additional file 1: Table S1).

Variants in conserved regions

Sequences that are under strong negative or purifying selection are thought to have important biological functions [31]. In previous studies, oncogenes or tumor suppressor genes are found to experience higher intensity of selective pressure than other disease-related and non-disease genes [32]. Cross-species genome comparison is a powerful approach to identify evolutionary conserved sequences. For example, GERP [33] is developed to estimate the position-specific evolutionary constraint; sequences that are shared across species are defined as ultra-conserved elements [31]. Meanwhile, human population-level constrained regions are identified from 1000 Genomes [26,34] using depletion of common polymorphisms. We combined these data to detect potential deleterious variants in noncoding sequences. Each variant will be annotated with its corresponding conservation information.

Linking regulatory elements with likely target genes

Functional genomic studies have characterized biological functions of a large number of genes. Linking regulatory variants with coding genes, especially those well-known cancer driver genes, will help us understand the regulatory mechanisms that govern oncogenesis and potentially benefit therapeutic treatment. Positioned distant to their target genes, regulatory elements regulate gene expression through long-range interactions [35]. The linkages between regulatory elements and genes remain elusive. To explore likely functional consequences of regulatory variants, we comprehensively define regulatory element-target gene pairs through correlating various epigenetic modifications with gene expression levels. We consider the enhancer marks H3K4me1 and H3K27ac as two types of activity signals, and DNA methylation as an inactivity signal. Using ChIP-Seq and RNA-Seq data from the Roadmap Epigenomics Mapping Consortium (REMC), for each regulatory element-candidate target gene pair, we computed the correlations of H3K4me1 and H3K27ac and the anti-correlations of DNA methylation at the regulatory element with gene expression levels across 20 tissue types (Material and methods). In total, we identified approximately 2,448 K significant associations involving roughly 1,279 K regulatory elements and 17 K genes (Material and methods; Additional file 1: Figure S1). All noncoding variants in these regulatory elements can be associated with potential target genes (with various association confidences). To incorporate the ever-increasing amounts of genomic data, we include a pipeline for users to extend the data context with their own data. For example, users can input annotation regions or chromatin marks to find novel associations between regulatory elements and coding genes (Additional file 1).

Network analysis of variants associated with genes

Unlike germline mutations, somatic alterations are not expected to be under organism-level evolutionary selection pressure and thus are more likely to affect functional centers in gene interaction networks [36]. Network studies have found that cancer genes possess high topological centralities, even higher than essential genes [26,36]. We used the regulatory element-target gene pairs to connect noncoding variants to a variety of networks: protein-protein interaction, regulatory and phosphorylation networks [26,35,37]. For each noncoding variant, we calculated the scaled network centrality (the percentile after ordering centralities of all genes in a particular network) of the associated gene in each network (Material and methods). Among the different network centralities, we used the maximum centrality as the network disruptive measure of the variant. The higher this value, the more likely the variant is to be deleterious. We make the scheme flexible so it can integrate user networks in addition to the pre-collected networks.

Recurrent elements across cancer samples

One criterion to identify cancer driver genes is to examine their mutational recurrence across multiple samples [3]. We extended the concept to noncoding regulatory elements, such as transcription factor binding sites. Our method is able to detect sites, genes, and regulatory elements that are mutated in at least two samples from user-input (Material and methods).

When having small sample sizes, comparisons with available tumor genomes would be useful to discover recurrent mutations. Related to this, we have created a recurrence database (Recurrence DB) (including regulatory elements, coding genes and the same-site mutations) with publicly available cancer whole-genome sequencing data. Currently, we have identified recurrent loci or sites from 570 samples of 10 tumor types [38–40] and from COSMIC [41] (Table 1; Material and methods). Variants in user-input tumor genomes are compared to the recurrence database and the results in different cancer types are reported in the output. The use of the database along with our framework would provide higher confidence in prioritization of regulatory drivers (Figure 2). The database will be updated with newly available dataset.

Weighted scoring scheme to prioritize variants

All of the above features are used to annotate and score variants (Figure 3). To integrate the various features to predict ‘high-impact’ somatic alterations, we developed a weighted scoring scheme, taking into account the relative importance of each feature (Figure 2). In general, features can be classified into two classes: discrete and continuous. Discrete features are binary, such as in ultra-conserved elements or not. For continuous features, taking the ‘motif-breaking score’ as an example, the values would be the changes in PWMs. We weighted each feature based on the mutation patterns observed in natural polymorphisms (Material and methods; Additional file 1: Table S3). Constrained by selective pressure, natural variations tend to arise in functionally unimportant regions. Thus, features that are frequently observed are less likely to contribute to the deleteriousness of variants and are weighted less. We calculated the information content to denote the importance of each feature. For each cancer mutation, we scored it by summing the information contents of all its features (details in Material and methods). Variants ranked on top of the output are those with higher scores and are most likely to be deleterious.

Open in a separate window

Figure 3

Weighted scoring scheme. (A) Features used in the weighted scoring scheme. Features can be classified into two classes: discrete and continuous. Discrete features are binary, such as in ultra-conserved elements or not. For continuous features, taking the ‘motif-breaking score’ as an example, the values would be the changes in PWMs. * - the feature depends on the user (see Figure 2); (B) We weighted each feature based on the mutation patterns observed in natural polymorphisms. Features that are frequently observed are less likely to contribute to the deleteriousness of variants and are weighted less (entropy-based method, details described in Material and methods). For a continuous feature, such as the ‘motif-breaking scores’, we calculated a weight for each observed value. The x-axis is the observed motif-breaking score and the y-axis is the corresponding weight. The black line shows the values observed in natural polymorphisms. We then fit a smooth curve (the red dashed line) to obtain continuous weights for all possible motif-breaking scores.

Highlighting variants using prior knowledge of genes and user annotations

Interpretation of the functional impact of noncoding variants can be greatly enhanced if the function of its target protein-coding gene is known. Many cancer genes are known to play a crucial role in cell proliferation and DNA repair. We integrated prior knowledge of genes, such as known cancer-driver genes [2,42], genes involved in DNA repair [43] and actionable genes (‘druggable’ genes) [44] to highlight noncoding variants that are likely to be involved in cancer development and growth or their associated genes could be used as drug targets. In addition, user-specific gene lists can be easily input (Figure 1 and Additional file 1).

Variants in regulatory elements may disrupt the expression of coding genes. Differential expression of target genes in cancer samples indicates the potential effect of noncoding variants. We provide a ‘differential gene expression analysis’ module (Figure 1 and Material and methods) to detect differentially expressed genes in cancer samples (relative to matched normal) from RNA-Seq data. Lists of differentially expressed genes can be generated and used to annotate variants.

We also provide an option for users to incorporate their own annotations. Impact of variants in regulatory regions is generally restricted to cis-acting effects that control the spatial and temporal patterns of gene expression [13]. Activation of regulatory elements is associated with the underlying epigenetic or open chromatin landscape, which is largely cell-line specific [45]. For example, enrichment of H3K27ac may indicate an active state of enhancers in a particular sample. Therefore, it would be useful to incorporate sample-specific epigenetic or open chromatin profiles, if available, to highlight variants in activated or inactivated regulatory sequences.

All features used in our method and corresponding details are described in Additional file 1: Table S2.

Performance on regulatory cancer somatic variants and germline pathogenic variants

Recurrent somatic variants

Currently, only two known regulatory variants are thought to act as cancer drivers. Hence, to evaluate the performance of our scoring scheme, we used recurrence to approximate the deleteriousness of somatic variants. Recurrence is considered as one potential sign of positive selection among tumors and is more likely to be associated with driver events [3]. We examined recurrence from two perspectives: recurrence at the exact same site and recurrence in the same regulatory element. First, we classified COSMIC regulatory somatic variants [41] as same-site recurrent or non-recurrent (Material and methods) [25]. Our method scores recurrent variants higher than non-recurrent ones (Wilcoxon rank-sum test: P value <2.2 e-16; Figure 4A). Variants that occur in more than two samples have higher scores than those that are in two samples. Next we evaluated variants in recurrent regulatory elements using a separate dataset. We ran our pipeline on 119 breast cancer samples [38] and classified variants as occurring in recurrent elements or not (Material and methods). We found that variants in recurrent elements get significantly higher scores (Wilcoxon rank-sum test: P value <2.2e-16; Figure 4B) than variants elsewhere. Similar patterns are observed with other cancer types (Additional file 1: Figure S5). Results from CADD and GWAVA are shown in Additional file 1: Figure S6.

Open in a separate window

Figure 4

Application to pathogenic germline and cancer somatic regulatory variants. (A) Score distribution of variants based on their recurrence in COSMIC. Variants are classified into three categories based on the number of recurrent samples; (B) Score distribution of variants based on recurrent regulatory elements in 119 breast cancer samples. Recurrent regulatory elements are defined as elements mutated in more than one sample, such as the same DNase I hypersensitive site; (C) Prediction scores for regulatory variants in HGMD and controls; (D) ROC curves comparing HGMD variants with controls. All scores in this figure are ‘core scores’ (Figure 2).

We note that cancer is a very heterogeneous disease and distinct molecular subtypes may involve unique oncogenic mechanisms. Thus, tumor samples from different patients may involve different driver events. These unique drivers would not show recurrence across samples. Furthermore, in the absence of large sample sizes, it might be impossible to detect recurrence of mutations. Our method would be especially useful in such scenarios, since it has the ability to prioritize deleterious variants in each tumor genome. Moreover, the functional relevance and hence the biological mechanism by which drivers act is largely unknown. Our method provides an in-depth annotation of such variants, including the relative contribution of each feature to its deleteriousness. This knowledge would greatly help understand the potential oncogenic mechanisms.

A case study: somatic variants from an individual tumor genome

High recurrence of the TERT promoter mutations indicates their important roles in tumorigenesis [7]. Among the 570 cancer samples we collected, seven samples contain the TERT promoter mutation (chr5: 1295228). We used one Medulloblastoma sample as an example to prioritize regulatory variants from whole-genome sequencing. Among the 2,183 somatic single nucleotide variants, the TERT promoter mutation ranks second (0.09%). Our method further suggests potential functional impact of this variant. As shown in Table 2, this mutation occurs in ENCODE regulatory regions, creates a novel ETS binding motif and potentially affects a highly connected and known cancer gene - TERT. It is also found in another five liver samples and 54 COSMIC samples in our recurrence database. Besides DNA sequences, epigenome or transcriptome could also be altered in cancer genomes. These data provide sample-specific activation or inactivation signatures of regulatory sequences. If provided, our framework is flexible in integrating those data into our annotation scheme (refer to Additional file 1 for details).

We also applied CADD and GWAVA on the Medulloblastoma sample. CADD ranked the TERT promoter mutation 224th (10.3%) and GWAVA ranked it 25th (1.15%) with the matched region model (Additional file 1: Figure S10). However, only our method shows that the mutation corresponds to gain-of an ETS binding motif in the promoter of a cancer driver gene. Results of the additional six samples are shown in Additional file 1: Table S4.

High-impact variants in motifs: nucleotide resolution effect

User-input variants are first filtered against natural polymorphisms based on user-defined minor allele frequency (MAF) threshold to get rid of unlikely somatic variants (hg19). Currently, SNVs and small indels (<=20 bp) will be analyzed.

Associating regulatory elements with likely target genes

We define both proximal and distal associations. For proximal associations, we assign variants in gene promoters, introns, or UTRs to their nearby genes. For distal associations, in addition to those identified in [27], we further expanded the method to all ENCODE regulatory elements and identified roughly 2,448 K significant associations between 1,279 K regulatory elements and 17 K genes (see below). The distributions of regulatory element-gene associations are shown in Additional file 1: Figure S1. The median number of associations is 26 per gene and 1 per regulatory element, respectively. The association confidence is reported in the output for each regulatory element - target gene pair.

Correlating histone modifications with gene-expression data to identify likely target genes of regulatory elements

Differential gene expression analysis

We incorporated a module to detect differentially expressed genes in cancer samples (relative to matched normal) from RNA-Seq data. When provided with gene expression files, our module calls NOISeq [54] when having RPKM values and DESeq [55] with raw read counts (from reads-mapping tools) to detect differentially expressed genes. Genes that are up- or down-regulated with FDR <0.05 (with biological replicates) and FDR <0.1 (without replicates) in cancer samples are identified and annotated in the output.

Network analysis of variants associated with genes

For each variant associated with genes, we examined its topological properties in various networks. For each network, we calculated the centrality position (cumulative probability after ordering centralities of all genes increasingly) of the associated gene. If one variant is associated with multiple genes or the associated gene participates in multiple networks, the maximum cumulative probability is used as the continuous value for centrality score. Scripts are provided to calculate network centralities (Additional file 1). Users can easily incorporate other networks in this analysis.

Recurrence database from whole-genome sequencing

With the increasing number of cancer samples being whole-genome sequenced, we are able to study recurrence patterns in regulatory sequences. We analyzed somatic alterations from 570 samples of 10 cancer types to create the recurrence database [38–40], similar to the cancer recurrent gene database in cBio [56]. For each cancer type, recurrent regulatory elements, coding genes, and the same-site mutations are stored as entries in the database. We also incorporated same-site recurrent regulatory variants from COSMIC (version 68) into our database. Recurrent elements are defined as identified in whole-genome sequencing and observed in at least two samples.

Weighted scoring scheme

Noncoding scoring scheme (weighted scoring scheme)

Features used to score noncoding variants are shown in Additional file 1: Table S3. In general, features can be classified into two classes: discrete and continuous. Discrete features are binary, such as in ultra-conserved elements or not. Continuous features: (1) GERP score; (2) motif-breaking score is the difference between germline and mutated alleles in PWMs; (3) motif-gaining score is the sequence score difference between mutated and germline alleles; (4) network centrality score (the cumulative probability, see ‘Network analysis of variants associated with genes’). If one variant has multiple values of a particular feature (for example, breaking multiple motifs), the largest value is used.

We weighted each feature based on the mutation patterns observed in the 1000 Genomes polymorphisms. We randomly selected 10% of the 1000 Genomes Phase 1 SNPs (approximately 3.7 M) and ran them through our pipeline. For each discrete feature d, we calculated the probability p_d that overlaps a natural polymorphism. Then we computed 1-Shannon entropy (1) as its weighted value w_d. The value ranges from 0 to 1 and is monotonically decreasing when the probability is between 0 and 0.5c below 0.5).

1

The situation is more complex for continuous features, as different feature values have different probabilities of being observed in natural polymorphisms. Thus, one weight cannot suffice for varied feature values. For a continuous feature c, which is associated with a score v_c (for example, motif-breaking score), we calculated feature weights for each v_c. In particular, we discretized at each v_c and computed 1-Shannon entropy using (2). Then we fitted a smooth curve for all v_c to obtain continuous . Now, when we come to evaluate the continuous feature c for a particular variant, we calculate its weighted value (on the curve) using the actual v_c corresponding to the variant.

2

Taking ‘motif-breaking score’ as an example (Figure 3B), for each score v, we calculated the probability of observing motif-breaking scores ≥v in polymorphism data, then used (2) to fit the smooth function. ‘nls’ function in R is used to fit curves.

The criterion of ‘GERP >2’ has been commonly used to define conserved bases [15]. For the GERP score, we chose to use a sigmoid transformation to transform scores to the range 0 to 1. The parameters we chose make the sigmoid curve sharp at ‘GERP = 2’ (Additional file 1: Figure S2). The sigmoid transformation preserves the ‘GERP > 2’ cutoff and makes the score continuous at the same time. We calculated (1) treating ‘GERP >2’ as a discrete feature. Then we used w_d*sigmoid transformed value to assign weighted value for each continuous GERP score.

Finally, for each cancer variant, we scored it by summing the weighted values of all its features (3). If a particular feature is not observed, it is not used in the scoring. Considering the situation that some features are subsets of other features, to avoid overweighting similar features, we took into account feature dependencies when calculating the summed scores. As shown in Additional file 1: Table S3, when having leaf features, the weighted values of root features are ignored. For example, when a variant occurs in sensitive regions, the score of ‘in functional annotations’ is not used in the sum-up. Leaf features are assumed independent. Variants ranked on top of the output are those with higher scores and are most likely to be deleterious.

3

Application to regulatory pathogenic and cancer somatic variants

All scores in this section are ‘core scores’ described in Figure 2.

This work was supported by the National Institutes of Health, AL Williams Professorship, and in part by the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center (Grant Number RR029676-01). JB acknowledges support from the National Science Foundation - Graduate Research Fellowship Program (Grant Number 1346837). Funding for open access charge: National Institutes of Health. We also thank Zhixiang Lin, Curran Oi, and Yanwei Cai for useful discussion and proofreading the manuscript.