WebGWAS: A web server for instant GWAS on arbitrary phenotypes

2.1 Overview of WebGWAS

We built WebGWAS.org, a free, public, efficient web app that can rapidly compute approximate GWAS summary statistics. WebGWAS allows users to define arbitrary phenotypes, compute GWAS results, visualize them, and download them. No protected health information (PHI) is stored on the WebGWAS server—only GWAS summary statistics, phenotypic covariance matrices, and anonymized phenotypic data with k = 10 anonymity. No user data are collected beyond which phenotype definitions are submitted. At present, WebGWAS contains data for ICD-10 codes and quantitative measurements in the White British subset of the UK Biobank. Figure 1 shows an overview of WebGWAS.

• Download figure
• Open in new tab

Figure 1: Graphical overview of WebGWAS.

A: Three artifacts are needed for WebGWAS: GWAS summary statistics, a covariate-adjusted phenotypic covariance matrix, and anonymized phenotypes. No protected health information (PHI) is included in these outputs, which can be uploaded to the WebGWAS server. B: Computations on WebGWAS involve five main steps. First, a user defines a phenotype using the web interface. Second, that phenotype definition is evaluated in anonymized data to produce a vector of values. Third, that vector is approximated by the anonymized features using linear regression. Fourth, the fit quality of that linear regression is computed. Fifth, the GWAS summary statistics for the user’s phenotype definition are computed using indirect GWAS. This process only uses the phenotype approximation coefficients, the covariance matrix of the features, and the pre-computed GWAS summary statistics for the features.

2.2 Indirect GWAS method

WebGWAS depends on a statistical approach that we term indirect GWAS. The core idea of indirect GWAS is that GWAS summary statistics for a linear combination of phenotypes can be computed without individual-level data, provided certain summary statistics are available for all phenotypes included. For example, a GWAS on the sum of systolic and diastolic blood pressure can be computed given only summary statistics about systolic and diastolic blood pressure, separately. What follows is a high-level overview of this method. A full derivation is provided in the methods.

Consider a genetic association test for a single quantitative phenotype (y) and a single genetic variant (g), y ~ g. Assuming both have zero mean, the coefficient and standard error are the following:

Suppose y can be written as a linear combination of m feature phenotypes x₁ … x_m, (i.e. y = ∑x_ip_i, for some coefficients p₁ … p_m). Define two m-vectors, b and s, that hold the coefficient and standard error estimates for regressions of the genetic variant against the feature phenotypes (i.e. b_i = g^Tx_i/g^Tg and s_i = SE(b_i)). Let C be the phenotypic covariance matrix of the features. Then the coefficient and standard error for y can be written in terms of feature summary statistics as

Var(g) can be estimated using the coefficients, standard errors, and variances from each feature phenotype and taking the mean across features.

To summarize, the association test summary statistics for phenotype y = ∑ x_ip_i can be computed using only the coefficients and standard errors for all x_i, the covariance matrix among x_i, and the coefficients p_i. In other words, given appropriate summary statistics about a set of phenotypes, the GWAS summary statistics for any linear combination of those phenotypes can be computed without individual-level data. In the methods we show how this approach is applicable to any linear GWAS method, and how covariates can be included.

2.3 Validation of indirect GWAS

To evaluate whether this derivation is correct, we compared to the standard method of direct GWAS on the evaluated linear combination. We used PhenotypeSimulator [7] to simulate genetic and phenotypic data for 10,000 individuals, generated 100 random linear combinations of traits, and compared GWAS summary statistics between direct and indirect approaches to validate our derivation and implementation. The results of this evaluation show that indirect GWAS is mathematically correct in the case of linear models and a good approximation for linear mixed models (Figure 2).

• Download figure
• Open in new tab

Figure 2: Validation of indirect GWAS using five GWAS methods.

GWAS chi-squared statistics are compared between direct and indirect approaches for 100 randomly projected phenotypes, 10,000 genetic variants using simulated data. Each point represents one genetic variant and one phenotype (1,000,000 total points in each facet). R² values indicate goodness of fit between the direct and indirect statistics.

