Secure large-scale genome-wide association studies using homomorphic encryption

A association study (GWAS) evaluates one single-nucleotide polymorphism (SNP) at a time for association to a phenotype or outcome. In the disease case/control setting, this is typically performed through a goodness-of-fit test or logistic regression, which report association odds ratios, standard errors, and P values. The results from a GWAS have two broad downstream uses: First, variants that pass a statistical threshold are reported as genome-wide significant and evaluated for functional mechanisms; second, all variants can be integrated into polygenic risk score analyses to predict phenotypes in held-out samples.

A critical challenge for GWASs is the dependence on individual-level data that typically have strict privacy requirements, creating an urgent need for methods that preserve the individual-level privacy of participants (1, 2). There are two main approaches to privacy-preserving GWASs: secure multiparty computation (MPC) and homomorphic encryption (HE). The MPC approach typically uses a protocol invented by Yao in the 1980s called the garbled circuit solution (3, 4). In this protocol, two clouds, each owned by a different hospital or corresponding to two noncollaborating servers within one hospital, hold part of the genomic data to be analyzed. The alternative approach is based on fully HE (FHE), a novel secure encryption method developed in 2008 by Gentry (5), which is much less communication intensive, and is secure even if the servers collaborate. HE allows performing secure computations over encrypted sensitive data without ever decrypting them.

Recent work has focused on secure MPC solutions to facilitate individual-level privacy-preserving GWASs (3, 6). The work of Jagadeesh et al. (3) addressed diagnosis of monogenic diseases while preserving participant privacy using MPC. Due to the communication and computationally intensive nature of the garbled circuit solution, GWASs beyond monogenic diseases were not addressed, and the patient cohort was small. Jagadeesh et al. estimated that, even for the monogenic example, garbled circuits would be at least 5,000 times faster than FHE. Cho et al. (6) followed by successfully computing a GWAS by dividing data among multiple servers and computing the GWAS via multiparty secure protocol among the servers, subsets of which are trusted not to collaborate against other servers, else privacy is lost. Here, we no longer need to resort to this trust assumption. We are successfully using HE to encrypt the genomic sequences of study participants while enabling GWAS computations without the ability to decrypt, and scaling to hundreds of thousands of samples (Fig. 1).

