FUMA covers various features of existing tools

As a variety of bioinformatics tools have been developed to obtain insights in GWAS results^23–25, we compared the list of features available in FUMA with the features available in other tools, and describe these further below (Table 1).

LD calculation is the first step to characterize risk loci of GWAS by computing population specific LD structure, so called clumping which identifies independent significant SNPs and defines the genomic risk loci. PLINK²⁶ is the most widely used software for this task which takes GWAS summary statistics (requiring a reference panel) or genotype data as input. In FUMA, this task is automated by using pairwise LD (r ²) of SNPs in the reference panel (1000 genomes project phase 3²⁷) pre-computed by PLINK, resulting in a list of independent significant SNPs, lead SNPs and genomic risk loci based on the GWAS input file. FUMA also adds SNPs to the identified risk loci that do not have a P-value (i.e., they were not available in the GWAS input file), but that are LD proxies of the identified lead SNPs, as these SNPs might be causally relevant. Alternatively, users can pre-compute lead SNPs or risk loci and upload these to FUMA.

Variant Annotation is required to obtain information on biological consequences of SNPs in the risk loci. There are several tools such as ANNOVAR¹² and VEP²⁸ which annotate functional consequences on genes, and variant scores such as deleteriousness and phylogenetic conservations (extensive review is available in Hou and Zhang²⁹). Particularly for non-coding SNPs, SCAN³⁰, RegulomeDB¹⁴ and HaploReg³¹ annotate regulatory information, such as eQTLs, enhancer/promoter regions, and transcription factor binding sites (see Tak and Farnham³² for extensive overview). Although SCAN and HaploReg correct for LD, the input of the tools mentioned above is a list of SNPs of interest which does not take genetic associations into account and thus requires pre-processing of GWAS results by the user. FUMA performs annotation of SNPs that are in LD of independent significant SNPs in a single flow, and does not require additional data preformatting.

Gene-based test/gene-set analyses are methods that enable to summarize SNP associations at the gene level and associate the set of genes to biological pathways. For instance, VEGAS performs permutation based simulation^33,34, MAGMA employs multiple linear regression³⁵ and Pascal computes sum and maximum of chi-squared statistics³⁶ to obtain gene-based P-values. Additionally, there are several tools that perform not only gene-based test but also gene-set analyses using full distribution of genetic associations (e.g., MAGMA³⁵, MAGENTA³⁷, INRICH³⁸, and DEPICT³⁹). FUMA implements MAGMA gene-based analysis and gene-set analysis on the full GWAS input data. In addition, genes prioritized by SNP2GENE or by the user are also tested for overrepresentation in various gene sets in GENE2FUNC process.

Visualization is one of the essential features that allows (quick) insights into the GWAS results, e.g., summarizing annotated information of SNPs and genes. LocusZoom is one of the most widely used visualization tool for GWAS results which plots LD structure of a risk locus, gene locations as well as SNP association values⁴⁰. LocusTrack is an extension of LocusZoom which also plots additional information together such as Chip-seq and chromatin state⁴¹. 3D Genome Browser is a recently developed web application which contains comprehensive 3D chromatin interaction datasets such as Hi-C and ChIA-PET⁴², though it does not integrate with GWAS summary statistics. These tools are primarily focused on visualization of a subset of functionally relevant data sources. FUMA integrates results from multiple lines of evidence and provides interactive visualization of results, facilitating rapid interpretation.

The current lack of a single platform that integrates all possible resources for post-GWAS annotation hampers our understanding of GWAS results, as different GWAS studies may use a different selection of queried resources rendering their post-GWAS interpretation incomplete and difficult to compare. FUMA provides a central place for a wide variety of post-GWAS annotation strategies and to our knowledge is the most versatile tool in doing so.

Application to GWAS of body mass index

