ENCODE consortium hotspots ²⁶ were obtained from ftp://ftp.ebi.ac.uk/pub/databases/ensembl/encode/integration_data_jan2011/byDataType/openchrom/jan2011/combined_hotspots/. Roadmap Epigenome DNase1 sequencing tag alignments were obtained from http://www.genboree.org/EdaccData/Current-Release/experiment-sample/Chromatin_Accessibility/. The files used correspond to that part of the Gene Expression Omnibus (GEO) accession GSE18927 beyond the data use embargo date. These alignments were processed by the Hotspot ( http://www.uwencode.org/proj/hotspot/) ^{28, 29} method with the default parameters to give hotspot and peak files. For this analysis we choose to use the hotspots which are regions of general DNase I sensitivity rather than peaks which are more similar to DNase I hypersensitive sites because, although the method works with peaks, hotspots reveal more tissue specific signal (data not shown). Cell and tissue assignments for each of the datasets were made using the decodings available from the ENCODE Data Coordination Center tables ( https://genome.ucsc.edu/encode/cellTypes.html) or from the BioSamples database ( http://www.ebi.ac.uk/biosamples/) sampleGroup SAMEG31306. A list of samples used is provided in Supplementary Table S1.

The complete collection of SNPs discovered in GWAS studies curated in the NHGRI GWAS Catalog ⁵ were downloaded from http://www.genome.gov/gwastudies/ ³⁰ (Accessed 3rd September 2014). SNPs were grouped according to the annotation provided in the Disease/Trait field and only sets with 5 or more non-redundant SNPs were retained. See Supplementary Table S2 for list of SNP sets analysed. A set of files of the SNPs included in analysis for each phenotype is available in the data directory of the GitHub release, https://github.com/iandunham/Forge. For FORGE analysis SNP sets were further filtered by LD pruning removing all but one SNP from a set of SNPs at r ² >= 0.8 in the 1000 genomes data (see below) and were analysed for both ENCODE and Roadmap Epigenome DNase I hotspots, selecting background SNP sets from the default GWAS genotyping array SNPs.

