A Fast and Accurate Algorithm to Test for Binary Phenotypes and Its Application to PheWAS

Abstract

The availability of electronic health record (EHR)-based phenotypes allows for genome-wide association analyses in thousands of traits and has great potential to enable identification of genetic variants associated with clinical phenotypes. We can interpret the phenome-wide association study (PheWAS) result for a single genetic variant by observing its association across a landscape of phenotypes. Because a PheWAS can test thousands of binary phenotypes, and most of them have unbalanced or often extremely unbalanced case-control ratios (1:10 or 1:600, respectively), existing methods cannot provide an accurate and scalable way to test for associations. Here, we propose a computationally fast score-test-based method that estimates the distribution of the test statistic by using the saddlepoint approximation. Our method is much (∼100 times) faster than the state-of-the-art Firth's test. It can also adjust for covariates and control type I error rates even when the case-control ratio is extremely unbalanced. Through application to PheWAS data from the Michigan Genomics Initiative, we show that the proposed method can control type I error rates while replicating previously known association signals even for traits with a very small number of cases and a large number of controls.

Keywords: GWAS; PheWAS; rare variants; saddlepoint approximation; single-variant test; unbalanced case-control.

Figure 1

The Projected Computation Times for Testing 10 Million Variants across 1,500 Phenotypes by Various Tests with MAFs Sampled from the MAF Distribution of the MGI Data The computation times are based on testing 10,000 simulated variants on an Intel i7 2.70GHz processor and then projecting them onto a PheWAS with 10 million variants and 1,500 phenotypes.

Figure 2

Type I Error Comparison between Score, fastSPA-2, and Firth’s Test for Variants Simulated with MAFs Sampled from the MAF Distribution of the MGI Data Type I error rates were estimated on the basis of 10⁹ simulated datasets. From left to right on the x axis, the plots consider case-control ratios 10,000:10,000 (balanced), 2,000:18,000 (moderately unbalanced), and 40:19,960 (extremely unbalanced). The top and bottom panels show empirical type I error rates at α = 5 × 10⁻⁵ and 5 × 10⁻⁸, respectively.

Figure 3

Type I Error Comparison at Different MAFs between Score, fastSPA-2, and Firth’s Test The top and bottom panels show empirical type I error rates at α = 5 × 10⁻⁵ and 5 × 10⁻⁸, respectively. From left to right, the plots consider case-control ratios 10,000:10,000 (balanced), 2,000:18,000 (moderately unbalanced), and 40:19,960 (extremely unbalanced). In each plot, the x axis represents MAF with the expected minor allele count (MAC) in parentheses, and the y axis represents empirical type I error rates. Empirical type I error rates were estimated on the basis of 10⁹ simulated datasets. 95% confidence intervals at different MAFs are also presented.

Figure 4

Empirical Power Curves for Score, fastSPA-2, and Firth’s Test The top and bottom panels consider MAF = 0.05 and 0.01, respectively. From left to right, the plots consider case-control ratios 10,000:10,000 (balanced), 2,000:18,000 (moderately unbalanced), and 40:19,960 (extremely unbalanced). In each plot, the x axis represents genotype odds ratios, and the y axis represents the empirical power. Empirical power was estimated from 5,000 simulated datasets at the test-specific α levels where their empirical type I errors were equal to 5 × 10⁻⁸.

Figure 5

Q-Q Plots for Score, fastSPA-2, SPA-2, and Firth’s Test on 5 × 10⁶ Simulated Variants with MAF Randomly Sampled from the MAF Distribution of the MGI Data The top, middle, and bottom panels show Q-Q plots in the balanced (case-control ratio = 10,000:10,000), moderately unbalanced (case-control ratio = 2,000:18,000), and extremely unbalanced (case-control ratio = 40:19,960) case-control scenarios, respectively. In each plot, the x axis represents –log₁₀ expected p values, and the y axis represents –log₁₀ observed p values.

Figure 6

Manhattan Plots for Four Different Phenotypes from MGI Data All imputed variants with MAF > 0.001 and all directly genotyped variants were included in this analysis. From left to right, the three panels show associations based on fastSPA-2, Firth’s test, and Score. The red line represents the genome-wide significance level α = 5 × 10⁻⁸.

Figure 7

Q-Q Plots for Four Different Phenotypes from MGI Data From left to right, the three panels show the Q-Q plots based on fastSPA-2, Firth’s test, and Score. The plots are color coded according to different MAF categories. 95% confidence bands are presented in gray to signify the deviance from the uniform distribution.