Functional annotation of putative regulatory elements at cancer susceptibility Loci.

Most cancer-associated genetic variants identified from genome-wide association studies (GWAS) do not obviously change protein structure, leading to the hypothesis that the associations are attributable to regulatory polymorphisms. Translating genetic associations into mechanistic insights can be facilitated by knowledge of the causal regulatory variant (or variants) responsible for the statistical signal. Experimental validation of candidate functional variants is onerous, making bioinformatic approaches necessary to prioritize candidates for laboratory analysis. Thus, a systematic approach for recognizing functional (and, therefore, likely causal) variants in noncoding regions is an important step toward interpreting cancer risk loci. This review provides a detailed introduction to current regulatory variant annotations, followed by an overview of how to leverage these resources to prioritize candidate functional polymorphisms in regulatory regions.

Introduction

Genome-wide association studies (GWAS) have successfully identified hundreds of cancer susceptibility loci, and have been particularly powerful in the discovery of distal regulatory elements such as enhancers. Understanding the mechanisms underlying GWAS findings requires the functional characterization of causal alleles and is therefore an important next step for gaining deeper insight into cancer susceptibility. If the single-nucleotide polymorphism (SNP) with the strongest statistical association (ie, the “index SNP”) was guaranteed to be causal, then functional validation would likely occur more often. However, SNPs that are represented on genotyping platforms are explicitly chosen as markers for “bins” of SNPs having highly correlated genotypes, with the consequence that SNPs within the same bin will have similar levels of statistical association, even when the association is driven by a single functional variant. The present generation of genotyping platforms contain roughly one million markers that tag a total of nearly six million SNPs, for a mean bin size of approximately six SNPs. The median number of SNPs within a bin is somewhat smaller, but still greater than three. Under the assumption that most associated bins contain only one functional variant, a large proportion of index SNPs with reported associations are not functional (ie, they are not causally linked to the phenotype). Thus, winnowing the list of associated variants to those with hypothesized function is an important step prior to laboratory-based validation.

The large number of GWAS loci in noncoding regions1,2 has illuminated two important facets of genetic influences on sporadic cancer: (1) many of the loci contributing to cancer susceptibility are likely to impact gene regulatory elements, and (2) the individual effects of common risk variants are associated with small quantitative differences in expression rather than qualitative changes in protein structure.3 Given these findings, it is not surprising that the rapid discovery of genetic associations with cancer outcomes has not been easily translated into improved clinical practice. In part, the identification of causal variants within known loci will likely improve genetic risk prediction and thereby enable more tailored prevention strategies. Furthermore, improved understanding of the mechanisms underlying genetic associations might enable development of interventions that mimic the effects of the protective allele. As noted, testing variants for functional follow-up in the laboratory can be labor- and resource-intensive, making the selection of candidate variants a critical next step.

When studying familial cancer, the vast majority of findings were germline mutations that disrupt the coding sequence of tumor suppressor genes, DNA repair genes, or oncogenes. In protein-coding regions, there are many annotation resources for assigning functional information. For example, PMD,4 PhenCode,5 and HumDiv/HumVar/Polyphen26 are catalogs constructed from one or more variant databases or from the literature with the intent of annotating specific variants for functional effect. Similarly, nonsynonymous SNPs (nsSNPs) can be assigned function based on their impact on structural stability through databases like ProTherm7 and ProNit.8 In contrast, noncoding SNP annotations are constrained by a dearth of knowledge pertaining to the phenotypic impact of regulatory variation. However, the number of genome annotations relevant to regulation is increasing rapidly, and is finally maturing to the point that appropriate annotations can be leveraged by researchers to enrich lists of variation for candidates to test in the laboratory.

Most regulatory annotations attempt to define, or segment, noncoding regions into short stretches of the genome that associate with biochemical properties exhibited by a small number of experimentally validated regulatory elements. These segments are defined by: (1) biofeatures such as epigenetic remodeling of chromatin, (2) particular sequence motifs like conserved microRNA binding sites, (3) gene expression data, and (4) evolutionary constraint (see Fig. 1). This review will cover the existing tools that are available to annotate cancer loci with regulatory features (see Table 1).

Figure 1

Tools for bioinformatic annotation.

Note: *Combined evidence across histone modification, open chromatin, ChIP-seq protein and motif annotations provides finer demarcation of functional elements and stronger evidence for the regulatory potential.

Table 1

Description of bioinformatics tools used for functional follow-up of noncoding regions.

