VSEAMS: a pipeline for variant set enrichment analysis using summary GWAS data identifies IKZF3, BATF and ESRRA as key transcription factors in type 1 diabetes

Abstract

Motivation: Genome-wide association studies (GWAS) have identified many loci implicated in disease susceptibility. Integration of GWAS summary statistics (P-values) and functional genomic datasets should help to elucidate mechanisms.

Results: We extended a non-parametric SNP set enrichment method to test for enrichment of GWAS signals in functionally defined loci to a situation where only GWAS P-values are available. The approach is implemented in VSEAMS, a freely available software pipeline. We use VSEAMS to identify enrichment of type 1 diabetes (T1D) GWAS associations near genes that are targets for the transcription factors IKZF3, BATF and ESRRA. IKZF3 lies in a known T1D susceptibility region, while BATF and ESRRA overlap other immune disease susceptibility regions, validating our approach and suggesting novel avenues of research for T1D.

Availability and implementation: VSEAMS is available for download (http://github.com/ollyburren/vseams).

Fig. 1.

The VSEAMS pipeline; mandatory inputs are shaded grey; a dashed border indicates that one or the other input is required. VSEAMS takes as input either two lists of genes or two lists of regions for comparison. Given genes, regions are defined by taking gene coordinates ± 200 kb around the TSS. GWAS summary statistics (P-values) for SNPs in those regions are extracted. The observed Wilcoxon rank sum test statistic is compared with its null distribution determined by its theoretical mean and a variance derived by simulating null P-values with a correlation structure matching the underlying genotype structure. Caching of pregenerated LD matrices reduces computation time. A full description of each step is available in the Supplementary Information

Fig. 2.

A comparison of Z-scores generated using permuted phenotype method (10 000 permutations) versus using summary P-values and VSEAMS (10 000 simulations) for T1DGC study, over 100 randomly generated gene sets

Fig. 3.

A runtime comparison of simulation using multivariate normal (black) versus permutation (grey) over 1000 randomly selected LD blocks. In both plots the y-axis is the median execution time over 10 iterations, and lines indicate the fitting of a linear model. Specifically, (a) details the effect of sample size on median execution time over 14 753 SNPs summed over all randomly selected LD blocks. (b) Shows the effect of SNP count on execution time for 4000 cases and controls for all 1000 randomly selected LD blocks

Fig. 4.

T1D susceptible SNP enrichment (excluding major histocompatibility complex (MHC)) within transcription factor perturbed gene sets from Cusanovich et al. (2014) SNPs are pruned on the basis of r² threshold $\ge$0.95. A positive Z-score indicates enrichment, labels denote associated P-values. Black bars indicate that the knocked-down transcription factor overlaps a known autoimmune susceptibility locus curated in ImmunoBase

Fig. 5.

Comparison of VSEAMS and permuted phenotype methods with differing sample size, for example, gene sets, where enrichment is present (IKZF3) and absent (YY1). (a) Shows difference in Z-scores between both methods with 10 000 simulations and a variable sample size, with an equal number of cases and controls. (b) Shows how the correlation between Z-scores over a variable number of permutations varies with respect to sample size. The coloured lines represent a locally estimated scatterplot smoothing (LOESS) fitted model for each sample size