To validate the utility of FUMA, we applied it to summary statistics of the most recent GWAS for body mass index (BMI; 236,231 individuals)⁴³. FUMA identified 95 lead SNPs (from 223 independent significant SNPs) across 77 genomic risk loci (Fig. 2 and Supplementary Data ¹–³), in accordance with the original study. We first conducted positional mapping of deleterious coding SNPs and eQTL mapping (Methods) which prioritized 151 unique genes; 23 genes with deleterious coding SNPs (positional mapping), and 144 genes with eQTLs that potentially alter expression of these genes (eQTL mapping) including 16 genes that had both deleterious coding SNPs and eQTLs (Supplementary Data ⁴). The 151 genes consist of 55 genes that were also reported in the original study⁴³ and 96 novel genes implicated by FUMA, including 45 genes which are located outside the risk loci. These novel candidates have shared biological functions with the 55 previously known candidate genes such as “metabolism of carbohydrate”, “metabolism of lipid and lipoprotein”, “immune system”, and “calcium signaling” (Supplementary Data ⁵). In addition, FUMA results showed that, although several genomic loci for BMI included multiple prioritized genes, a single gene was prioritized in 22 out of 43 loci which contain at least one prioritized gene (Supplementary Fig. ²), suggesting that these 22 genes have a high probability of being the causal gene in that region. The 22 “highly likely causal genes” include several well-known genes for BMI such as NEGR1, TOMM40, and TMEM18. The strongest GWAS association signal for BMI was on 16q.12.2 where three genes were prioritized; FTO, RBL2, and IRX3 (Fig. 3). These three genes were only prioritized by eQTL mapping as the positional mapping showed no deleterious coding SNPs located in these genes. The original study⁴³ only mentioned FTO, because the associated SNPs were located in this gene, however none of the associated SNPs have a potential direct affect such as coding SNPs on FTO. Two of the genes prioritized by FUMA (RBL2 and IRX3) are physically located outside the genomic locus and are missed when using conventional approaches that prioritize genes located in the locus of interest based on LD around the top SNP. Although the IRX3 gene was not reported in the original study⁴³, recent functional work has indeed validated this as the causal gene whose expression is affected by SNPs in the 16q.12.2 locus⁴⁴.

Open in a separate window

Fig. 2

Overview of prioritized genes from BMI GWAS by FUMA. Starting from the BMI GWAS summary statistics, boxes represent results of the SNP2GENE process. The annotated SNPs include all independent lead SNPs and SNPs which are in LD with these lead SNPs. Prioritized genes are divided into three categories; genes that are implicated by deleterious coding SNPs (colored pink), by eQTLs for these genes (colored blue), or by chromatin interactions (colored green). The prioritized genes are further categorized into previously reported genes (blue) and novel genes (red) prioritized genes by FUMA. ^*These genes were not prioritized by FUMA since they do not have either deleterious coding SNPs, eQTLs or chromatin interactions, although they are located within GWAS risk loci

Open in a separate window

Fig. 3

Regional plot of the locus 16q.12.2 of BMI GWAS. a Extended region of the FTO locus, which includes prioritized genes RBL2 and IRX3. Genes prioritized by FUMA are highlighted in red. b Zoomed in regional plot of FTO locus with, from the top, GWAS P-value (SNPs are colored based on r ²), CADD score, RequlomeDB score and eQTL P-value. Non-GWAS-tagged SNPs are shown at the top of the plot as rectangles since they do not have a P-value from the GWAS, but they are in LD with the lead SNP. eQTLs are plotted per gene and colored based on tissue types. In the plots of CADD score, RegulomeDB score and eQTLs, SNPs which are not mapped to any gene are colored gray