DATASET	GENOMIC CLASS	DESCRIPTION	DATA SOURCE/PROGRAM
1	Non-synonymous coding	Exonic positions leading to amino acid replacement	Nonsense-Ensembl,61 Missense-Polypen2,13 GERP++,12 SIFT,74 SnpEff,75 LS-SNP,20 FoldX76
2	Promoter	1kb regions upstream of annotated transcription start sites	RefSeq, UCSC Genome Browser50: ENCODE Histone Modifications9 and Uniformly processed signal from Roadmap Project tracks
3	TFBS	Transcription factor binding sites (TFBS) predicted in promoter & non-promoter regulatory elements	UCSC Genome Browser50: ChIP-seq Transcription factor,77 PWM-scana JASPAR,17 CONSITE,18 HaploReg19
4	Non-coding RNA	All types of experimentally supported non-coding RNA, including microRNAs	RNAdb 2.078 & miRBase 17.079
5	MicroRNA target site	Computationally predicted microRNA target sites within 3′ UTRs	TargetScanS 5.280
6	Enhancer element	Experimentally supported enhancer elements in any tissue	VISTA Enhancer Browser UCSC Table Browser81: ENCODE Histone Modifications9
7	Candidate non-specific regulatory element	Open chromatin loci in at least one human cell type, as assessed by DNase I hypersensitivity (DHS) mapping	UCSC Table Browser81: Duke and UW DNase I HS data from > 50 cell types77
8	Insulator elements	CTCF binding sites assessed by ChIP-seq technology	UCSC Table Browser81: ChIP-seq TFBS77
9	eQTL	Allele-specific differences in expression levels	GTEx eQTL Browser,82 TCGA51
10	Conserved Element		UCSC Table Browser81: PhastCons 46-way conservation83
11	Splice Site		BDGP

Note:

PWM-scan was applied using positional weight matrices (PWMs) from the Transfac database.

Defining Regulatory Elements

There are now many biological datasets that have been made available through the extensive efforts of the Encyclopedia of DNA Elements (ENCODE)9 and the NIH Roadmap Epigenomics10 projects. These datasets provide information that enables functional annotation of noncoding regions harboring association signals.11 The ENCODE and Roadmap projects were created with the intention of cataloging the genome-wide regulatory landscape across many cell types9 and tissues.10 These experimental datasets can be summarized into three general biofeatures: (1) chromatin structure as detected through DNaseI hypersensitivity (DHS) assays, (2) histone modifications detected through chromatin immunoprecipitation sequencing (ChIP-seq), and (3) protein-binding detected through ChIP-seq. In contrast to tissue-specific features, we can use known binding motifs to predict variants that disrupt poly-A signals or binding of post-transcriptional regulators like miRNAs. Comparison of sequence alignment across species can also be used to identify conserved sequences reflective of regions that are important to organismal fitness.12 PhastCons13 is a particularly useful conservation measurement for detection of regulatory elements as it factors in the constraint of neighboring nucleotides. Taken together, these datasets can identify specific remodeling and general accessibility of chromatin for the binding of regulatory machinery (eg, transcription factors [TFs]) or for other post-transcriptional mechanisms.

Since regulatory elements such as enhancers often facilitate cell-type specific expression,14 it is helpful to look for evidence in a variety of tissues or cell lines in addition to those specifically related to the cancer of interest. Examination of multiple cell types allows differentiation between tissue-specific regulatory mechanisms likely to contribute to cancer susceptibility and regulatory elements shared by many cell types. Epigenetic chromatin modifications are useful for identifying cell-specific regulatory elements. The methylation and acetylation of histone proteins changes chromatin accessibility for transcription, and such marks can serve as a powerful tool for identifying enhancers, promoters, and other regulatory regions like insulators. There are three histone modification marks that are particularly informative for the identification of most active enhancer and promoter regulatory regions, namely H3k4me1, H3k27ac, and H3k4me3. The H3k4me1 histone mark is associated with enhancers downstream of transcription start sites, and the H3k27ac signal is similarly thought to enhance transcription. Alternatively, the H3k4me3 mark is associated with active promoters. In addition to these three epigenetic signals, H3k27me3 is useful to distinguish between active and poised promoters/enhancers.