2.4 Generating summary statistics for arbitrary phenotypes

Research phenotypes are not exclusively defined as linear combinations; they often include boolean operations like “and” and “or”, and inequalities like “greater than”. For example, one could define diabetics as anyone with either an ICD-10 E11 diagnosis or an HbA1c measurement ≥ 6.5%. While arbitrary phenotypes may be nonlinear, indirect GWAS requires them to be linear, so we first linearize each phenotype using linear regression.

Linear approximations perform well for many nonlinear phenotypes, including Phecodes [8] and Boolean operators (Figure 3A). Using ICD-10 code data from the UK Biobank, we generated 140 Phecodes and 300 random nonlinear Boolean combinations, and we approximated these using linear regression against ICD-10 codes. The quality of linear approximations varied greatly between phenotypes. Pairs of codes combined like “x and not y” or “x or y” were the best approximated, with median R² values of 0.96 and 0.98, respectively. Phecode performance ranged widely, with regression R² values ranging from 0.15 to 0.98, but overall they performed well with a median of 0.82. We found consistently poor performance among phenotypes defined as “x and y”, with R² values between 0.01 and 0.39. In short, many nonlinear phenotypes could be approximated with reasonable performance using a linear approximation (74% (324/440) above 0.5, 66% (289/440) above 0.75).

• Download figure
• Open in new tab

Figure 3: Indirect approach maintains GWAS fidelity for many nonlinear phenotypes.

A: Phenotype fit quality predicts GWAS fidelity. Nonlinear phenotypes (140 Phecodes and 300 boolean definitions) were approximated using linear combinations of ICD-10 codes via OLS. GWAS fidelity (Pearson R between direct and indirect chi-squared statistics) is worsened by linearization alone, and it is highly correlated (Pearson correlation 0.99) with phenotype fit quality (R² between exact and linearly approximated phenotypes). Outliers above the line of unity are inflated due to exceptionally high association statistics (see supplementary materials). B: Anonymization degrades performance but phenotype fit quality is still highly correlated with GWAS fidelity (Pearson correlation 0.76). ICD-10 code data were first anonymized using MDAV with k = 10. Nonlinear phenotypes were evaluated on the anonymized data, then the best linear combinations were estimated in anonymized rather than original data. GWAS fidelity is worsened both due to linearization and anonymization. C: Increasing privacy degrades GWAS fidelity. ICD-10 data were anonymized using MDAV at various degrees of anonymity. In each anonymized dataset, all nonlinear phenotype definitions were evaluated, linearized, and used to perform indirect GWAS. GWAS fidelity in this subplot refers to the Pearson R between anonymized and non-anonymized indirect GWAS chi-squared statistics. This definition intentionally differs from the previous subplots to isolate the effect of anonymization.

GWAS fidelity is the key metric for evaluating WebGWAS. We quantified GWAS fidelity using the Pearson correlation between direct and indirect GWAS chi-squared statistics. While GWAS fidelity is the metric of greatest interest, it cannot be computed in WebGWAS, because it requires direct GWAS summary statistics. However, we found that phenotypic fit quality (R² as above) is an excellent predictor of GWAS fidelity (Pearson correlation 0.99; Figure 3A).

This analysis led us to draw two conclusions. First, many real phenotypes (e.g. Phecodes) can be well-approximated using linear regression (e.g. against ICD-10 codes). Second, phenotypic fit quality—which can be computed in WebGWAS—is an acceptable gauge for GWAS fidelity. In summary, indirect GWAS is applicable to arbitrary phenotypes, and phenotype fit quality, which can be computed cheaply, provides a quality gauge for the results.

2.5 Protecting privacy by anonymizing phenotype data

Phenotypic data are used for one purpose in WebGWAS: to find the linear combination of features that best predicts the user-defined phenotype. We hypothesized that anonymized data would give similar coefficients to non-anonymized data, while reducing the risk of exposing sensitive information. To anonymize phenotype data, we chose maximum distance to the average vector (MDAV), an anonymization method that ensures every record is at least k-anonymous in the dataset.