The FORGE tool utilises either an SQLite (command line tool, http://www.sqlite.org) or MySQL (web tool, http://www.mysql.com) database of the overlaps of every 1000 genomes project ( http://www.1000genomes.org) ³² SNP with the ENCODE and Epigenome Roadmap DNase1 hotspots. To prepare this database, SNPs from the 1000 genomes phase 1 integrated call data set were downloaded from ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/analysis_results/integrated_call_sets and compared to indexed DNase I hotspots using tabix from the SAMtools package ( http://samtools.sourceforge.net/tabix.shtml) ³³ using a distributed approach on the EBI compute farm. The overlaps for each SNP were stored in a single large indexed table of SNP location and identifier with binary strings representing the presence (1) or absence (0) of overlap in each sample for each of the hotspot datasets (ENCODE or Roadmap).

To prepare sets of background SNPs matched to the test SNP set, FORGE matches SNPs based on G+C content (GC), minor allele frequency (maf) and transcription start site (TSS) distance, and repeats the overlap analysis for each of 1000 background sets. Overall population mafs were obtained from the 1000 genomes project phase 1 integrated call dataset at ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/analysis_results/integrated_call_sets. To control for the processes involved in selecting SNPs for genotyping, only 1000 genomes phase 1 SNPs that had been included on one of the common genotyping platforms as described in Ensembl ( http://www.ensembl.org/info/genome/variation/data_description.html#variation_sets) were considered further. This left either 1875813 SNPs across various platforms (Affy GeneChip 100K Array, Affy GeneChip 500K Array Affy SNP6, HumanCNV370-Quadv3, HumanHap300v2, HumanHap550v3.0, Illumina Cardio Metabo, Illumina Human1M-duoV3, Illumina Human660W-quad) or 2231212 SNPs across the Illumina HumanOmni2.5 array. TSS distance was determined for each remaining SNP relative to the TSS defined by the Gencode project ^{34, 35} given in ftp.ebi.ac.uk/pub/databases/ensembl/encode/integration_data_jan2011/byDataType/gencode/jan2011/Gencodev7_CAGE_TSS_clusters_June2011.gff.gz using Bedops closest-features ³⁶ . GC was determined for a 100 bp window centred on the SNP at base 50. The SNPs were then sorted into 1000 bins partitioned by deciles for each of GC, maf and TSS. For each SNP in a test set, the corresponding bin is identified based on its GC, maf and TSS distance, and background selections are made from that bin.

A set of SNPs can be presented to FORGE as a list of RSIDs or by genome location on human genome build GRCh37 in either VCF or Bed formats. If RSIDs are not provided in one of these formats, the genome coordinates are used to identify the RSID. SNPs not present in the 1000 genomes phase 1 integrated call dataset are excluded from the analysis. With LD pruning selected a single SNP (the first in the file) is chosen from LD clusters within either r ² >= 0.8 or r ² >= 0.1 as specified. For each analyzable SNP in the test set, overlaps are retrieved from the FORGE database, and a count of total hotspot overlaps is recorded for each DNase I sample (cell) for the test SNP set. One thousand matching background SNP sets containing the same number of SNPs as the test SNP set are selected, matched for GC, maf and TSS distance by decile bins as described above. Overlaps for each of the SNPs in each of the background SNP sets are also retrieved from the database and an overlap count for each background set in each DNase I sample is recorded. For each test SNP set, the background probability of overlap is determined from the 1000 background set overlap counts and the probability of the observed test result under a binomial distribution is calculated. The P value thresholds of 0.05 and 0.01 are corrected for multiple testing by division by the number of tissue groupings tested, and the corrected threshold is used. The use of tissue as the unit for sample grouping is consistent with the groupings obtained by hierarchical clustering of samples by DNase I data (results not shown). The corrected thresholds are therefore more stringent than established by the random trials.

FORGE generates both tabular and graphic descriptions of the enrichment of overlap for the test SNPs for each DNase I hotspot sample. A tab-separated values (TSV) file is output including columns for the binomial P value, cell, tissue, filename of the sample hotspots, SNPs that contribute to the enrichment, and the GEO accession for each sample. These data are also provided as an interactive table produced using the Datatables ( https://datatables.net/) plug-in for the jQuery Javascript library accessed through the rCharts package ( http://ramnathv.github.io/rCharts/).

Each of the graphic outputs presents the –log ₁₀ binomial p by cell sample. A pdf graphic is generated using base R graphics ( http://www.r-project.org‎). The interactive Javascript graphic is generated using the rCharts package ( http://ramnathv.github.io/rCharts/) to interface with the dimple d3 libraries ( http://dimplejs.org). In both cases cells are grouped alphabetically by tissue, and for the pdf alphabetically by cell. The interactive graphic stacks replicate samples at the same x coordinate. In each of the graphics the colouring of results by P value is consistent, blue (p > 0.05 equivalent after correction), pink (0.05 => p < 0.01), and red (p <= 0.01). The corrected P value threshold is given on the pdf output.

To estimate false positive rates, 1000 sets of SNPs at each of a series of SNP counts between 5 and 300 SNPs were randomly chosen from the 1000 genome phase 1 integrated SNP set. FORGE analysis was run for each set across the ENCODE and Roadmap Epigenome data, and the number of tests with P values less than thresholds ranging from 0.05 to 0.001 were recorded. These represent the false positives from 1000 trials at each of 424 samples i.e. 424,000 tests, and were used to calculate false positive rates at each significance threshold.

The hierarchical clustering solution was obtained using a multi-scale bootstrap resampling approach. We first computed a binary regulatory signature for each cell type classifying each DNase I site as active or inactive in each cell type sample. Hierarchically clustering of the binary regulatory signatures was by Euclidean distance with Ward's agglomerative method using the pvclust R/CRAN package with default values ( http://cran.r-project.org/web/packages/pvclust/index.html). Finally, we identified clusters supported by the data at a bootstrap probability p value < 0.01.

FORGE analysis provides a method to view the tissue-specific regulatory component of a set of variants. In its current implementation, FORGE analysis takes a set of single nucleotide polymorphisms (SNPs), such as those SNPs reported above the genome wide significance threshold (p < 5e-8) in a GWAS study, optionally filters the SNPs to remove all bar one SNP from a region in high LD (“LD pruning”) and determines whether there is enrichment for overlap with putative regulatory elements compared to a matched background of SNP sets. Initially the elements considered are DNase I hotspots generated from either the ENCODE ²⁶ or Roadmap Epigenomics projects DNase I data by the Hotspot method ^{28, 29} , because of both the comprehensiveness of the sites identified and the broad range of cell types for which DNase I data were available. DNase I hotspots can be regarded as regions of general DNase I sensitivity.

For each set of test SNPs, an overlap analysis is performed against the DNase I hotspots for each available cell sample separately (125 samples for ENCODE, 299 for Roadmap, described in Supplementary Table S1), and the number of overlaps is counted. Major potential confounders in this analysis are the many biases of GWAS SNP distribution on the genome. To account for this a background distribution of the expected overlap counts for this SNP set is obtained by identifying 1000 matched background SNP sets of the same number of SNPs, matching each test SNP with an equivalent SNP by decile bin for each of G+C content (GC), minor allele frequency (maf) and distance to the nearest transcription start site (TSS). The matched background SNPs sets are overlapped with the DNase I hotspots and the background distribution of overlap counts is determined. The enrichment of the test SNP set for each sample is expressed as the binomial P value of the test SNP set given the background overlap distribution. The FORGE results are presented in interactive and static graphical and tabular forms by cell type. Enrichments above the background distribution with binomial P values less than 0.01 corrected for multiple testing are considered significant and are highlighted in red in the graphical output. Enrichments with p <= 0.05 are also highlighted in pink ( Figure 1). As the DNase I patterns are not independent between cell types, we conducted simulation experiments with randomly selected input SNPs. We chose 1000 random test SNP sets for each of a series of SNP counts ranging between 5 and 100 SNPs and conducted FORGE analysis on both ENCODE and Roadmap data. The false positive rate was determined as the number of cell type enrichments identified greater than the significance thresholds used by FORGE expressed as the proportion of the total number of sample overlap tests performed (424,000) for each SNP count. This analysis showed that the P value thresholds are reasonably well calibrated to false positive levels around 0.5–0.75 %. In practice we correct the thresholds further for multiple testing across different tissues so as to be more stringent and so expect typical false positive levels to be less. As discussed below, many of the GWAS SNP sets do not reveal any enrichment, consistent with low false positive rates.

We have implemented FORGE as a web tool available at http://browser.1000genomes.org/Homo_sapiens/UserData/Forge. The interface accepts a list of SNPs by dbSNP RefSNP identifier (RSID) or by genomic location on human genome build GRCh37 in either Variant Call Format (VCF, http://www.1000genomes.org/wiki/Analysis/Variant Call Format/vcf-variant-call-format-version-41) or Bed format (Personal Genome SNP format, http://genome.ucsc.edu/FAQ/FAQformat.html#format10), and allows specification of the background selection from two common sets of GWAS SNP typing microarrays. LD filtering is achieved at either r ² >= 0.8 or r ² >= 0.1 using 1000 genomes project population data. The outputs of the analysis are an interactive graphic for exploration of the analysis, a static pdf for printing or publication ( Figure 1), and a table of enrichments in either an interactive or standard tab separated format.

In addition the code is available to download from https://github.com/iandunham/Forge, with the required database files available at ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/browser/forge. Installed on a Macbook Pro with core i7 processor, 16Gb RAM and a solid state hard drive, a typical FORGE analysis of a set of 55 SNPs with 1000 background tests is accomplished in around 30s, and in 35s with LD filtering.

We ran FORGE analysis on 260 phenotypes in the NCBI GWAS catalog NCBI GWAS catalog ^{5, 30} with a reported associated SNP count of 5 or more after LD pruning (see Supplementary Table S2 for list of phenotypes analysed and references) at genome-wide significance. Complete tables of results for all SNP sets analysed are included in the data directory of the Github release, https://github.com/iandunham/Forge. Thirty-five and 60 out of 260 SNP sets had at least one significant enrichment at the P value thresholds of 0.01 and 0.05 after correction for multiple testing, respectively ( Table 1). A set of example positive outputs from this analysis is available from http://www.1000genomes.org/forge-gwas-catalog-example-gallery11 (PDF format). Removing SNPs that directly alter a protein coding exon from the GWAS catalog sets did not substantively alter the patterns of enrichments (data not shown).

Figure 1 shows a series of example FORGE analyses for autoimmune disease studies on the Roadmap Epigenome samples (references for the studies are provided in Supplementary Table S2). In each case there is a clear signal for enrichment of overlap with DNase I hotspots in the blood-derived samples including cells of immune function. In more detail, for those phenotypes where there is involvement of T cell activation or invasion in the aetiology (e.g. Crohn’s disease, Multiple sclerosis) there is enrichment in the CD3, CD4 and CD8 positive samples containing T cells as well as enrichment in the CD56 positive sample including NK cells. In addition, these disease SNPs overlap with hotspots present in the fetal thymus samples, consistent with the location of T cell maturation. Further signals specific to the individual aetiologies are also identified. Crohn’s disease SNPs show enrichment of overlap with hotspots in the fetal small and large intestine samples, as well as fibroblasts and skin cells. For inflammatory bowel disease SNPs there is a much more general enrichment, in addition to the specific immune cells, which may be consistent with the more generalized inflammation. In contrast, for autoimmune diseases where the primary involvement is a B cell response (Rheumatoid arthritis (RA), Systemic lupus erythematosus (SLE)), the most pronounced overlap enrichment is in CD19 positive samples characteristic of B cell activation or circulating plasma cells. In rheumatoid arthritis there is also some overlap enrichment for samples characteristic of T cells and thymocytes, but it is relatively less, and this is much less pronounced for SLE. Thus, FORGE analysis reveals tissue-specific enrichment of overlap for GWAS SNPs with regulatory regions indicative of known tissue involvement in the disease aetiology.

The tissue-specific enrichment of overlap is not specific to just autoimmune diseases (see results gallery at http://www.1000genomes.org/forge-gwas-catalog-example-gallery11). For instance, for QRS duration the GWAS associated SNPs are strongly enriched for overlap with fetal heart samples. GWAS SNPs associated with pulmonary function measured by spirometry are enriched for overlap with hotspots in fetal lung cells and lung cell lines. For red blood cell traits and platelet count the major overlap enrichment signal is in CD34 positive hematopoietic progenitor cells consistent with their role in both red blood cell and platelet development. In contrast for GWAS SNPs involved in height, the overlap enrichment is not tissue-specific but is more general over many tissues and cell lines. There are further examples displayed in the results gallery at http://www.1000genomes.org/forge-gwas-catalog-example-gallery11, in most cases consistent with expected disease aetiologies.

FORGE is available through a web interface at http://browser.1000genomes.org/Homo_sapiens/UserData/Forge.

The source code for FORGE is available on GitHub at https://github.com/iandunham/Forge with the Forge.db sqlite database and background selection hash tables available at ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/browser/forge_11.

FORGE has been successfully been installed and run on Mac OSX 10.8.4 and Red Hat Linux.

https://github.com/iandunham/Forge

https://github.com/F1000Research/Forge/releases/tag/v1.1

http://dx.doi.org/10.5281/zenodo.13900 ³⁸

GNU GPL