Dear Editor,

We read with interest the recent work by Jerome et al. [1] using a phenome-wide association study (PheWAS) approach to predict potentially deleterious effects of proprotein convertase subtilisin/kexin type 9 (PCSK9) inhibitors. We agree with the authors’ argument that human genetics provides valuable “experiments of nature” to evaluate drug target safety, and they provide thoughtful discussion of the strengths and caveats of this approach. However, the specific associations that the authors highlight between the PCSK9 loss-of-function variant R46L (rs11591147) and spina bifida and fracture phenotypes are not statistically significant when correcting for multiple hypothesis testing and do not replicate in a larger dataset.

When scanning for statistical associations in a large database, it is critical to correct the p value threshold for multiple hypothesis testing [2]. Genome-wide association studies (GWAS) commonly use a Bonferroni-corrected p value threshold of p < 5 × 10−8, because approximately 106 independent loci are tested [3]. Similarly, PheWAS should be interpreted with a corrected p value threshold based on the number of independent phenotypes [4]. The authors do not state the total number of phenotypes tested; however, a recent study that used the same BioVU database for PheWAS on human leukocyte antigen variants [5] reported that 1368 phenotypes were tested. Assuming the same number of phenotypes was tested in this analysis, performing a Bonferroni correction on a threshold of p < 0.05 when testing a variant against 1368 phenotypes would result in a threshold of p < 3.7 × 10−5. Therefore, even the strongest non-cardiovascular result from the Jerome et al. [1] PCSK9 study (spina bifida p = 2.7 × 10−4) is statistically consistent with a chance occurrence, assuming that many phenotypes were tested (for example, p < 2.7 × 10−4 would be expected by chance 5% of the time if at least 186 independent phenotypes are tested).

To attempt to replicate the result of Jerome et al. [1] in an independent dataset, we inspected an online portal (the Global Biobank Engine) that provides PheWAS analysis of 337,208 individuals and 1837 phenotypes with genotype–phenotype data from the UK Biobank cohort [6,7,8]. Using the PheWAS report page for PCSK9 R46L [9], we inspected all associations that surpassed a phenome-wide significance threshold of p < 2.7 × 10−5 and found that all results associated PCSK9 R46L with cardioprotective phenotypes, as previously reported (Table 1). We then found the phenotypes in the UK Biobank that most closely matched the non-cardiovascular associations that Jerome et al. [1] reported, and we tested these for associations (Table 2). Although all of these phenotypes have a much larger sample size in the UK Biobank, none of the associations replicated at p < 0.10.

Table 1 Phenome-wide associations with PCSK9 R46L in the UK Biobank as reported in the Global Biobank Engine
Full size table
Table 2 Replication of associations reported in Jerome et al. [1] in the UK Biobank as reported in the Global Biobank Engine
Full size table

Jerome et al. [1] then examined the medical records of six individuals with both PCSK9 R46L loss-of-function mutation and a classification of spina bifida to define which clinical features had led to this PheWAS classification. These individuals are discussed as being homozygous and having spina bifida, which is very surprising; given the allele frequency, observing even a single homozygous individual with a rare disease such as spina bifida would be extremely unlikely unless the mutation had a highly penetrant, recessive effect on the disease (Table 3). Alternatively, if these six case studies were in fact heterozygous for PCSK9 R46L loss-of-function mutation, this would be consistent with the language the authors use in their Table 2, which reports six affected carriers; this would also be consistent with the authors’ reported odds ratio of 5.90—the mutation is more common in the disease but not strongly penetrant (six spina bifida patients carry the mutation vs. one patient expected to carry the mutation; Table 3). Clarity from the authors in this regard would be helpful to further evaluate these cases and the penetrance of the claimed effect.

Table 3 Expected prevalence of people with heterozygous and homozygous R46L genotypes in the Jerome et al. [1] study as a whole and among subsets of the spina bifida cases reported, assuming Hardy–Weinberg equilibrium and no genetic association; and observed and expected number of R46L homozygotes in two populations in the Genome Aggregation Database (gnomAD), based on Hardy–Weinberg equilibrium and no association between R46L homozygosity with severe childhood disease
Full size table

A highly penetrant recessive effect of PCSK9 R46Lis unlikely, given previous descriptions of homozygotes [10]. We looked for evidence of underrepresentation of R46L homozygotes in the Genome Aggregation Database (gnomAD) [11, 12], which excludes individuals with severe pediatric disease, and observed no deficit of homozygotes (Table 3), consistent with PCSK9 R46L homozygosity causing no severe pediatric disease.

In summary, the report by Jerome et al. [1] is an important demonstration of how PheWAS can be used to assess drug safety. However, the links they report between PCSK9, spina bifida, and fracture phenotypes are likely false positives and do not replicate in an independent database. Furthermore, evidence from case reports and population genetics does not suggest that PCSK9 R46L homozygosity causes any severe childhood disease.