HaploReg v4: systematic mining of putative causal variants, cell types, regulators and target genes for human complex traits and disease.

More than 90% of common variants associated with complex traits do not affect proteins directly, but instead the circuits that control gene expression. This has increased the urgency of understanding the regulatory genome as a key component for translating genetic results into mechanistic insights and ultimately therapeutics. To address this challenge, we developed HaploReg (http://compbio.mit.edu/HaploReg) to aid the functional dissection of genome-wide association study (GWAS) results, the prediction of putative causal variants in haplotype blocks, the prediction of likely cell types of action, and the prediction of candidate target genes by systematic mining of comparative, epigenomic and regulatory annotations. Since first launching the website in 2011, we have greatly expanded HaploReg, increasing the number of chromatin state maps to 127 reference epigenomes from ENCODE 2012 and Roadmap Epigenomics, incorporating regulator binding data, expanding regulatory motif disruption annotations, and integrating expression quantitative trait locus (eQTL) variants and their tissue-specific target genes from GTEx, Geuvadis, and other recent studies. We present these updates as HaploReg v4, and illustrate a use case of HaploReg for attention deficit hyperactivity disorder (ADHD)-associated SNPs with putative brain regulatory mechanisms.

INTRODUCTION

Phenotype-associated loci from genome-wide association studies (GWAS) are usually non-coding, and functionally interpreting them is a challenge due to linkage disequilibrium (LD) and our almost complete inability to predict regulatory function directly from non-coding sequence. Therefore, regulatory genomic data such as maps of enhancers and transcription factor binding sites are essential to interpreting GWAS, developing mechanistic hypotheses, and ultimately understanding the genetic architecture of complex traits and disease (1–3). For human geneticists, these regulatory data can be unwieldy to translate from a genome browser to insights about a set of genomically-dispersed disease variants. HaploReg (4) integrates regulatory genomic maps together in the context of haplotype blocks, allowing researchers to intersect regulatory elements with genetic variants to quickly formulate functional hypotheses, both through dissection of multiple variants within a haplotype block and through global enrichment analysis of a set of associated loci. HaploReg annotation of GWAS has successfully been applied for haplotype fine-mapping (5–9) and enrichment analysis (7,10,11).

DATA AND INTERFACE UPDATES

HaploReg has been expanded substantially since it first launched in 2011. Here we describe the updates that have been incorporated in Haploreg v4 in response to new research in regulatory genomics and feedback from users.

Catalog of variants

HaploReg v4 defines a core set of 52 054 804 variants, consisting primarily of single-nucleotide polymorphisms (SNPs) using all refSNP IDs, hg19 positions and alleles from dbSNP release b137 (12). Corresponding hg38 coordinates for these variants were obtained from dbSNP release b141. This core set of dbSNP variants was integrated with other data sets either by rsID (for GWAS, eQTL and 1000 Genomes data) or by intersecting intervals by coordinate using the BEDTools software package (13) (for all other functional tracks.)

Linkage disequilibrium was calculated using phased low-coverage whole-genome autosomal sequences for four ancestral super-populations (AFR, AMR, ASN and EUR) from the 1000 Genomes Project Phase 1 release (14), using a search space of all variants within 250 kilobases of each other. Allele frequencies were also obtained for each population.

Location of variants relative to genes was calculated using BEDTools and both GENCODE (15) and RefSeq (16).

Genome-wide association studies

GWAS were obtained from the EBI-NHGRI GWAS Catalog (17) (downloaded 30 October 2015). When there were multiple GWAS for the same trait, a trait-wide pruning was performed to retain only the strongest (lowest P-value) GWAS result from all studies on that trait, when two results from different studies were overlapping or within one megabase of each other.

Sequence conservation

Mammalian evolutionarily constrained elements are defined as originally reported, using both SiPhy elements (18) and GERP elements (19). Both of these comparative genomics studies report both base-level conservation scores as well as discretized elements; we chose to report discretized elements resulting from the authors’ algorithms for the sake of simplicity and interpretability. A colored cell represents that that the element is conserved according to the algorithm.

Regulatory protein binding

Protein-binding sites from a variety of cell types and experimental conditions was obtained from the ENCODE Project ChIP-Seq data (20), processed by the narrowPeak algorithm.

Reference epigenomes

Epigenomic data from the Roadmap Epigenomics project (11) for the following data sets were included: ChromHMM states corresponding to enhancer or promoter elements, from the 15-state core model and 25-state model incorporating imputed data (21); histone modification ChIP-seq peaks using the gappedPeak algorithm for H3K27ac, H3K9ac, H3K9me1 and H3K9me3; and DNase hypersensitivity data peaks using the narrowPeak algorithm.

Expression quantitative trait loci

Expression QTL (eQTL) results were obtained from the GTEx pilot analysis v6 (22), the GEUVADIS project (23) and 12 other studies (10,24–34) in order to annotate variants with their putative regulatory target genes and the tissue(s) in which genotype has been associated with gene expression level. A wide range of QTLs, including eQTLs and other molecular QTLs such as metabolite QTLs, were also extracted from the GRASP database, build 2.0.0.0 (35,36).

Regulatory motifs

A library of position weight matrices from commercial, literature and motif-finding analysis of the ENCODE project (37) was used to score the effect of variants on regulatory motifs using the position weight matrix (PWM)-scanning process described previously (4).

Enrichment analysis

For a given set of lead SNPs from a GWAS or user-input SNPs, the overlap of SNPs with predicted enhancers in each reference epigenome is assessed. Users have four different options for defining enhancers, available in the option panel: using the 15-state core model, using the 25-state model incorporating imputed epigenomes, using H3K4me1 peaks and using H3K27ac peaks. The overlap with enhancers in each cell type is compared to two background models to assess enrichment: all 1000 Genomes variants with a frequency above 5% in any population and all independent GWAS catalog SNPs. The enrichment relative to these background frequencies is performed using a binomial test and uncorrected P-values are reported in an enrichment table underneath the haplotype views.

USE EXAMPLE

To become acquainted with HaploReg, use the GWAS drop-down menu to select ‘Attention deficit hyperactivity disorder (Lesch KP, 2008, 26 SNPs)’ and select ‘Submit’. Notice that the first two haplotype blocks from this study (38) are driven by lead SNPs with the same P-value = 1 × 10−8. Go to the second haplotype result, for lead SNP rs864643 (Figure 1). Note that the top row in the haplotype block shows the SNP rs561543, and that it has LD of r2 = 0.81 and D′ = 0.95 with the lead variant rs864643. It overlaps with an HMM-predicted enhancer in four major tissue types; hover over ‘4 tissues’ in that row to see a variety of enhancer tissues, including brain. Note that there is also an experiment with HNF4 protein bound by ChIP-Seq, 9 QTL results and an HNF4 motif disruption.

Figure 1.

Example of haplotype summary view.

Notice the enrichment results at the bottom of the page below the haplotype results. Note that the strongest enrichment for enhancers (as defined by the 15-state core ChromHMM model) is in the angular gyrus sample from brain, with binomial P = 2.0 × 10−6 relative to all common SNPs.

Then go to the entry on the block for the lead SNP itself, rs864643. Click on the rsid, which is colored red because it is the lead SNP. Note that in the full table of epigenomic information from Roadmap Epigenomics (11), there is a cluster of enhancer activity in brain, and that it is classified as a genic enhancer by the 15-state core model and transcribed 3′ enhancer by the 25-state model (Figure 2). Note that H3K4me1, H3K27ac and H3K9ac all contribute to the chromatin state assignment at this locus. Black cells on the right hand of this part of the table indicate that DNase was not assayed by Roadmap in these tissues.

Figure 2.

Example of epigenome details in a SNP detail view.

Go to the bottom of the detail page for rs864643. Note that the SNP has been correlated with MOBP expression in two brain tissues (29), MPRL15 expression in blood (39) and serum ratio of allantoin to quinate (40); all three of these studies were curated by GRASP and found by cross-referencing this SNP to its database (35,36) (Figure 3). Looking at studies individually curated by HaploReg, notice that the SNP has been associated with differential expression of a single exon of RPSA in lymphoblastoid cells by the GEUVADIS study (23). In the motif table, note that the SNP changes the match to the p300 PWM, ATTAYRWCA, with the alternate allele changing a match to the fourth A to a G. Hover over the ‘p300_disc’ ID to see that the motif was discovered using the Trawler algorithm on a p300 ChIP-Seq experiment in HeLa cells from the ENCODE dataset (37).

Figure 3.

Example of eQTL and motif alteration details from a SNP detail view.

These lines of evidence suggest regulatory mechanisms by which the SNPs from this GWAS may affect the complex phenotype of ADHD. While individually each piece of evidence is relatively weak, they offer ways in which molecular biologists could proceed with further experiments that would more definitively establish mechanisms. For example, the GWAS-wide enrichment suggests global differential gene regulation in angular gyrus, which has been associated with hyperactivation in ADHD by fMRI (41) and suggests a tissue to study gene expression directly in animal models. ChIP-seq and motif data suggest specifically testing HNF4 binding differentially to the alleles of rs561543, and the strong motif coupled with eQTL data suggest looking at whether p300 binds differentially to rs864643 in a brain tissue model. Finally, MOBP eQTL evidence suggests experiments to dissect the mechanism of MOBP differential expression, perhaps modulated by p300 at rs864643 and suggests that it may be useful to perform ADHD-relevant behavioral assays of MOBP-deficient mice, which do not show an overt behavioral phenotype (42).