Anonymization involves a loss of fidelity. We evaluated whether this loss is empirically tolerable for WebGWAS. We anonymized ICD-10 code data from the UK Biobank to k = 10, applied the previous phenotype definitions (i.e. Phecodes and Booleans) to the anonymized data, regressed anonymized nonlinear phenotypes against anonymized ICD-10 codes, then used those coefficients to compute indirect GWAS using GWAS summary statistics from the non-anonymzied ICD-10 codes. Both sources of noise (linearization and anonymization) are present in these GWAS summary statistics. Anonymization worsened the phenotype fit quality (R²) of linearized phenotypes from 0.71 to 0.37 (Figure 3B). As in the non-anonymized analysis, phenotype fit quality predicted GWAS fidelity (Pearson correlation 0.76), and only phenotypes with fit quality above 0.5 lead to GWAS fidelity above zero. Overall, GWAS fidelity was high for a large fraction of nonlinear phenotypes (65% (285/440) above 0.5, 60% (265/440) above 0.75), indicating that indirect GWAS is applicable even with the added noise from anonymization.

WebGWAS uses k = 10 to provide reasonable performance with minimal privacy risk. We evaluated how our other choices of k affect GWAS fidelity. This comparison used the same sets of phenotypes (300 Boolean combinations and 140 Phecodes) and features (340 ICD-10 codes) as before. Here, we compared GWAS fidelity of anonymized, linearly approximated data relative to non-anonymized, linearly approximated data in order to isolate the effect of anonymization. We found, unsurprisingly, that GWAS fidelity degraded as the degree of anonymity increased (Figure 3C). With k = 10, all phenotypes had a GWAS fidelity above zero, and this decreased with k across all phenotypes.

2.6 Interactive web app

WebGWAS provides two interfaces for defining phenotypes. The list interface allows users to define phenotypes using ICD-10 codes and two operators, AND and NOT (Figure 4A). The tree interface allows users to define arbitrary phenotypes using ICD codes and quantitative phenotypes, as well as 10 operators (e.g. less than, OR, etc.; Figure 4B). After a phenotype is submitted, a quality gauge appears which shows the phenotype approximation quality, which proxies GWAS fidelity. Finally, once the GWAS summary statistics are computed, an interactive Manhattan plot is shown, and a download link to the full summary statistics appears (Figure 4C).

• Download figure
• Open in new tab

Figure 4: Screenshots of the WebGWAS interface

A: The list interface allows users to pick series of ICD-10 codes that define a phenotype using “and” and “not” operations. B: The tree interface enables defining truly arbitrary phenotypes using 10 unary and binary operators, constants, and both continuous (e.g. height) and binary (e.g. ICD-10 codes) phenotypes. C: Summary information is displayed after indirect GWAS is complete. The phenotype quality gauge indicates the R² between the anonymized user-defined phenotype and the anonymized, linearized phenotype. An interactive Manhattan plot shows the results of the indirect GWAS, including clickable links to DBSNP for each variant and color highlighting to indicate which variants also reached genome-wide significance for the feature phenotypes included in the user’s definition.

2.7 Accelerating pan-biobank GWAS

In addition to WebGWAS, the underlying indirect GWAS method can be used to accelerate pan-biobank GWAS in exchange for a moderate reduction in result fidelity. To do so, we reduce the phenotype dimensionality with principal component analysis (PCA), perform direct GWAS on some fraction of the principal components (PCs), then reconstruct full-dimensional summary statistics using indirect GWAS. By reducing the number of direct GWAS that must be computed, this approach can save a large fraction of computation time.

We found that a large fraction of ICD-10 codes from the UK Biobank could be reconstructed using only a fraction of the PCs, and that performance increased as more PCs were used (Figure 5). Running this analysis in an even larger dataset, we found that GWAS sensitivity and specificity are largely preserved, even with a small fraction of the PCs (Table 1). In other words, while using fewer PCs diminishes overall GWAS fidelity (Figure 5), fidelity is largely maintained for the most extreme associations (Table 1), which are generally of greatest interest.

• Download figure
• Open in new tab