We then performed chromatin interaction mapping using Hi-C data of 14 tissue types (Methods). FUMA prioritized 310 genes (Supplementary Data ⁴), of which 61 genes are overlapping with the genes prioritized by positional and/or eQTL mappings and 232 genes are located outside of the genomic risk loci (Fig. 2). That resulted in a total of 400 prioritized genes by combining three mapping strategies including 330 novel candidates which were not reported in the original study (Table 2 and Supplementary Data ⁴). These novel candidates further supported shared biological functions with previously reported known genes, such as lipid and lipoprotein metabolism, homeostatic process and various metabolic pathways, with a greater number of genes compared to the mappings without Hi-C data (Supplementary Data ⁵). Out of 400 prioritized genes, 59 genes are mapped by both eQTLs and chromatin interactions including IRX3 on the 16q.12.2 locus (Fig. 4), which further supports the hypothesis that these genes are involved in the risk of BMI. Of the 48 loci that contained at least one prioritized gene from positional and eQTL mappings, chromatin interaction mapping identified candidate genes in additional 18 loci (Supplementary Fig. ²), including loci mapped to known genes associated with BMI such as MC4R, FOXO3, and ADCY9. The 400 prioritized genes showed enrichment in 9 GO terms, such as “response to zinc ion” and “oligopeptide binding” overlapping with multiple metallothionein and glutathione S-Transferase genes whose association with obesity risk has been reported^45,46 (Supplementary Data ⁶).

Open in a separate window

Fig. 4

Chromatin interactions and eQTLs of BMI risk loci on chr. 16. The most outer layer is the Manhattan plot displaying SNPs with P-value<0.05. Candidate SNPs are colored based on the highest r ² to one of the independent significant loci (red: r ²>0.8, orange: r ²>0.6). Other SNPs are colored in gray. rsID of top SNPs per locus are labeled. The outer circle is the chromosome coordinate and genomic risk loci are highlighted in blue. Genes mapped by either Hi-C or eQTLs are shown on the inner circle. Genes mapped by Hi-C, eQTLs are colored orange and green, respectively. Genes mapped by both are colored red. Chromatin interaction and eQTLs are shown as links colored orange and green respectively

Thus, using BMI summary statistics, FUMA confirmed known genes but also prioritized novel genes, including potential causal genes located outside the GWAS risk loci of BMI, which were missed in the original study.

Application to Crohn’s disease GWAS

To further illustrate its utility, we applied FUMA to the summary statistics of Crohn’s disease⁴⁷ (CD; 6333 cases and 15,056 controls). With FUMA, 95 lead SNPs from 184 independent significant SNPs across 71 genomic loci were identified for CD (Supplementary Fig. ³ and Supplementary Data ⁷–¹⁰). First, describing the results of positional mapping of deleterious coding SNPs and eQTL mapping, FUMA prioritized 95 unique genes from 32 loci (Supplementary Fig. ⁴), of which 39 genes were implicated by deleterious coding SNPs and 69 were implicated by eQTLs influencing expression of these genes (12 genes had both deleterious coding SNPs and eQTLs; Table 2 and Supplementary Data ¹¹). The prioritized 95 genes include 37 known candidate genes that were also reported in the original study⁴⁷ including well-known CD-related genes such as NOD2, IL23R, and SLC22A5, while 58 genes were novel (Supplementary Fig. ³; see Supplementary Note ⁴ and Supplementary Figs. ⁵–⁷ for detail results). These novel candidates include 18 genes that are physically located outside the GWAS risk loci, and the novel candidates mainly share immune system related biological functions with 37 previously known genes (Supplementary Data ¹²).

Chromatin interaction mapping using Hi-C data in small bowel and liver prioritized 199 genes of which 18 genes are overlapping with genes prioritized by positional and/or eQTL mappings and 149 genes are located outside of the genomic risk loci (Supplementary Data ¹¹). That resulted in a total of 276 prioritized genes including 215 novel candidates which were not reported in the original study (Table 2 and Supplementary Fig. ³). Of the 23 loci which are mapped to at least one gene by positional and eQTL mappings, additional 23 loci are mapped to candidate genes by chromatin interaction mapping, in which several of the genes prioritized from those loci are involved in immune system and cytokine signaling pathways (Supplementary Fig. ⁴ and Supplemental Data ¹²). One of these 23 risk loci, the 17q12 locus is mapped to six chemokine ligands by Hi-C in liver: CCL1, CCL2, CCL7, CCL8, CCL11, and CCL13. Additionally, prioritized genes include 11 cytokines (IL4, IL5, IL10, IL19, IL23R, IL24, IL27, IL33, IL1RL1, IL18R1, and IL18RAP) wherein IL18R1 and IL18RAP are also mapped by eQTLs in whole blood and IL23R and IL27 are also mapped by deleterious coding SNPs which further supports the involvement of these cytokine genes in CD. The role of these chemokines and cytokines in inflammatory disease has been widely studied⁴⁸ and yet, chromatin interaction mapping identified additional relevant candidates from the risk loci. The prioritized 276 genes showed enrichment in 123 canonical pathways such as immune system and cytokine related pathways which are known to be highly relevant to CD⁴⁹ (Supplementary Data ¹³).