We implement two common GWAS techniques—the allelic chi-square test for case control differences and a logistic regression approximation (LRA) with covariates—within our HE framework. The LRA algorithm utilizes a previously proposed semiparallel approach to efficiently iterate over each genetic variant without requiring repeated likelihood maximizations (7). Our HE LRA implementation of this approach was independently tested in the iDASH 2018 secure genome analysis competition (http://www.humangenomeprivacy.org/2018/) and received the first place.* We additionally present a highly efficient chi-square test that is faster than the LRA implementation by a factor of 40$\times$ and consumes 6$\times$ less memory at the cost of excluding covariates from the model. Our HE framework provides postquantum security and is based on several advances. First, we reformulated the compute models for both the chi-square and LRA algorithms to fully benefit from ciphertext packing, enabling the parallel execution of thousands of multiplications/additions using a single homomorphic multiplication/addition. Second, we introduced two types of data encoding to minimize the number of computationally expensive key switching operations, and developed several methods for converting between the encodings homomorphically (used in the LRA solution). Third, we applied multiple plaintext approximations for the LRA model. Fourth, we developed an efficient residue number system (RNS) variant of the Cheon–Kim–Kim–Song (CKKS) HE scheme (8), which naturally supports approximate number arithmetic. Finally, we applied more than a dozen cryptoengineering optimizations.

We apply our HE framework to a real GWAS of age-related macular degeneration (AMD) (9) of 12,461 cases and 14,276 controls (restricted to self-reported Europeans) genotyped on 263,941 total markers with minor allele frequency of $>$1%. For our gold standard, we computed association statistics in the clear using full logistic regression on each variant with sex, age, and age squared as covariates. We first compared the distributions of GWAS statistics on a subset of 16,384 SNPs and 5,000 individuals evaluated by the same statistical test with/without HE, where we expect essentially perfect concordance (Fig. 2 A and B). Both the chi-square and LRA tests produced HE statistics with an R² of 1.00 to the statistics in the clear, and a replication slope of 1.00 and 0.98, respectively, indicating negligible bias (see Materials and Methods). We next compared the HE GWAS statistics to the gold standard logistic regression statistics, with any differences now arising from both the statistical assumptions and the HE (Fig. 2 C and D). The LRA again produced highly accurate HE statistics, with an R² of 1.00 and a replication slope of 0.98. The HE chi-square statistics exhibited some loss of signal relative to the gold standard but remained highly robust, with an R² of 0.96 (replication slope 0.99) primarily due to noise at nonsignificant variants. Importantly, the chi-square test odds ratios remained highly accurate and nearly unbiased (R² = 0.95; SI Appendix, Fig. S1). We confirmed that accuracy was high across all variants by computing polygenic risk scores, wherein genetic risk values are predicted for each individual as the sum of risk alleles they carry weighted by the allelic effect size (see Materials and Methods, Fig. 2 E and F, and SI Appendix, Fig. S2). Both risk scores were highly correlated with scores from the gold standard logistic test (R² $>$ 0.99) with a slight downward bias in the absolute score for the LRA (replication slope 0.95) (Fig. 2 E and F). Finally, we examined individual genome-wide significant associations reported by the original GWAS, where we again observed highly concordant results with the gold standard using both statistical tests (Table 1).

We next downsampled the data to investigate the run time and scalability characteristics as a function of SNPs and sample size (Fig. 3 and SI Appendix, Tables S1–S6). We found the LRA computation to scale linearly in the number of markers and the number of individuals. At the largest evaluated sample size of n = 15,000 and M = 16,384, computation took 1.1 h, extrapolated to 7.7 h for a GWAS of n = 100,000 (M = 16,384), and extrapolated to 234 h for a GWAS of n = 100,000 and M = 500,000 run as a single job. By comparison, the chi-square test (which requires only simple mathematical operations) was more than an order of magnitude faster than the LRA, completing the same analysis of n = 15,000 and M = 16,384 in 98 s (41$\times$ faster than the LRA). This was extrapolated to 11 min for n = 100,000 (M = 16,384) or 5.6 h for a GWAS of n = 100,000 and M = 500,000. The chi-square solution also required approximately 6$\times$ less peak RAM and thus could be run on a full-scale cohort of n = 25,000 (M = 49,152) in 8 min within the memory constraints available to us (SI Appendix, Table S4). Detailed run time characteristics including encryption/decryption are available in SI Appendix, Tables S1, S2, S4, and S5. As both the LRA algorithm and chi-square test are natively parallel over the number of SNPs, the computations can be trivially distributed to multiple nodes, with each node working with 16,384 SNPs at a time (see Materials and Methods). This implies that a GWAS of n = 100,000 and M = 500,000 could be run in 11 min on 31 nodes running in parallel.

Our HE solution for the chi-square test is faster than the state-of-the-art MPC approach of Cho et al. (6) extrapolated to 100,000 individuals and 500,000 SNPs: 5.6 h (for HE) vs. 37 h (for MPC association tests only, without quality control or population stratification analysis; Phase 3 in figure 2a of ref. 6) or 193 h (for full MPC). The accuracy of both solutions is similar. Our LRA solution has a run time of 234 h for this scenario, while the previously published MPC approach “did not yield a practical runtime for a genome-wide application of logistic regression.” Both of our solutions are fully noninteractive, produce valid odds ratios in the analyses of real data, and natively parallelize over the number of SNPs, enabling their execution in distributed computing cloud environments. Although we did not implement it here, a hybrid approach where all SNPs are evaluated with the chi-square test and then the 5% most significant SNPs are retested by the LRA could also be used to achieve the same accuracy as LRA for significant associations, requiring only 17 h (excluding ciphertext repacking overhead, which would be relatively small).

Our approach has several limitations and areas of future work. First, unlike previous work (6), our model assumes that encrypted data have been fully processed and does not perform additional quality control or genetic ancestry inference, although such methods can be easily applied preencryption. In particular, Chen et al. (10) showed that fine-scale genetic ancestry is much more accurately inferred by projection from external population reference data than by principal component analysis directly on the target samples and leads to more effective correction for population stratification. High-quality population reference data are available for all major populations, and the preencryption data can be easily projected using these references to compute ancestry covariates [requiring a simple matrix–vector product, as, for example, implemented in the PLINK score function (11)]. Second, while the chi-square test requires no parameter tuning, the LRA relies on a learning rate parameter (see Materials and Methods) that may differ by study depending on size and relationship of covariates. This can be circumvented by tuning the parameter on subsets of the data in the clear, or by comparing to parameter-free solutions such as the chi-square or linear regression results, at the cost of some additional computation. Third, our approach does not prevent the HE Compute Cloud from colluding with the GWAS coordinator to decrypt the original data, which is also true for existing MPC solutions. This problem can be addressed by adding a secret sharing protocol or using a variant of threshold HE (12) described in the next paragraph.

Extensions to a multiparty scenario are possible using threshold HE (12), a protocol where many parties cooperatively generate a common public key using their individual secret keys (“secret shares”). In this setting, the joint secret key corresponding to the common public key is never seen by any party. In GWASs, the same genotypes and phenotypes can be transmitted from multiple participants and then combined together, or genotypes and phenotypes can be separately transmitted for the same individuals from different participants and then joined together. This extension does not add substantial computation overhead to our single-party HE solution (the computation itself is performed the same way). Our work here is thus a step toward enabling analyses of sensitive phenotypes that cannot be shared between groups/institutions and individual patient participation in research studies without risk to genomic privacy.

Many of our HE improvements are general-purpose and can be applied to other application domains where similar large-scale association and regression tools are used, including phenome-wide association studies from electronic medical record data (13), discovery of predictors of treatment response in clinical trials (14), and correlative studies of multimodal data such as expression/microbiome activity (15). The tests developed here can also be extended to richer machine learning models, including decision and gradient boosted trees.