Figure 5: Indirect GWAS with PCA accelerates pan-biobank GWAS

ICD-10 phenotypes from the UK Biobank were compressed to various dimensionalities using PCA. GWAS on the reduced dimensions were used to reconstruct GWAS summary statistics for the original ICD-10 codes using indirect GWAS. Using more PCs improves both phenotype fit quality and GWAS fidelity, and it improves correspondence between the two.

We found that indirect GWAS can substantially reduce computation time for pan-biobank GWAS (Table 1). The time costs for PCA and indirect GWAS are small compared to the runtime savings from computing fewer direct GWAS. In fact, the computation time savings are roughly in line with the latent space fraction used (Table 1), meaning that roughly 50% of computation time can be saved by excluding 50% of PCs. Overall, these results show that indirect GWAS can accelerate pan-biobank GWAS in cases where some degree of approximation is acceptable, and the time savings can be immense.

4.1 Indirect GWAS method

Indirect GWAS is a method to exactly reproduce GWAS summary statistics for a phenotype defined as a linear combination of other phenotypes without directly running a new GWAS. For example, given m feature phenotypes x₁,…, x_m, define a linear combination y =∑ x _ip_i. Indirect GWAS can reproduce the GWAS results for y given only information about the features (x_i) and coefficients (p_i). To do so, our method needs to compute the following for each genetic variant: coefficient estimate, coefficient standard error, p-value, and sample si ze. We take the indirect sample size to be the minimum of the sample sizes available for each feature trait, as this is the number of samples that would be available in the equivalent direct approach. We will focus on the coefficient estimate and its standard error, as the p-value can be computed from these statistics.

Consider a regression involving a single genetic variant, a single phenotype, and N samples. Let y be the N-vector of phenotype values, g be the N-vector of genotype values (coded in any way), and Z be the N × C matrix of covariates that includes an intercept and is assumed to be full-rank.

Above, α is the C-vector of covariate fixed effects, β is the scalar effect of the genotype on the phenotype, and ε is the N-vector of errors. Covariate effects can be removed using a residual projection matrix, P = I_N − Z(Z^TZ)⁻¹Z^T. Removing covariates results in residualized phenotype , genotype , and error vectors, and it allows writing the regression more compactly.

Equation 4 is a univariate least squares regression. Let d be the appropriate degrees of freedom for this analysis. The coefficient estimate for this regression is , and the standard error is

Suppose y can be written as a linear combination of m feature phenotypes. Let X be an N × m matrix of feature phenotypes and let p be an m-vector of the projection coefficients.

Suppose that for every phenotype in X, the same genotype-phenotype regression has been performed. Let b and s be m-vectors holding coefficient estimates and standard errors, respectively, from these regressions (i.e. and ). To estimate the standard error for , we first need an estimate of , which we obtain by rearranging equation 5 using b and s.

This equation for the genotypic partial variance is written in terms of variables that can be obtained from a GWAS study—the phenotypic partial variance, the estimated coefficient, t he e stimated c oefficient st andard er ror, an d the number of degrees-of-freedom (d = N − C − 1). Every input phenotype i may result in a slightly different estimate for this value, so we use the mean across phenotypes as the final estimate for genotype partial variance . The final term to be estimated is , which can be computed using the partial covariance matrix of the feature traits, C.

Finally, we can write the coefficient and standard error estimates for y as follows:

4.2 Linear mixed models

Linear mixed models GWAS methods can also be evaluated indirectly. However, the specifics of these operations vary depending on the specifics of the method to be used. In general, linear mixed models are only treated differently when computing the phenotypic partial covariance matrix and when estimating the genotype partial variance. Mixed effect terms should be treated as additional covariates that can be residualized out before computing the phenotypic covariance matrix. This process is simplified for methods like Regenie, which produce intermediate results that include phenotypes residualized against the mixed effect terms. However, some methods (like Regenie) treat chromosomes differently (i.e. they use leave-one-chromosome-out computations), meaning that indirect GWAS needs a different phenotypic covariance matrix for each chromosome. There is no general approach to including linear mixed models, as they all use different approximations and report different intermediate results.