Characterization of genomic risk loci based on GWAS

To define genomic loci of interest to the trait based on provided GWAS summary statistics, pre-calculated LD structure based on 1000G of the relevant reference population (EUR for BMI, CD and SCZ) is used. First of all, independent significant SNPs with a genome-wide significant P-value (<5e-8) and independent from each other at r ²<0.6 are identified. For each independent significant SNP, all known (i.e., regardless of being available in the GWAS input) SNPs that have r ²≥0.6 with one of the independent significant SNPs are included for further annotation (candidate SNPs). These SNPs may thus include SNPs that were not available in the GWAS input, but are available in the 1000G reference panel and are in LD with an independent significant SNP. Candidate SNPs can be filtered based on a user-defined minor allele frequency (MAF, ≥0.01 by default).

Based on the identified independent significant SNPs, independent lead SNPs are defined if they are independent from each other at r ²<0.1. Additionally, if LD blocks of independent significant SNPs are closely located to each other (<250kb based on the most right and left SNPs from each LD block), they are merged into one genomic locus. Each genomic locus can thus contain multiple independent significant SNPs and lead SNPs.

Besides using FUMA to determine lead SNPs based on GWAS summary statistics, users can provide a list of pre-defined lead SNPs. In addition, users can provide a list of pre-defined genomic regions to limit all annotations carried out by FUMA to those regions.

Annotation of candidate SNPs in genomic risk loci

Functional consequences of SNPs on genes are obtained by performing ANNOVAR¹² (“gene-based annotation”) using Ensembl genes (build 85). Note that SNPs can be annotated to more than one gene in case of intergenic SNPs which are annotated to the two closest up- and down-stream genes. CADD scores, RegulomeDB scores and 15-core chromatin state are annotated to all SNPs in 1000G phase 3 by matching chromosome, position, reference, and alternative alleles. eQTLs are also extracted by matching chromosome, position and alleles of all independent significant SNPs and SNPs which are in LD with one of the independent significant SNPs for each user-selected tissue type, wherein SNPs can have multiple eQTLs for distinct genes and tissue types (Supplementary Note ²). Information on previously known SNP-trait associations reported in the GWAS catalog is also retrieved for all SNPs of interest by matching chromosome and position.

Gene mapping

Gene annotation is based on Ensembl genes (build 85). To match external gene IDs, ENSG ID is mapped to entrez ID yielding 35,808 genes which consist of 19,436 protein-coding genes, 9249 non-coding RNA, and other 7123 genes (e.g., pseudogenes, processed transcripts, immunoglobulin genes, and T-cell receptor genes).

Positional mapping is performed based on annotations obtained from ANNOVAR¹². Two optional filters are provided to control the maximum distance from SNPs to genes and select specific functional consequences of SNPs on gene. When the former option is defined, FUMA maps SNPs to genes based on ANNOVAR annotation and a user-defined maximum distance is applied for intergenic SNPs. When the latter option is provided, FUMA maps only SNPs which have selected annotations annotated by ANNOVAR (e.g., coding or splicing SNPs).

For eQTL mapping, all independent significant SNPs and SNPs in LD of them are mapped to eQTLs in user-defined tissue types. By default, only significant SNP–gene pairs (false discovery rate (FDR)≤0.05) are used. Optionally, eQTLs can be filtered based on a user-defined P-value. eQTL mapping maps SNPs to genes up to 1Mb apart (cis-eQTLs).

