A Fast and Accurate Method for Genome-Wide Time-to-Event Data Analysis and Its Application to UK Biobank

Abstract

With increasing biobanking efforts connecting electronic health records and national registries to germline genetics, the time-to-event data analysis has attracted increasing attention in the genetics studies of human diseases. In time-to-event data analysis, the Cox proportional hazards (PH) regression model is one of the most used approaches. However, existing methods and tools are not scalable when analyzing a large biobank with hundreds of thousands of samples and endpoints, and they are not accurate when testing low-frequency and rare variants. Here, we propose a scalable and accurate method, SPACox (a saddlepoint approximation implementation based on the Cox PH regression model), that is applicable for genome-wide scale time-to-event data analysis. SPACox requires fitting a Cox PH regression model only once across the genome-wide analysis and then uses a saddlepoint approximation (SPA) to calibrate the test statistics. Simulation studies show that SPACox is 76-252 times faster than other existing alternatives, such as gwasurvivr, 185-511 times faster than the standard Wald test, and more than 6,000 times faster than the Firth correction and can control type I error rates at the genome-wide significance level regardless of minor allele frequencies. Through the analysis of UK Biobank inpatient data of 282,871 white British European ancestry samples, we show that SPACox can efficiently analyze large sample sizes and accurately control type I error rates. We identified 611 loci associated with time-to-event phenotypes of 12 common diseases, of which 38 loci would be missed within a logistic regression framework with a binary phenotype defined as event occurrence status during the follow-up period.

Keywords: Cox proportional hazards regression model; GWAS; PheWAS; UK Biobank; electronic health record; saddlepoint approximation; survival analysis; time-to-event data.

Figure 1

Quantile-Quantile (QQ) Plots for Standardized Statistics and −log₁₀(P) of Score and Wald Tests (A and B) Normal QQ plots for standardized statistics (A) and QQ plots for -log10(P) of score and Wald tests (B). Standardized statistics were calculated as raw statistics $\text{S}$ divided by the estimated variance $\widehat{\text{Var}}\left(\text{S}\right)$. p values were calculated from a normal approximation. We simulated 2$\times$10⁵ replications under three event rates (ERs) of 1%, 10%, and 50%. The sample size was 4,000, and we considered common variants (MAF = 0.3, expected MAC = 2,400) and low-frequency variants (MAF = 0.01, expected MAC = 80). MAF, minor allele frequency; MAC, minor allele counts.

Figure 2

Projected Computation Time for a Genome-wide Time-to-Event Data Analysis of 20 Million Variants The projected time is based on computation time for 10,000 variants on an Intel Xeon Platinum 8176 CPU at 2.10 GHz. For example, suppose it takes $\text{α}$ h to analyze 10,000 variants, then the projected time for 20 million variants is 2,000 $\times \text{α}$ h. Solid and dashed lines represent ERs of 1% and 50%, respectively. The MAFs are randomly generated from the MAF distribution of UK Biobank, and we considered 10 covariates.

Figure 3

Empirical Type I Error Rates of SPACox, SPACox-NoSPA, Wald, Firth, and Score Tests From left to right, the plots considered five event rates (ERs) of 0.2%, 1%, 10%, 20%, and 50%. Top and bottom plots are for empirical type I error rates at $\text{α}=5\times {10}^{-5}$ and $5\times {10}^{-8}$, respectively. Sample size $\text{n}$ = 100,000. For each pair of MAF and event rate, we simulated ${10}^9$ replications.

Figure 4

Empirical Powers of SPACox, Firth, Wald, Score, and SPACC Tests When $\text{γ}$ Is Positive From left to right, the plots considered three MAFs of 0.01, 0.05, and 0.3. From top to bottom, the plots considered five ERs of 0.2%, 1%, 10%, 20%, and 50%. Empirical powers were evaluated at the significance level $5\times {10}^{-8}$. Sample size $\text{n}$ = 100,000. For each pair of MAF and event rate, we simulated 1,000 replications.

Figure 5

Manhattan Plots for 12 Phenotypes from UK Biobank Manhattan plots were based on p values calculated from the SPACox method. The red line represents the genome-wide significance level $\alpha =5\times {10}^{-8}$.

Figure 6

p Values of SPACC and SPACox for 38 Highlighted SNPs from UK Biobank We highlight 38 loci that are significant on the basis of SPACox but not significant on the basis of SPACC. The red lines represent the genome-wide significance level $\alpha =5\times {10}^{-8}$.