4.3 Validation

To confirm that indirect GWAS is mathematically correct, we compared direct vs indirect GWAS summary statistics using simulated data. First, we used PhenotypeSimulator [7] to simulate the following for 10,000 individuals: genotypes for 10,000 genetic variants, 10 phenotypes (x₁ … x₁0), and 5 covariates. Next, we simulated 10 random linear combinations of the phenotypes (e.g. y₁ = 0.5 × x₁ + … + 1.1 × x₁₀). We performed association tests for all 10 phenotypes against all genetic variants using several different methods: univariate ordinary least squares (OLS; only genotype and intercept), multivariate OLS (genotype, intercept, age, sex, 10 genetic PCs), FastGWA [10] (with covariates), SAIGE [3] (with covariates), and Regenie [11] (with covariates). Finally, we performed association tests on the random projections using both direct and indirect GWAS approaches and compared GWAS chi-squared statistics to verify that indirect GWAS produces equivalent results.

4.4 Implementation

Indirect GWAS is mathematically straightforward but challenging to implement efficiently. The method requires potentially thousands of input files, including GWAS summary statistics for every input phenotype, a phenotypic covariance matrix, and a phenotype projection matrix. As the total size of these files can be very large, we provide a high-performance implementation that uses multithreading and chunked processing to provide results efficiently. Our implementation is easily installable, written in Rust, and available on GitHub [12]. We also developed a fast implementation of the data anonymization method (MDAV) [13]. The complete code powering the web application is available on GitHub [14].

4.5 UK Biobank data processing

WebGWAS currently contains only data from White British individuals in the UK Biobank. We are working to expand this as soon as possible.

To create the input data to WebGWAS, we first removed individuals whose data were flagged for various reasons as being potentially flawed or erroneous, specifically individuals whose genetic sex was mismatched with their selfreported sex, individuals who were outliers for heterozygosity or missingness (defined by the UK Biobank using the outlier detection algorithm, abberant [15]), individuals with ten or more third-degree relatives in the UK Biobank, individuals with sex chromosome aneuploidy, and we restricted to individuals used in the computation of genetic principal components by the UK Biobank. Finally, we restricted to individuals with both binary and quantitative phenotype data available, described below.

For binary phenotypes, we considered all International Statistical Classification of Diseases and Related Health Problems, 10th revision (ICD-10) codes. We mapped and included ICD-10 codes from all six sources available in the UK Biobank: hospital inpatient ICD-9 codes, hospital inpatient ICD-10 codes, self-reported non-cancer illness codes, primary cause of death, secondary cause of death, and general practitioner outpatient diagnoses. Where necessary, we used mappings provided by the UK Biobank to convert each coding to ICD-10. Only codes with at least 100 observations were retained.

For quantitative phenotypes, we considered all quantitative data fields and excluded fields that were not relevant to this study. Specifically, we chose fields meeting the following criteria: only field data (exclude all bulk data, e.g. imaging), unisex fields only, only fields with continuous or integer values, only fields with at least 75,000 samples (per the UK Biobank data dictionary), no strictly genetic fields (e.g. exclude genetic principal components (PCs), exclude pre-computed polygenic scores), no fields dealing with sample quality control (QC) or calibration (e.g. genotyping batch, Affymetrix QC metrics), no fields dealing with home location (e.g. amount of traffic around home location), and no fields dealing strictly with data collection (e.g. the number of samples taken). We manually filtered fields using these criteria.

Many fields have both instance and array indices. Instances correspond to visits to the UK Biobank assessment center. Array indices indicate that multiple measurements or values result from a single investigation instance. We kept array indices as separate fields and took the per-individual mean across instances for each field. Finally, we transformed all quantitative phenotypes using the inverse rank normal transformation, leading to approximately standard normal distributions. Applying the above QC and filtering procedure resulted in 1378 quantitative phenotypes, 1238 binary phenotypes, and 342,350 samples. To save computation time, we restricted our analysis to HapMap3 SNPs, resulting in 1,166,145 SNPs in our final dataset.