Chromatin interaction mapping is performed by overlapping independent significant SNPs and SNPs in LD of them with one end of significantly interacting regions in user-selected tissue/cell types. These SNPs are then mapped to genes whose promoter regions (250bp up- and 500bp down-stream of transcription start site by default) overlap with another end of the significant interactions. Optionally SNPs can be filtered for those overlapping with predicted enhancer regions of the user-selected epigenomes. Similarly, mapped genes can also be filtered for having promoter regions overlap with predicted promoter regions of the user-selected epigenomes.

Optional filtering of SNPs based on functional annotations obtained in step 2 of SNP2GENE (i.e., CADD score, RegulomeDB score, 15-core chromatin state) can be performed for positional, eQTL and chromatin interaction mappings separately. When any of these filters is activated, candidate SNPs are filtered primary to gene mapping. Note that this filtering of SNPs based on functional annotations for a certain mapping does not affect other mappings, e.g., when SNPs are filtered by CADD score in positional mapping but not in eQTL mapping, SNPs are filtered prior to positional mapping but eQTL mapping uses the original set of candidate SNPs.

For mapped genes, two scores of intolerance to functional mutations are annotated; probability of being loss-of-function intolerant (pLI)⁵⁸ and non-coding residual variation intolerance score (ncRVIS)⁵⁹.

MAGMA for gene analysis and gene set analysis

FUMA uses input GWAS summary statistics to compute gene-based P-values (gene analysis) and gene set P-value (gene set analysis) using the MAGMA³⁵ tool. For gene analysis, the gene-based P-value is computed for protein-coding genes by mapping SNPs to genes if SNPs are located within the genes. For gene set analysis, the gene set P-value is computed using the gene-based P-value for 4728 curated gene sets (including canonical pathways) and 6166 GO terms obtained from MsigDB v5.2. For both analyses, the default MAGMA setting (SNP-wise model for gene analysis and competitive model for gene set analysis) are used, and the Bonferroni correction (gene) or FDR (gene-set) was used to correct for multiple testing. 1000G phase 3²⁷ is used as a reference panel to calculate LD across SNPs and genes.

GTEx gene expression data set

Normalized gene expressions (reads per kilo base per million, RPKM) of 53 tissue types were obtained from GTEx (Supplementary Table ³). A total of 56,320 genes was available in GTEx, which we filtered on an average RPKM per tissue greater than or equal to 1 in at least one tissue type. This resulted in transcripts of 28,520 genes, of which 22,146 were mapped to entrez ID (see “Gene Mapping” section for details). In the GENE2FUNC, the heatmap of prioritized genes displays two expression values; (i) the average log2(RPKM+1) per tissue per gene, in which RPKM is winsorized at 50, allowing comparison of expression level across genes and tissue types and (ii) the average of the normalized expression (zero mean of log2(RPKM+1)) per tissue per gene allowing comparison of expression level across tissue types within a gene.

To obtain differentially expressed gene sets (DEG; genes which are significantly more or less expressed in a given tissue compared to others) for each of 53 tissue type, the normalized expression (zero mean of log2(RPKM+1)) is used. Two-sided Student’s t-tests are performed per gene per tissue against all other tissues. After the Bonferroni correction, genes with corrected P-value <0.05 and absolute log fold change ≥0.58 are defined as a DEG set in a given tissue, i.e., for these gene expression in the given tissue had the largest discrepancy with expression in all other tissues. In addition, we distinguish between genes that are upregulated and downregulated in a specific tissue compared to other tissues, by taking the sign of t-score into account. In GENE2FUNC, genes are tested against those DEG sets by hypergeometric tests to evaluate if the prioritized genes (or a list of genes of interest) are overrepresented in DEG sets in specific tissue types.

Gene set enrichment test