There are just over 350 cell lines assayed through ENCODE and 264 phenotypically normal tissue samples assayed through Roadmap. The “Histone Modification” tracks within the ENCODE Integrative Analysis Data Hub to the Genome Browser public hubs page (http://genome.ucsc.edu/cgi-bin/hgHubConnect) and “Uniformly Signal” within the Roadmap Epigenomics Data Complete Collection at Wash U VizHub (http://vizhub.wustl.edu/VizHub/RoadmapReleaseAll.txt) can be used to query data from multiple cell lines/tissues that have been processed through a uniform quality pipeline. Importantly, these datasets have been taken through a uniform processing pipeline developed by the ENCODE and Roadmap Analysis Working Groups (AWG) to reduce cross-lab differences and ensure comparability between datasets (http://genome.ucsc.edu/ENCODE/qualityMetrics.html).

DHS demarks regions of open chromatin by partial digestion of intact chromatin.9 Chromatin bound by TFs tends to be more accessible than chromatin bound to histones. As such, DHS can be used to identify nonspecific protein binding at finer resolution than histone modifications. ChIP-seq methods provide an independent mapping of TF occupancy, which can reveal the identity of the TF contributing to a DHS, and also can narrow peaks of histone modification to a relatively small region (a few hundred base pairs).15 Although the binding of TFs creates a DHS peak, protein occupancy blocks nuclease cutting at the nucleotides specifically bound by protein. Therefore, at nucleotide resolution a TF-bound region would be expected to have low DHS signal flanked by two strong DHS peaks (see Fig. 3). This pattern, referred to as a DHS footprint, has been used to provide nucleotide resolution of evidence for TF binding. Algorithms have been designed to detect this signature for genome-wide detection of TF occupancy in the absence of appropriate antibodies for ChIP-seq experiments.16 However, since DHS signals are an aggregate of many cells, for which protein occupancy is inconsistent, the boundaries of these two flanking peaks are not always clearly defined and can diminish the DHS valley in protein-binding regions. Thus, although this can be a powerful alternative for detecting novel protein-binding regions, the sensitivity of DHS may be imperfect.

Figure 3

DHS footprinting annotation.

Notes: The rs6983267 variant is a laboratory confirmed functional variant in a distal MYC enhancer associated with colorectal cancer risk. Within the DHS peaks (shown in green) and broad ChIP-seq signal for TCF4 (shown in red), we find a slight valley between two higher DHS peaks. This valley corresponds to a region that actually binds TCF4 and can be detected through DHS footprinting (shown in blue). Motif analysis reveals that rs6983267 falls within the conserved binding motif of TCF4. It is important to note that because there are a number of TFs that bind to this region, as is the case with many enhancers, the valley is less pronounced than would be expected of a strong regulatory element binding a single TF.

ChIP-seq TF binding site (TFBS) data provide evidence for regions of the genome that bind a specific TF in a given cell line or tissue sample. Using ChIP-seq datasets, sequences of several hundred base pairs are identified, each nominally containing one or more motifs that bind to a specific TF. After ChIP-seq identification of motifs that bind to a particular region, position weight matrix (PWM) motifs can be used to screen for candidate regions that might bind a given TF, or to identify allelic motifs that might perturb TF binding. Motif recognition databases trained from ChIP-seq and similar methods include JASPAR,17 ConSite,18 and HaploReg19 PWM. Although the identification of altered affinities for allelic binding sequences positioned within a specific protein’s ChIP-seq signal provides strong supporting evidence for a variant having function, many ChIP-seq peaks lack a clear copy of the consensus motif. It is possible that the lack of known motifs reflects TF binding to weak motifs influenced by cooperative binding of multi-TF complexes or the existence of unknown binding motifs. It is important to note that many consensus TF motifs throughout the genome are not necessarily ever accessible for binding to the relevant TF. As such, the co-occurrence of both epigenetic remodeling and protein binding from independent assays is considered stronger evidence of a regulatory element than the presence of either signal alone. Each of these annotation features (histone marks, DHS, and ChIP-seq) captures a different snapshot of regulatory potential of a region, making the combination of annotations stronger evidence of function than any single annotation alone.

Segmentation through machine learning techniques

For nsSNPs there exist a number of annotation tools that exploit existing annotation data for their predictions, such as LS-SNP,20 SNPs&GO,21 PMUT,22 and SNAP.23 In noncoding regions, similar techniques can be used to aggregate biofeatures and predict regulatory elements. Several groups have independently applied machine learning methodologies to combine epigenetic datasets with chromatin structure, such as those available through ENCODE and the RoadMap project, into discrete annotations of regulatory regions.24,25 Segway implements a dynamic Bayesian network method using the aforementioned biofeatures (ChIP-seq and DHS signals) to identify patterns associated with transcription start sites, gene boundaries, enhancers, and other transcription regulators in an unsupervised approach.24 Software and genome browser tracks using ENCODE data are available at http://noble.gs.washington.edu/proj/segway/. Alternatively, ChromHMM first labels each biochemical assay as high or low signal in 200-bp bins across the genome and then runs a 25-state multivariate Hidden Markov Model that can similarly use the presence or absence of chromatin signals to map chromatin states.25 ChromHMM and Segway data for ENCODE cell lines can be downloaded from https://genome.ucsc.edu/ENCODE/downloads.html under the “Genome Segmentations” link within the ENCODE Analysis Hub. Similarly, segmentations for the Roadmap data are available from the University of California Santa Cruz (UCSC) Table Browser within the “Roadmap Epigenomics Data Complete Collection at Wash U VizHub” under the “Roadmap ChromHMM” track. The primary difference between these two annotations is resolution. While segments in published ChromHMM segmentation25 have a mean of 4,862 bp and a median of 800 bp,25 Segway segments have a mean length of 168 bp and a median of 124 bp.24 Although these assays could use a large number of labels, the number of states with easily meaning (such as transcription start sites, promoters, enhancers, or untranslated regions [UTR]) becomes equivocal after approximately 25 states because the clustering of histone modifications cannot readily be linked with a known chromatin state. However, an advantage of machine learning techniques is the capacity to elucidate new putatively functional epigenetic states, in addition to summarizing complex data for easier interpretation (see Fig. 2). It should be acknowledged that our understanding of the true underlying biology is incomplete, and given the paucity of gold standard data, assigning broader categories reflecting more general biological properties is usually advisable for functional follow-up.

Figure 2

Bioinformatics annotation of a promoter element using UCSC Genome Browser.

Notes: Figure 2 Bioinformatics annotation using the UCSC Genome Browser for the 8q24 index SNP rs6983267. In the first panel, the index SNP is shown to be located in an intergenic region over 200 kb upstream of MYC. The Roadmap tracks in the second panel zoom in on 8q24 showing the genomic position chr8:128,406,826–128,419,067. The chromosomal positions of three variants in LD (r2 > 0.8) with rs6983267 are shown in the first custom bed file track. The next two plots are biological replicates of normal rectal enhancer ChIP-seq signal. The following two show DHS peak enrichment in fetal intestine. Shown below the DHS plots are genomic footprinting tracks and below that is aggregate ChIP-seq signal from 161 TFs. Genomic segmentation by ChromHMM in normal colon and rectal tissue labels the region containing rs6983267 as an enhancer. In the third panel we see strong TCF7L2 and CTCF binding across the rs6983267 variant. While rs6983267 has a strong bioinformatic functional annotation, the other two variants in strong LD, rs10505477 and rs12682374, do not align to predicted regulatory elements. In vitro evaluation of this region has shown that the rs68983267 is a functional variant that disrupts the binding of TCF7 and expression of the target gene MYC.56,84

In addition to classifying segments of the genome, a new score, referred to as the combined annotation dependent depletion (CADD),26 has been developed to annotate the genome at single nucleotide resolution using a support vector machine (SVM) on features identified from previously mentioned genome-wide annotations. In contrast to Segway and ChromHMM, the C-score from CADD is one of the first attempts to assign a measure of “deleteriousness” to all possible alleles across the human genome. CADD shows great promise in the identification of Mendelian alleles (as in Thalassemia or Kabuki syndrome), with total penetrance and large effects on fitness. However, it is unclear how sensitive CADD is in detecting functional alleles with more subtle effects. By definition, deleterious alleles are functional, but functional alleles might not always cause a significant reduction in evolutionary fitness. In particular, alleles associated with late-onset conditions may not exhibit evolutionary constraint as they have modest, quantitative effects in disease risk. Thus, the negative predictive value of modest CADD (C-score < 5) scores will need to be evaluated in panels of known functional variants for common, late-onset diseases.

Chromatin interaction annotations

In general, chromatin immunoprecipitation (ChIP) is a powerful technique for studying protein–DNA interactions. In addition to detecting histone modifications and TF binding, various ChIP methods can be used to detect long-range interactions, such as the interaction of distal enhancers with target gene promoters. A detailed review of ChIP-based methods has previously been published by Fullwood and Ruan.27 In short, chromosome conformation capture (3C)28 methods are able to detect chromatin architecture, but are expensive, burdensome to implement, are unable to provide insight at a genome-wide scale, and are difficult to interpret because of a large noise-to-signal ratio. In addition, these data are often nonspecific, resulting from the stochastic properties of chromatin interactions. Indeed, two points closely positioned on the chromosome are more likely to collide at random than those at greater physical distance.29 This randomness can lead to the detection of a large amount of nonfunctional interactions relative to true functional signal in all chromatin capture techniques. Higher-throughput methods such as 4C,30,31 which can detect one regulatory element interacting with other targets, or 5C,32 detection of many regulatory elements interacting with targets, are still not implemented at a genome-wide scale.

Adding ChIP to chromosome capture methods helps to reduce noise by requiring the presence of a particular protein. However, chromosome capture methods still restrict the scope of detection because they depend on specific sites, and the addition of ChIP further limits the analysis to a single factor. Under the assumption that functional interactions between chromosomal regions are tethered, it was proposed that the use of sonication rather than the restriction enzyme digestion can help reduce noise by removing weak interactions associated with the random proximity of two regions. Chromatin interaction analysis with paired-end tag sequencing (ChIA-PET) is the first of these techniques to move toward genome-wide scale. Unlike other chromatin capture techniques, ChIA-PET introduces a linker sequence to the DNA fragments that are anchored together by a protein. Introduction of the linker allows this approach to be conducted genome-wide. The addition of both ChIP and sonication to ChIA-PET is a powerful method to detect long-range interactions across the genome associated with particular TFs.33

Expressed quantitative trait loci (eQTL) annotation

The majority of cancer GWAS do not identify plausible coding variants obviously connected to a single gene. ENCODE chromatin conformation capture experiments in modest subsets of the genome suggest that many intergenic regulatory elements skip over the closest promoter to interact with more distal genes.15,34 As a result, ad hoc identification of the causal gene(s) (eg, reporting the gene that is physically closest to the index SNP) that underlie an association signal remains a significant problem that can limit the biological interpretability of disease association study results. However, eQTLs are by definition associated with a specific gene, and trait-associated variants that are also eQTLs can generate candidate causal genes for further functional studies. The Genotype-Tissue Expression (GTEx) database is a useful resource for accessing eQTL data and can be accessed through http://www.ncbi.nlm.nih.gov/gtex/GTEX2/gtex.cgi. HaploReg19 is an additional tool for identifying variants that are positioned in eQTLs for various traits. Pairing eQTL analysis with association results to identify target genes is a powerful tool for functional follow-up because nonspecific enhancer activity assays using weak promoters are generally meaningful, but less compelling than activity assays measuring the specific interactions between an enhancer and the natural target promoter. As such, knowledge of the target gene allows for stronger experimental design before entering the lab.

Resources for Annotation of Variants

In addition to the difficulty of predicting regulatory effects of sequence variation, follow-up of cancer risk loci is complicated by the statistical correlation between nearby variants, a phenomenon referred to as linkage disequilibrium (LD). Since genetic architecture differs among populations, the resolution of GWAS results is defined by the LD structure of the studied population. Since the vast majority of cancer loci were originally discovered in populations of European descent, there are often a very large number of potential functional variants tagged by the index variant. As such, the first step in bioinformatic follow-up of functional loci is to identify all variants tagged by the index variant in the population from which the variant was discovered.

There are now many databases that catalog human variation across populations including 1000 Genomes Project,35 UK10 K project (www.uk10k.org), and a similar effort initiated in Saudi Arabia to sequence the genomes of up to 100,000 individuals (Saudi human genome program http://rc.kfshrc.edu.sa/sgp/). These catalogs can be used to identify variants that are in LD with the initial index SNP, and LD can be retrieved through various resources including Haploview36 or HaploReg19 from the Broad Institute, and SNAP.23 In addition to obtaining LD information, Haploview36 and SNAP23 allow visualization of LD across different populations through two different regional plots. Among these resources, HaploReg provides LD information for the following super populations: Europeans (EUR), Africans and African Americans (AFR), East Asians (ASN), and Hispanics (AMR) from the 2011 release of the 1000 Genomes Project. If a more specific population is desired, then Haploview is the most appropriate among these three. SNAP calculates LD in CEU, YRI, CHB/JPT populations but is restricted to 500 kb. These resources should be consulted for the identification of the bin of tagged SNPs for a given GWAS association. Early maps of common variation (the HapMap) may only have cataloged a portion of the variants within a bin, so it is advisable to consult the more recent sequencing data in order to more comprehensively identify the panel of candidate functional variants within a bin.

Annotation of variant files

As described above, there are many different datasets that can be used to annotate variants. In addition, there are many resources that have been developed to take lists of candidate functional SNPs (variant files) and integrate a large set of functional annotations for all SNPs in the file. Annotation files are typically formatted as bed track files (http://genome.ucsc.edu/FAQ/FAQformat.html#format1) for the desired functional annotations and in variant call format (VCF) (http://samtools.github.io/hts-specs/VCFv4.1.pdf) for the variants to be annotated. The functional annotation bed files can be downloaded directly from the UCSC Table Browser37 or through the MySQL site (genome-mysql.cse.ucsc.edu). There are many tools designed to annotate, visualize, or manipulate VCF files such as VCFtools,38 Annovar,39 and AnnTools.40 Other programs, such as HaploReg,19 can take lists of SNP identifiers to annotate variants with regulatory annotations in the attempt to streamline bioinformatics follow-up of GWAS loci. HaploReg takes the input variant(s), generates a list of correlated variants with the input variants, and then annotates all listed variants with pre-loaded ENCODE, Roadmap, or expression data. Recently, a bioinformatics tool called Enlight (http://enlight.usc.edu/) was developed to help identify causal variants by taking regional plots from GWAS results and overlaying biofeatures such as epigenetic modification, DHS, and TFBS. One limitation of both HaploReg and Enlight is that the meta-collection available from these sites reflects only a limited proportion of the available functional annotations, which may not be relevant to the disease of interest, and may or may not be the most updated version of the dataset. Furthermore, to make these resources usable, the LD calculation may not be in the most relevant population, and the calculation may not encompass a large enough region of the genetic architecture at that particular locus.

A relatively new annotation tool, GWAS3D (http://jjwanglab.org/gwas3d), is a web server that combines chromatin interaction data (5C, Hi-C, ChIA-PET) with histone modification, TF binding (ChIP-seq), and ChromHMM segmentation in 16 cell types from ENCODE to predict the probability that a genetic variant affects regulatory pathways. To predict functional variants, GWAS3D then looks to differences in binding affinity for the predicted regulatory factors. This information is then used to annotate GWAS signals, taking into consideration both LD structure of super populations and effect estimates, to aid in the interpretation of results. Input data can be VCF or PLINK (http://pngu.mgh.harvard.edu/purcell/plink/) files, and the implementation is straightforward. A potential limitation of this resource is that it is not currently made tissue specific. As such, it is possible that this resource could miss regulatory elements that are tissue specific and related only to certain cancers. However, this restriction is reflective of the scarcity of chromosomal interaction data for many tissues and cell types.

Cancer-specific annotations

The annotations described thus far have focused on the annotation of noncoding germline variation to predict causal variants within susceptibility loci. These annotations can be generalized to functional annotation of associations for any disease or trait. However, there are also a few cancer-specific resources. COSMIC41 and TCGA (https://tcga-data.nci.nih.gov/tcga/) focus on acquired somatic variants in cancer and how to detect those that play a significant role in tumor development, proliferation, or treatment response.42,43 Genes that are frequently amplified, deleted, dysregulated, or somatically mutated in the cancer of interest can be prioritized as potential regulatory targets for nearby index SNP associations. Similar to eQTL and chromosome capture analyses, knowledge of the target gene for a candidate regulatory enhancer can be used to design stronger functional analysis experiments, as shown for the MYC region and rs6983267 in colon cancer.44,45 In order to identify novel driver mutations, several large-scale sequencing projects of tumor tissue have commenced, such as TCGA. Many of the available annotations for identifying driver mutations relate somatic mutations to expression data, or leverage existing knowledge about functional mutations, cancer genes, or associated pathways. A comprehensive list of available programs and algorithms for distinguishing driver mutations was provided in a recent publication by Zhang et al.46 Many of the tools focus on annotation of the variants themselves using measurements of deleteriousness (eg, PolyPhen2). However, these tools also include methods that compare the genes harboring aberrations to known cancer genes in order to predict novel drivers. More sophisticated algorithms consider mutational profiles or background rates for passenger mutations or use machine learning methods to train a classifier on a known set of driver and passenger mutations, which can be collected from resources like the COSMIC database,41 and then apply those predictors to a new sequencing dataset (eg, Driver Mutation Identification [DMI],47 Cancer-Specific High-throughput Annotation of Somatic Mutations [CHASM],48 and Screening for Non-Acceptable Polymorphism [SNAP]).49 Many of the annotations that specifically try to identify driver mutations are based on the hypothesis that a driver gene would accumulate a greater number of mutations than passenger genes. Major limitations of these tools include differences in predictions based on the method used, and that passenger mutations can assume the role of driver mutations given a change in the environment, such as targeted chemotherapy.46 In addition, it is now appreciated that there is a great deal of tumor heterogeneity, both intratumor (within a single tumor) and intertumor (between tumors) that make it difficult to generalize parameters trained on a limited number of examples.

Visualization of variant annotations

Databases used to visualize annotated variants include the UCSC Genome Browser,50 the UCSC Cancer Genome Browser,51 and as described earlier, Enlight. The UCSC Cancer Genome Browser currently displays a growing catalog of data, including 201 datasets from 22 TCGA (The Cancer Genome Atlas) cancers as well as data from Cancer Cell Line Encyclopedia and Stand Up To Cancer. This database allows samples to be examined by common clinical features such as therapeutic response or by genomic signatures that predict drug response. Genomic data can be viewed by genes, which allows users to easily see functional changes to the genome as well as examine trends across pathways of genes. Several statistical tools are available to test quantitative differences in expression as well. Additionally, the Tumor Image Viewer, based on Google Maps, allows users to interactively view slides of tumor tissue samples. The UCSC Genome Browser contains the ENCODE, Roadmap, and many other functional data-sets, which can be used to identify variants that could impact splicing, transcription, translational, and post-translational regulation of gene expression. After identifying a list of strong functional candidates it is often useful to upload a bed file containing potentially functional variants underlying GWAS signal to visualize their proximity to genes and their positioning relative to segments of the genome likely to have regulatory potential.

Proposed Bioinformatics Framework

Given the volume of annotations and resources presently available, it is useful to establish a bioinformatics framework to guide functional evaluation of noncoding variants associated with cancer risk. As such, a framework is presented below as a general approach for conducting such analyses.

To follow up an index SNP that has been associated with cancer, it is first necessary to generate a list of variants tagged by the index SNP in the studied population. There are now a number of public resources available to achieve this (such as HaploReg). After obtaining a list of correlated variants, it is useful to create a custom bed file to visualize these variants in the UCSC Genome Browser. In Figure 2 the bin of variants associated with an index SNP (rs6983267) for colorectal cancer risk (P = 5 × 10−11)52 is smaller than most containing only three SNPs. This bin of variants is positioned 300 kb upstream of the protooncogene MYC. It is not unusual for variants to be located in intergenic regions that are tens of thousands of base pairs away from the closest gene (see Table 2) making it difficult to predict the potential target gene. The cancer genome browser (https://genome-cancer.ucsc.edu/)51 and the Online Mendelian Inheritance in Man (OMIM)53 are useful resources for identifying biologically plausible genes. For instance, an OMIM search of MYC reveals that its gene product plays a key role in cell proliferation, differentiation, transformation, and apoptosis. Although this locus is now well established as an enhancer for MYC, in many instances the target is less clear. For example, two different breast cancer loci have been followed up in laboratory, and functional variants, rs691357854 and rs4784227,55 were found to impact reporter protein expression but have unknown targets. Further chromosome capture in relevant tissue could provide substantial insight to help translate these loci into better clinical targets.

Table 2

Examples of variants with laboratory follow-up.

VARIANT	REGULATORY REGION	GENE	TRAIT	PUBLICATION	YEAR	BEFORE GWAS?
rs6983267	Distal Enhancer 300kb downstream	MYC	Colorectal Cancer	Pomerantz et al., Nat Genet56	2009	Follow-up
rs16888589	Enhancer 20kb downstream	EIF3H	Colorectal Cancer	Pittman et al., PLoS Gen85	2010	Follow-up
rs10822013	Intronic	ZNF365	Breast Cancer	Cai et al., Hum Mol Genet86	2011	Follow-up
rs2735940	Promoter	TERT	Telomere Length	Matsubara et al, Bichem Biophys Res Commun	2006	Pre-GWAS
rs2735940	Promoter	TERT	Acute Lymphoblastic Leukemia (ALL)	Sheng et al., Carcinogenesis69	2013	Follow-up
rs2736108	Promoter	TERT	Ovarian and Breast Cancer	Beesley et al., PLoS One87	2011	Follow-up
rs1512268	Enhancer 14kb downstream	NKX3.1	Prostate Cancer	Akamatsu et al., Hum Mol Genet88	2010	Follow-up
rs6913578	Regulates unknown target	Unknown	Breast Cancer	Cai et al., Cancer Res54	2011	Follow-up
rs4784227	Possible enhancer of TOX3, 18.4 kb upstream	Unknown	Breast Cancer	Long et al., PLoS Genet55	2010	Follow-up
rs2239632	Promoter	CEBPE	ALL	Ryoo et al., J Hum Genet89	2013	Follow-up
rs11730582	Promoter	OPN	Gastric Cancer	Zhao et al., BMC Cancer90	2012	Follow-up
rs11730582	Promoter	OPN	Melanoma	Schultz et al., Mol Carcinog91	2009	Follow-up
rs4590952	Intronic	KITLG	Testicular Cancer	Zeron-Medina et al., Cell92	2013	Follow-up
rs944289	Enhancer	PTCSC3 (lincRNA)	Papillary thyroid carcinoma (PTC)	Jendrzejewski et al., Proc Natl Acad Sci93	2012	Follow-up
rs1859961	Distal Enhancer 1Mb upstream	SOX9	Prostate Cancer	Zhang et al., Genome Res94	2012	Follow-up
rs12194974	Promoter	LIN28B	Ovarian Cancer	Permuth-Wey et al., Cancer Res95	2011	Follow-up
rs8506	miRNA binding in exon	lincRNA-NR_024015	Gastric Cancer	Fan et al., PLoS One96	2014	Follow-up
rs10993994	Promoter	MSMB	NA	Buckland et al., Hum Mutat70	2005	Pre-GWAS
rs10993994	Promoter	MSMB	Prostate Cancer	Chang et al., Hum Mol Genet71	2009	Follow-up
rs10993994	Promoter	MSMB	Prostate Cancer	Lou et al., Proc Natl Acad Sci72	2009	Follow-up

Bioinformatics follow-up of loci positioned in intronic or intergenic regions can be conducted within the UCSC Genome Browser, or another annotation tool. After aligning variants of interest to the reference genome, it is useful to annotate the region with biofeatures (ChIP-seq and DHS) assayed in relevant tissues and cell lines from Roadmap and ENCODE, respectively. Figure 2 illustrates that unlike the other two tagged variants (rs10505477 and rs12682374), the index SNP (rs6983267) aligns to a region of open chromatin with promoter histone modifications in normal colon and rectal mucosa. Furthermore, rs6983267 is located within a TFBS for TCF4 making this a likely functional candidate for laboratory follow-up. Figure 3 illustrates that although the ChIP-seq signal was quite broad, DHS footprinting identifies a much smaller region that demarks the true underlying binding site for TCF4 and that rs6983267 falls in the center of that footprint. ChromHMM segmentation in normal colon and rectal mucosa illustrates the ease of interpretation of the annotation, labeling the region harboring rs6983267 as an enhancer in rectal mucosa. However, the variant misses this region in colon tissue. This variant has now been confirmed as causal in the laboratory, with the risk allele decreasing the binding affinity of TCF4.56 For variants within UTRs one can look for miRNA binding sites using the TargetScanS 5.2 track within the UCSC Genome Browser. Variants falling near splice sites can be annotated using a number of different tools including BDGP, ESE finder,57 RescueESE,58 ESRSearch,59 PESX,60 or Ensembl.61

Limitations of Existing Resources

In noncoding regions, the majority of functional predictions to date inform whether a polymorphism is positioned in a region that is likely to be biologically relevant based on the presence of one or more biochemical signals. As such, with the exception of CADD and PWM motif analysis, noncoding annotations do not predict the effects of a specific allele. Furthermore, even when a variant has a particularly strong bioinformatics annotation relative to other tagged SNPs, laboratory follow-up of multiple variants is often required. For instance, Figure 4 outlines an example where the variant with the strongest bioinformatics annotation (rs16892766) is not the same variant shown to have allelic effects on gene expression in vitro. The index tagSNP rs16892766, which is strongly associated with colorectal cancer (1 × 10−18),62–67 is located in an intergenic region over 140 kb telomeric of the promoter for the closest gene EIF3H. The bioinformatics annotation of this locus shown in Figure 4 suggests that rs16892766 is positioned in a poised enhancer in colon and rectal tissues and falls within a DHS footprint that corresponds to ChIP-seq binding sites for ESR1, FOXA1, and NR3C1. Alternatively, rs16888589 is in strong LD with rs16892766 (r2 = 0.89 in 1000 Genomes Project EUR population), has similar enhancer signals, but does not fall within a binding site for an assayed TF. Fine-mapping of both variants through in vitro and in vivo analysis suggested that in human colorectal cancer cell lines, the ancestral allele rs16888589-A significantly repressed gene expression while there were no allelic effects observed with rs16892766. The authors also reported that rs16888589 interacts with the promoter of EIF3H, as shown through chromatin capture techniques. Interestingly, increased expression of this gene was shown to increase colorectal cancer growth and invasiveness. Thus, despite the relatively stronger bioinformatic annotation for the index SNP, another predicted functional SNP, rs16888589, was shown to be the likely causal variant. This illustrates the importance of laboratory follow-up for multiple variants, particularly when there is more than one putative functional variant predicted.

Figure 4

Limitation of bioinformatics annotation.

Notes: Figure 4 outlines an example where the variant with the strongest annotation (shown in red) is not the underlying functional variant (shown in green). The rs16892766 locus is located in an intergenic region over 20 kb downstream of EIF3H. Although both rs16888589 and rs16892766 fall in a poised enhancer based on histone modification, DHS, ChromHmm segmentation, and footprinting evidence, only the index rs16892766 falls in the highly conserved ChIP-seq binding site for the three shown TFs. Laboratory follow-up of this variant revealed that the true underlying causal variant is rs16888589 and not the index, although the index had the stronger bioinformatic annotation illustrating the importance of laboratory follow-up, particularly when there exists more than one predicted functional variant.

Although higher throughput functional assays might increase the production of training data to improve predictions on the functionality of noncoding variants, it is possible that a great deal of pre-GWAS knowledge of functional variation is lost in the literature. Table 2 outlines several examples of variants that were confirmed to be functional through in vitro analysis after the initial association was detected in GWAS. However, this table also reveals that several known functional variants were subsequently investigated as a follow-up of a GWAS finding without reference to the original experiment. For instance, the TERT promoter variant, rs2735940, was shown to have functional effects on gene expression in 2006,68 and then again after a recent GWAS finding for acute lymphoblastic leukemia in 2013.69 Similarly, the MSMB promoter variant, rs10993994, which is associated with prostate cancer, has recently been investigated in the laboratory three times, while the original experiment showing the allelic effects of rs10993994 was conducted in 2005.70–72 A significant challenge for extracting functional evidence from the literature is that prior to GWAS, promoters (rarely enhancers) may have been investigated because the gene was interesting for a trait that may be unrelated to the subsequent disease or trait association detected through GWAS. These examples indicate that there is a need to more effectively link the laboratory investigations of allelic variants to GWAS findings in order to avoid redundant efforts.

Future Directions

Advances in sequencing technology have enabled genetic association studies of complex diseases to test rare variant hypotheses in both coding and noncoding regions. Unlike traditional genome-wide approaches, which make use of markers to capture bins of correlated SNPs, sequencing-based approaches identify a much larger number of rare variants, significantly reducing the power of a single variant approach. To address this loss in power, aggregate association tests such as burden and variance-component analyses have been developed to effectively collapse the genome into smaller units of analysis, and decrease the multiple testing penalty. Each method collapses a large number of infrequent observations into a much smaller set of testing units. In the exome, a gene serves as a natural functional unit for combining rare variants. However, segmenting intergenic space into functional units is more difficult because enhancers work in concert to regulate one or more target genes. Although chromatin capture technology can provide substantial insight into how distal regions interact, these assays are only available for a limited number of tissues and cell types because of the associated cost, and there remains a large number of nonfunctional predicted interactions. Since rare variants are intermittently scattered across intergenic space, aggregation within a single enhancer unit, as defined through segmentations like ChromHMM or Segway, is not feasible since each element is likely to contain only a single variant. Furthermore, the resolution of ChIP-seq is on the order of thousands of nucleotides and is therefore far too broad to distinguish between variants closely positioned on the chromosome. Although the additional requirement of DHS overlap can refine this signal to a few hundred nucleotides, the underlying binding motif for functional machinery (TFs) is about 10 nucleotides (5–31 nucleotides, with mean 9.9 nucleotides for eukaryotes).73 Furthermore, these annotations do not provide insight into whether the minor allele of the polymorphism would impact the binding of the TF. As such, in aggregate testing methods, it is difficult to prioritize likely functional variants over those less likely to be functional in the same testing unit. Since this is an important component for improving power in aggregate testing, this is a significant limitation that will need to be addressed.

Moving forward it will be important to develop more robust and high-throughput methods to identify both functional and nonfunctional allelic changes in motif sequences and their relative impact on the expression of a target. The generation of these data will lead to improved training data-sets, as is currently available for nsSNPs. However, it will also be important to create a better reporting system to link this information into a publicly available resource to avoid redundant efforts and maximize the utility of these more onerous experiments. Although there is now a wealth of information on the functional potential of noncoding regions, there remains much to be discovered on how these regions interact, and the impact of regulatory variation.

Given the time and expense associated with laboratory follow-up, in silico prediction remains a critical step for generating testable hypotheses. The recognition of these functional, regulatory polymorphisms is becoming feasible as datasets annotating evidence for functional noncoding regions across the genome accrue. Hopefully, this will facilitate the identification of the actual variant(s) in each region associated with disease. In the short term, genetic risk models based on functional variants rather than index tagSNPs should yield more precise genetic risk models, as well as facilitate the translation of GWAS associations into therapeutic targets. In the longer term, as methods for the prediction of regulatory function improve, these methods might be used to identify less frequent regulatory polymorphisms associated with cancer, even when such variants are too rare for statistical correlations in reasonably sized populations, and thereby account for some of the missing heritability of these diseases.