Comment on "Using Human 'Experiments of Nature' to Predict Drug Safety Issues: An Example with PCSK9 Inhibitors".

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-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

Phenotype (code)	N	OR/beta (95% CI)	p value
High cholesterol (HC269)	43,054	0.57 (0.53–0.61)	4.5 × 10−58
Statin (MED1000000001)	59,198	0.63 (0.60–0.67)	1.9 × 10−52
Simvastatin (MED1140861958)	59,046	0.69 (0.64–0.73)	1.3 × 10−28
Atorvastatin (MED1141146234)	15,983	0.60 (0.53–0.68)	2.3 × 10−14
Heart attack/myocardial infarction (HC326)	12,138	0.73 (0.65–0.81)	6.9 × 10−8
Angina (HC132)	18,473	0.79 (0.72–0.87)	3.4 × 10−7
Ezetimibe (MED1141192736)	2993	0.43 (0.30–0.61)	2.8 × 10−6
Lipitor 10 mg tablet (MED1141146138)	3788	0.50 (0.37–0.68)	7.9 × 10−6

CI confidence interval, OR odds ratio, PCSK9 proprotein convertase subtilisin/kexin type 9

Table 2

Replication of associations reported in Jerome et al. [1] in the UK Biobank as reported in the Global Biobank Engine

Phenotype (code)	N	p value
Spina bifida (HC87)	187	> 0.10
Fracture lower leg/ankle (HC339)	5881	> 0.10
Fracture foot (HC76)	1270	> 0.10
Bone disorder (HC134)	1459	> 0.10
Osteoporosis (HC22)	9135	> 0.10
Fracture vertebra/crush fracture/vertebral collapse (HC14)	294	> 0.10

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

Study and group	N (%)	Notes
Jerome et al. [1]
Total individuals in study	25,302 (100)
R46L carriers (het + hom) expected	800.2 (3.2)	Based on a = 1.594%.Ncarriers = Ntotal (a2 + 2ab)
R46L homozygotes expected	6.4 (0.025)	Based on a = 1.594%.Nhom = Ntotal (a2)
Spina bifida cases	34 (0.13)
Spina bifida R46L carriers (het + hom) expected	1.1 (0.0042)	Assuming no association
Spina bifida R46L homozygotes expected	0.0086 (0.000034)	Assuming no association
gnomAD European (non-Finnish)
Total individuals sequenced	48,202 (100)	Based on allele number = 96,404
R46L homozygotes expected	11.0 (0.023)	Based on a = 1.511%.Nhom = Ntotal (a2)
R46L homozygotes observed	10 (0.021)
gnomAD Finnish
Total individuals sequenced	10,244 (100)	Based on allele number = 20,488
R46L homozygotes expected	18.6 (0.18)	Based on a = 4.266%.Nhom = Ntotal (a2)
R46L homozygotes observed	25 (0.24)

het heterozygous, hom homozygous, N number of people in study, N expected number of R46L carriers in population, N expected number of R46L homozygotes in population

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.