To test for overrepresentation of biological functions, the prioritized genes (or a list of genes of interest) are tested against gene sets obtained from MsigDB (i.e., hallmark gene sets, positional gene sets, curated gene sets, motif gene sets, computational gene sets, GO gene sets, oncogenic signatures, and immunologic signatures) and WikiPathways, using hypergeometric tests. The set of background genes (i.e., the genes against which the set of prioritized genes are tested against) is 19,283 protein-coding genes. Background genes can also be selected from gene types as described in the “Gene Mapping” section. Custom sets of background genes can also be provided by the users. Multiple testing correction (i.e., Benjamini–Hochberg by default) is performed per data source of tested gene sets (e.g., canonical pathways, GO biological processes, hallmark genes). FUMA reports gene sets with adjusted P-value ≤0.05 and the number of genes that overlap with the gene set >1 by default.

FUMA parameters for application to GWAS summary statistics

In the described applications, three mapping strategies were applied to GWAS summary statistics with the following settings: positional mapping was performed by selecting exonic and splicing SNPs with CADD score ≥12.37 (defined by Kircher et al.¹³) to restrict the mapping to deleterious coding SNPs. eQTL mapping was performed using GTEx eQTLs with FDR<0.05. Chromatin interaction mapping was performed using Hi-C data from Schmitt et al.¹¹ and interactions were filtered by FDR<1e-6. Tissue types used for eQTLs and chromatin interaction mappings are described in the following section for each of three phenotypes. Other parameters not mentioned here were kept as default (Supplementary Table ²).

Application to BMI GWAS

Parameters were set as described in the above section and we used eQTLs in 44 tissue types from GTEx. For chromatin interaction mapping, Hi-C data of 14 tissue types (Adrenal, Aorta, Bladder, Dorsolateral Prefrontal Cortex, Hippocampus, Left Ventricle, Liver, Lung, Ovary, Pancreas, Psoas, Right Ventricle, Small Bowel and Spleen) from GSE87112 was used. Indels were excluded. rsID was mapped to dbSNP build 146 and chromosome and positions were extracted based on human genome hg19 reference. Only protein-coding genes were used in gene mapping and enrichment of DEG in 53 tissue types, Canonical Pathways and GO terms were tested.

Application to CD GWAS

We set parameters as described above and we used eQTLs in five tissue types from GTEx which are relevant to CD, i.e., Small Intestine, Colon Sigmoid, Colon Transverse, Stomach, and Whole Blood. Chromatin interaction mapping was performed using Hi-C data of two tissue types; Liver and Small Bowel from GSE87112. The MHC region and indels were excluded from the analysis. Since the input GWAS summary statistics only contained results from the discovery phase, we manually submitted the 71 reported lead SNPs to FUMA in addition to the independent lead SNPs that were identified as described above (Supplementary Data ⁷). Only protein-coding genes were used in mappings and enrichment of DEG in 53 tissue types, Canonical Pathways and GO terms were tested.

Application to SCZ GWAS

Parameters were set as described above and eQTLs in 10 brain tissues from GTEx. Chromatin interaction mapping was performed using Hi-C data of two brain regions; hippocampus and prefrontal cortex. The extended MHC region (25–34Mb), Chromosome X and indels were excluded from this analysis. The input GWAS summary statistics are based on the discovery phase and not all reported lead SNPs from the combined results of discovery and replication phases reached genome-wide significance. To include all reported lead SNPs, 111 non-indel lead SNPs were provided to FUMA and additional independent lead SNPs were identified at P<5e-8 (Supplementary Data ¹⁴). Only protein-coding genes were used in mappings and enrichment of DEG in 53 tissue types, Canonical Pathways and GO terms were tested.

Code availability

Source code of FUMA web application is available through a git repository at https://github.com/Kyoko-wtnb/FUMA-webapp/.

This work was funded by The Netherlands Organization for Scientific Research (NWO VICI 453-14-005) and Ingrosyl. We thank the GIANT consortium, WTCCC and PGC for providing GWAS summary statistics and GTEx Portal for RNA-seq and eQTL data. We also thank Prof Patrick Sullivan for discussion of 3D chromatin interaction data.