Genetics of 35 blood and urine biomarkers in the UK Biobank.

Clinical laboratory tests are a critical component of the continuum of care. We evaluate the genetic basis of 35 blood and urine laboratory measurements in the UK Biobank (n = 363,228 individuals). We identify 1,857 loci associated with at least one trait, containing 3,374 fine-mapped associations and additional sets of large-effect (>0.1 s.d.) protein-altering, human leukocyte antigen (HLA) and copy number variant (CNV) associations. Through Mendelian randomization (MR) analysis, we discover 51 causal relationships, including previously known agonistic effects of urate on gout and cystatin C on stroke. Finally, we develop polygenic risk scores (PRSs) for each biomarker and build 'multi-PRS' models for diseases using 35 PRSs simultaneously, which improved chronic kidney disease, type 2 diabetes, gout and alcoholic cirrhosis genetic risk stratification in an independent dataset (FinnGen; n = 135,500) relative to single-disease PRSs. Together, our results delineate the genetic basis of biomarkers and their causal influences on diseases and improve genetic risk stratification for common diseases.

Introduction

Serum and urine biomarkers are frequently measured to diagnose and monitor chronic disease conditions. Understanding the genetic predisposition to particular biomarker states, and the factors that confound them, may have implications for disease treatment. While the genetics of some biomarkers have been extensively studied, most notably lipids[1,2,3], glycemic traits[4-6], and measurements of kidney function[7-9], the genetic basis of most biomarkers has not been queried in large population-scale datasets.

To this end, UK Biobank (UKB) has performed laboratory testing of >30 commonly measured biomarkers in serum and urine on a cohort of >480,000 individuals with extensive phenotype and genome-wide genotype data, including the unrelated individuals in this study (Supplementary Figure 1)[10].

Here, we 1) performed a systematic analysis of the genetic architecture and detailed fine-mapping of biomarker-associated loci in 363,228 individuals including protein-altering, protein-truncating (PTV), non-coding, human leukocyte antigen (HLA), and copy number variants; 2) built phenome-wide associations for implicated genetic variants; 3) evaluated causal relationships between biomarkers and 40 medically relevant phenotypes; and 4) constructed polygenic prediction models (Figure 1).

Figure 1.

Schematic overview of the study.

We prepared a dataset of 35 serum and urine biomarkers from 363,228 individuals in UK Biobank. We analyzed the genetic basis of these biomarkers, assessed their relationship to medically relevant phenotypes, and generated predictive models of disease outcomes from genome-wide data.

Results

Biomarker phenotype distributions

We first examined the consistency of the biomarker measurements[10,11]. After adjusting for statin usage (Supplementary Table 1a–c), we fit a regression model with multiple covariates (see Methods). For each biomarker, we measured the proportion of phenotypic variance explained by these covariates; this ranged from 1.7% (Rheumatoid factor) to 90% (Testosterone) depending on the biomarker (Supplementary Figure 2a–c, Supplementary Table 2). We evaluated body mass index as a confounder in associations, and there were minimal differences in genetic effects under this model (Supplementary Tables 3–4). Taking all the 35 lab phenotypes together, we recover several previously estimated phenotype correlations (Figure 2)[12,13].

Figure 2.

Genetics of 35 biomarkers. (top left inset)

Correlation of phenotypic (lower triangular matrix) and genetic (upper triangular matrix) effects plot between the 35 lab phenotypes, estimated using LD Score regression. The absolute heritability estimates with standard errors are in Supplementary Table 11a. (main panel) Fuji plot of lab phenotypes across the six categories provided by UK Biobank and genetic variant associations shown for LD independent variants with meta-analysis p < 5 × 10−9. Large-effect protein-truncating and protein-altering variants (labeled when abs(beta) >= 0.1 standard deviation [SD]) annotated with the category of association displayed (colored fill boxes) and highlighted if the loci were not previously reported in the comparison studies (Methods). Pleiotropic association and trait-specific association are shown by different sized circles. The p-values were from two-sided tests and were not corrected for multiple hypothesis testing.

Genetics of biomarkers

We performed association analysis between directly genotyped and imputed autosomal genetic variants, copy number variations (CNVs), and HLA allelotypes and 35 biomarkers in the unrelated individuals in UK Biobank across 5 population groups (N = 318,953 for White British, 23,582 for non-British White, 6,019 for African, 7,338 for South Asian, and 1,082 for East Asian) followed by meta-analysis of all but East Asian populations (N in meta-analysis = 355,891, Methods, Figure 2, Supplementary Figure 3). We stratified the genetic variants into three bins: 1) protein-truncating (27,816), 2) protein-altering (87,430), and 3) synonymous and non-coding variants (minor allele frequency [MAF] > 0.1% and INFO score > 0.3, imputed variants present in Haplotype Reference Consortium [HRC], 9,444,561[14]) (Figure 2). Comparison of effect sizes estimated across 42 other previously published study cohorts for 25 of the biomarkers showed overall high agreement (Supplementary Figure 4, Supplementary Table 5). This was true when comparing to previous studies of lipids[1,2,15,16], glycemic traits[17,18], kidney function tests[19,20], liver function tests[17], and other biomarker measurements[21,22].

We adjusted the nominal association p-values for multiple hypothesis testing and identified over 10,000 significant associations (Bonferroni-corrected meta-analysis p < 5 × 10−9 for assayed and imputed variants; Bonferroni-corrected p < 1 × 10−6 for non-rare [MAF > 0.1%] CNVs and CNV burden test for 23,598 genes; and Benjamini-Yekutieli [BY] adjusted p < 0.05 for HLA alleles, Methods, Supplementary Figure 5, Supplementary Tables 6–10). Linkage disequilibrium (LD) Score intercepts for single-variant association results were between 0.999 and 1.137 for all 35 phenotypes, consistent with anthropometric traits in UKB and suggesting that population structure in our analysis is well-controlled[23] (Supplementary Table 11a).

Global and local heritability of biomarkers

To characterize the heritability of the 35 biomarkers we first applied LD Score regression[24] to stratify heritability into 10 tissue types and 53 general genomic features (e.g. coding variants and regulatory variants) and further applied the Heritability Estimator from Summary Statistics (HESS)[25,26]. We found that both LD Score regression and HESS indicate common SNPs explain substantial heritability of some but not all biomarkers (0.6% [Lipoprotein A, also referred to as lipoprotein(a)] to 23.9% [IGF-1] using LD Score regression and 3.2% [Microalbumin in urine] to 57% [Total bilirubin] using HESS across the studied continuous phenotypes, Supplementary Tables 11a,b). Estimates were lower in LD Score regression than HESS for traits with lower polygenicity (e.g. Lipoprotein A, h2ldsc = 0.6% and h2HESS = 24%), as LD Score regression estimates polygenic heritability[24]. We compared the polygenicity of all 35 biomarkers by computing the fraction of total SNP heritability attributable to loci by the top 1% of SNPs. We found that three biomarkers have more than 50% of the SNP heritability explained by the top 1% of loci (Lipoprotein A 67.7%, total bilirubin 60.9%, and direct bilirubin 57.5%) while the remaining 32 phenotypes show patterns of moderate to high polygenicity (Supplementary Table 11b).

Associated variants prioritize therapeutic targets

We found 58 (43 rare, MAF < 1%, and 55 not reported in comparison study, Methods) PTV associations and 1,323 (306 rare, 1,079 not reported in comparison studies) protein-altering variant (PAV) associations outside the major histocompatibility complex (MHC) region (hg19 chr6:25,477,797–36,448,354; meta-analyzed p < 5 × 10−9). We found 19 non-MHC PTV associations (17 rare [MAF < 1%]) with large estimated biomarker-lowering effects (>0.1 sd) and 26 (24 rare) with biomarker-raising effects (>0.1 sd) across 31 (27 rare) PTVs and at least one biomarker phenotype, where the same PTV may have both increasing and decreasing associations across different biomarkers (Figure 2, Supplementary Table 6). Similarly, there were 240 (161 rare) and 182 (125 rare) non-MHC PAV associations with large estimated lowering and raising effects (>0.1 sd) across 241 (179 rare) PAVs and at least one biomarker phenotype, respectively (Figure 2, Supplementary Table 7). To assess whether the variants associated with biomarkers impact medically relevant phenotypes, we performed a phenome-wide association analysis (PheWAS) across 166 traits in UK Biobank, compared our findings with previously published literature, and sought independent replication in the FinnGen R2 cohort (Supplementary Tables 12–13, Methods). We found 57 phenotype associations (33 and 24 for increasing and decreasing disease risk, respectively) across 26 medically relevant phenotypes for 2 PTVs and 31 PAVs (p < 1 × 10−7), of which 31 associations were previously reported and 26 were novel (Supplementary Tables 13a, Methods).

For eight cardiovascular biomarkers (Supplementary Table 4a), we identified a stop-gain variant in PDE3B with documented protection against high cholesterol and a range of effects on increasing HDL cholesterol and Apolipoprotein A (0.40, 0.27 sd) and decreasing triglycerides and Apolipoprotein B (0.43, 0.27 sd)[2,27]; a stop gain variant in ANGPTL8, where we replicated a previously-reported effect on HDL cholesterol (0.06 sd in our dataset) and discovered a triglyceride-lowering effect (0.06 sd)[28]; two PTVs in LPA with lowering effects on Lipoprotein A levels (0.37, 0.42 sd), of which one is known to be associated with decreased risk of coronary artery disease (p = 3 × 10−11; OR = 0.89 [95% CI 0.86, 0.92])[29]; a 0.2% MAF missense allele in ACACB associated with LDL, triglyceride, ApoB, and alkaline phosphatase[30]; two independent missense alleles in PLA2G12A with increasing effects on triglycerides, sex-hormone binding globulin (SHBG), and testosterone, and lowering effects on HDL cholesterol, ApoA, and HbA1c levels (Supplementary Table 6); a splice region variant in CPT1A, with lowering effects on triglycerides; and a missense variant in PCSK6 with ApoB- and LDL-lowering effects (Supplementary Table 7).

For seven liver biomarkers (Supplementary Table 4a), we found a 0.05% MAF inframe deletion in GOT1 with a lowering effect on aspartate aminotransferase (2.6 sd); a 0.1% MAF missense allele in SLC30A10 with increasing effects on alanine and aspartate aminotransferases; four missense alleles in GPT with alanine aminotransferase lowering effects; a missense variant in ABCB4 with increasing effect on alanine aminotransferase and increased risk of gallstones in UK Biobank (p = 1.2 × 10−8, OR = 1.38 [95% CI: 1.23, 1.38]); an allelic series of 3 missense variants in SERPINA1 with pleiotropic increasing effects on albumin, aspartate aminotransferase, direct bilirubin, and gamma glutamyltransferase, and lowering effects on AST to ALT ratio, with one of these missense alleles associated with increased risk of gallstones (p = 8.1 × 10−17, OR = 1.36 [95% CI: 1.27, 1.47]) and cholecystitis (p = 1.6 × 10−8, OR = 1.26 [95% CI: 1.16, 1.37]) in UK Biobank; and two missense alleles in DGKD, with raising and lowering effects, respectively, on direct and total bilirubin (Supplementary Tables 7, 13a).

For 12 renal biomarkers (Supplementary Table 4a), we found a PTV in COL4A4 associated with an increasing effect on microalbumin in urine (0.77 sd) and an increased risk of kidney disease (p = 6.7 × 10−13, OR = 6.9 [95% CI: 4.06, 11.6]) in UK Biobank, which is defined using a combination of hospital in-patient record (ICD-10 codes: Q60 [Renal agenesis and other reduction defects of kidney], and its sub-concepts) and self-report kidney diseases (coded as 1405 [other renal/kidney problem] in UK Biobank)[31]; a frame-shift variant in SLC22A2 with strong lowering effects on eGFR (0.52 sd) and increasing effect on creatinine (0.52 sd); a stop-gain variant in SLC22A11 with raising effects on urate (0.14 sd; Supplementary Tables 6, 13a); a 0.1% rare missense allele in SLC34A3 with strong eGFR and phosphate lowering and serum creatinine, Cystatin C, and urea raising effects; missense alleles in SLC6A19, LRP2, ALDOB, and SLC7A9, and two missense variants in SLC25A45, all associated with creatinine lowering and eGFR raising, among other examples (Supplementary Table 7). Notably, the majority of these genes are known to have high expression levels in renal tissue[32].

For three bone and joint biomarkers (Supplementary Table 4a), we found an allelic series of two frame-shift variants and a missense variant in GPLD133, in addition to an allelic series of missense variants in ALPL. Similarly, we found an allelic series in CASR that is associated with both calcium increasing and lowering effects (Supplementary Tables 6, 7).

For glucose and HbA1C (biomarkers for diabetes, Supplementary Table 4a), we found a known missense variant association in ANKH (−0.11 and −0.17 sd for glucose and HbA1C, respectively), which we also replicated the documented protective effects to diabetes (p = 1.2 × 10−8, OR = 0.66 [95% CI: 0.57, 0.76]). We also found a splice-donor variant in RHAG that is strongly associated with lower HbA1c (0.80 sd) and allelic series containing 4 missense variants each in G6PC2 and TMC8 (Figure 2, Supplementary Tables 6, 7).

For three hormone biomarkers (Supplementary Table 4a), we found a PTV in ADH1C, MSR1, and NUBP2 affecting serum IGF-1 levels, and an allelic series including the hepatocyte growth factor genes HGFAC, HGF, and HNF4A with effects on SHBG. Among those, we identified novel associations with HNF4A alleles: a missense variant with MAF 0.02% was associated with increased risk for diabetic eye disorders (p = 3.1 × 10−8, OR = 9.60 [95% CI: 4.30, 21.4]) and diabetes (p = 4.7 × 10−8, OR = 3.8 [95% CI: 2.34, 6.09]) and another missense variant (MAF 3.1%) was associated with increased risk for cholecystitis (p = 2.2 × 10−13, OR = 1.27 [95% CI: 1.22, 1.38] in UK Biobank, and also replicated in FinnGen R2, p = 2.9 × 10−17, OR = 1.46 [95% CI: 1.34, 1.60]) (Figure 2, Supplementary Tables 6, 7, 13a).

These results suggest that the genetic underpinning of biomarker levels could aid in prioritizing and better understanding the mechanisms of disease-associated variants.

CNVs and HLA allelotypes influencing biomarkers

Copy number variations (CNV) constitute a significant fraction of all base pair differences between individuals. We found 13 unique associations across 10 individual CNVs (Bonferroni p < 1 × 10−6, MAF > 0.01%, Supplementary Table 10a)[34]. We performed aggregate rare (MAF < 0.1%) CNV burden tests, pooling CNVs in each gene, for 23,598 genes. We found 29 gene-level associations (Bonferroni p < 1 × 10−6; Supplementary Table 10a) including a burden of rare CNVs overlapping HNF1B associated with serum urea, eGFR, creatinine, and cystatin C (the least significant p = 8.8 × 10−13) estimated to have large effects (beta = 0.77, −0.90, 0.93, 0.98 sd, respectively; Supplementary Figure 6a). Previous studies have associated mutations in HNF1B with maturity onset diabetes of the young (MODY5) and altered kidney function[35].The rare CNVs overlapping HNF1B were associated with chronic kidney disease (p = 1 × 10−7; OR = 4.94, SE = 0.30; Supplementary Figure 6a)[36,37] in a diabetes-dependent fashion (Supplementary Table 10b). We found a rare duplication in the CST3 gene associated with increased levels of cystatin C, the protein it encodes, of opposite effect to a rare PTV at the same locus (Supplementary Figure 6b). These results highlight the value of CNV analysis with potentially large effects on lab measurements.

To identify HLA allelotype associations that are not driven by pervasive LD structure in the HLA region, we applied Bayesian model averaging (Methods) to the significant allelotype-trait pairs (BY adjusted p-value < 0.05). We found 58 associations across 28 biomarker traits and allelotypes (Supplementary Table 9).

Fine-mapping of common associated variants

To nominate potentially causal variants at loci with common (MAF > 1%) variant associations, we performed fine-mapping analysis. Specifically, focusing on the White British summary statistics, we applied FINEMAP[38,39]. From over 9,000 biomarker-associated loci, we identified 27,853 distinct signals in 5,363 regions across 35 traits. In the identified credible sets, 17,696 signals were fine-mapped to 50 or fewer variants with posterior probability of including the causal variant > 0.99; at 2,547 biomarker-associated loci, we resolved the signal to a single nominated causal variant (Figure 3a, Supplementary Table 14a). Moreover, we identified 3,374 unique trait-variant associations with a posterior probability > 0.99 of being the causal variant (Figure 3b, Supplementary Table 14a). These explain between 0% (urine potassium) and 48% (Lipoprotein A) of the residual trait variance (Supplementary Table 14b).

Figure 3.

Summary of fine-mapped associations across 35 biomarker traits.

(a) FINEMAP analysis summary. (top) The number of identified distinct association signals (color gradient from green to blue) in each region with at least one genome-wide significant (UK Biobank meta-analysis p < 5 × 10−9) association and the number of regions are shown, such as a single signal at 33 regions and two to forty signals at 5,330 regions across 35 traits. (bottom) The number of identified candidate causal variants in the credible set with >= 99% posterior probability (color gradient from green to blue) and the number of signals are shown, such as 2,547 signals were mapped to a single variant in the credible set across 35 traits. (b) Breakdown of the number of fine-mapped associations with posterior probability greater than 0.95 or 0.99 across all biomarkers. Orange, posterior greater than 0.99, green, posterior between 0.95 and 0.99. The total variance explained for each trait is shown and in Supplementary Table 14b. (c) Allelic series showing combined missense, non-coding, and rare copy number variants at the SLCO1B1/SLCO1B3 on total bilirubin levels. Copy number variants annotated below axis and SNPs and short indels annotated above the axis. (d) Pleiotropic effects of fine-mapped rare coding (rs114303452, left) and common non-coding (rs59950280, right) variants at the HGFAC locus. Darker colors of purple indicate more significant associations. The p-values were from two-sided tests and were not corrected for multiple hypothesis testing. The error bars represent standard deviations.

Glycemic trait fine-mapping

We discovered fine-mapped associations for glycemic traits, including multiple variants at the TGFB1/AXL locus; rare missense variants in PFN1 and GYPC (previously implicated in a small GWAS of Mexican Americans)[40]; an intronic variant at the cytokine receptor IL6R; a downstream variant at VEGFA41; a missense variant at HFE, the gene responsible for hemochromatosis[42,43]; and an intronic variant at CD33 and 3’UTR variant at CD36 (Supplementary Table 14a). CD36 encodes for a well-studied fatty acid receptor and biomarker for type 2 diabetes (T2D)[44,45], and CD33 levels are known to be perturbed in T2D cases[46].

Allelic series at the SLCO1B locus

We discovered several alleles implicated in the genetic control of bilirubin levels at the SLCO1B locus (Figure 3c). We find several heterozygous deletion events, and single-nucleotide variants that we fine-mapped to two main signals: a missense variant in SLCO1B1 (rs34671512, marginal beta = −0.11 sd, p = 1.25 × 10−95) and a non-coding association in an intron of SLCO1B3 (rs11045598, marginal beta = 0.076 sd, p = 1.31 × 10−139). Despite two PTVs (one in SLCO1B1 and one in SLCO1B3) at the locus, neither had a conditionally independent effect on bilirubin levels. The diversity of variant types at this critical bilirubin and drug transporter suggest that large-effect loci can harbor variants with multiple independent genetic mechanisms contributing to their trait associations.

HGFAC pleiotropy

We scanned for loci with large effects across multiple biomarkers. The most prominent of these is HGFAC, the gene encoding hepatocyte growth factor activator. At this locus, we discovered two independent fine-mapped variants, rs114303452 (a missense variant with MAF = 1%) and rs59950280 (a non-coding variant with MAF = 34%). These two variants show significant associations with a number of diverse biomarker traits, including lipids, IGF-1, albumin, and calcium (Figure 3d). In addition, rs114303452 has been previously associated with serum HGF levels[47], supporting the role of HGFAC in control of a number of other serum biomarkers through regulation of hepatocyte growth factor.

Targeted phenome-wide association analysis

We conducted PheWAS of the fine-mapped imputed variants across 166 UK Biobank phenotypes and identified 14 and 263 coding and non-coding associations, of which 109 were not previously reported in literature (p < 10−7, Supplementary Table 12, 13b–c, Methods). For example, a common (MAF = 33%) intronic variant in DPEP1 has protective effects against skin cancers (OR = 0.88 [95% CI: 0.86, 0.90], 0.81 [0.77, 0.84], and 0.89 [0.87, 0.91] for skin cancer, malignant melanoma, and non-melanoma skin cancer, respectively), with replication in FinnGen R2 (p = 3.1 × 10−5, OR = 0.81 [95% CI: 0.74, 0.90] for malignant neoplasm of skin). An allelic series of two intronic variants in ABCG2 identified with increasing and lowering urate level associations that have risk-increasing (p = 2.8 × 10−67, OR = 1.38 [95% CI: 1.33, 1.44]) and protective (OR = 0.72 [95% CI: 0.69, 0.74]) associations with gout, respectively. Both of these associations with gout are also replicated in FinnGen R2 (p = 6.3 × 10−6, OR = 1.25 [95% CI: 1.13, 1.37] and p = 8.4 × 10−5, OR = 0.84 [95% CI: 0.78, 0.92]). Those two variants (r2 = 0.47) have low linkage with a known common protein-altering variant in ABCG2 (r2 = 0.22 and 0.11 in UKB White British for rs2231142 [Q141K]) which contributes to risk of gout[48]. These results indicate that variants with effects on biomarkers may have pleiotropic effects across medically relevant phenotypes.

Causal inference

Given the relevance of several of the biomarkers studied to disease conditions we estimated causal effects of biomarker levels on 40 medically relevant phenotypes (including 32 diseases; Supplementary Table 15) using two-sample Mendelian Randomization with the genome-wide significant variants for each biomarker as instrumental variables[49-52] (see Methods). We identified 51 significant causal relationships at FDR of 5% (Figure 4a, Supplementary Table 16). Many of these and their causal effects are well-described. We found genetic evidence supporting the protective effect of sex hormone binding globulin on diabetes (0.7 OR/SD) consistent with existing reports[53-54], of ApoA on fasting glucose (0.84 OR/SD, FDR-adjusted p = 0.02)[55-57], as well as an increasing effect of alanine aminotransferase levels on diabetes (1.53 OR/SD, FDR-adjusted p = 0.0018)[58,59]. There was a consistent effect of cystatin C on stroke (1.2 OR/SD, FDR-adjusted p = 8.7×10−4 for any stroke and 1.21 OR/SD, FDR-adjusted p = 2.8×10−3 for ischemic stroke)[60,61]. Finally, both HDL and ApoA were associated with increased risk of age-related macular degeneration[62], as was cystatin C[63-65].

Figure 4.

Causal inference, transferability of polygenic risk scores, and complex trait association in polygenic risk tails.

(a) Mendelian Randomization estimates causal links between biomarkers (blue nodes) and selected complex traits (red nodes). Association arrows are drawn based on effect direction (red decreasing, blue increasing). Associations were adjusted for FDR 5% cutoff across all tests (Methods, Supplementary Table 16). Edge width is proportional to the absolute causal effect size (log odds per standard deviation). (b) Summary of prediction accuracy of the snpnet polygenic scores across traits, evaluated on a held-out test set in White British as well as other 4 populations in UK Biobank. (c) (x-axis) Biomarker polygenic risk scores for the top 1%, top 10%, bottom 1%, and bottom 10% of individuals and their association to different diseases in UK Biobank, represented as the odds ratio of the disease in this group relative to the 40–60% quantiles. Traits without rows did not have any outcomes with FDR-adjusted significant associations.

We estimated a causal protective effect of testosterone on inflammatory bowel disease (0.70 OR/SD, FDR-adjusted p = 3.86×10−3)[66,67], and a protective effect of urate on breast cancer risk (0.87 OR/SD, FDR-adjusted p = 0.033)[68].

Polygenic prediction of biomarkers

The vast size of the UK Biobank cohort affords the opportunity to build predictive polygenic risk models of biomarkers from genotype data alone[69]. We constructed PRS for all 35 biomarkers using batch screening iterative lasso (BASIL) implemented in the R snpnet package[70,71]. Specifically, we split the White British individuals into 70% training, 10% validation (to identify the optimal sparsity parameter), and 20% test sets and evaluated the predictive performance (R2) in the held-out test set (n = 63,818) as well as in 4 populations in UK Biobank (see Methods). We found the mean predictive performance relative to the White British test set for these 4 populations were 93%, 70%, 51%, and 24%, respectively, suggesting these polygenic models have limited generalizability across populations (Figure 4b, Supplementary Figure 7, Supplementary Table 17)[72]. As an external validation, we found that the PRSs had high portability to self-identified white individuals from the MESA cohort (Supplementary Table 18)[73].

Multiple regression with PRSs

We hypothesized that the 35 biomarker PRSs may improve the prediction of higher-level traits and diseases in combination with the PRS for the trait or disease itself. To this end, we constructed multi-PRS models for traits by using multiple regression to predict the trait or disease from a) its own PRS, b) the PRSs for each of the 35 biomarkers, and c) the covariates age, sex, and principal components (Methods).

We selected disease endpoints for multi-PRSs analysis by considering the enrichment of disease prevalence at the tails of the distribution of the single-trait biomarker PRSs (Figure 4c, Supplementary Table 19, Supplementary Figure 8). We focused on traits with three or more associated biomarkers (Supplementary Figure 9), as we reasoned these would benefit most from the combination of multiple biomarker PRSs.

For chronic kidney disease, the multi-PRS stratified individuals according to disease status better than the snpnet PRS (Figure 5a–b, Supplementary Table 20). In contrast, the myocardial infarction snpnet PRS was equally stratifying as compared to the multi-PRS, with both explaining a substantial portion of trait heritability (AUC 0.58–0.59; Figure 5c). This trend held after including additional existing polygenic scores for type 2 diabetes and myocardial infarction as well (AUC 0.594 and 0.611 respectively; Supplementary Figure 10, Supplementary Tables 20–23). This suggests that the genetic basis of myocardial infarction, as previously reported[74], already captures the majority of the genetic component of serum lipids and other biomarkers. Similar weak effects of biomarkers were estimated in hypertension, angina, and gallstones, while alcoholic cirrhosis, gout, type 2 diabetes, and heart failure were better predicted with multi-PRS models (Figure 5c, Supplementary Table 20). Improved predictions relied on relevant and variable biomarkers across these traits (Supplementary Figure 11), including eGFR, creatinine, cystatin C, and bilirubin for CKD; creatinine, bilirubin, total and LDL cholesterol, cystatin C, and eGFR for heart failure; and bilirubin, GGT, eGFR, HDL cholesterol, and IGF-1 for alcoholic cirrhosis.

Figure 5.

Multiple regression with biomarker polygenic scores improve prevalent and incident disease prediction.

(a) (x-axis) quantiles of polygenic risk score, spaced to linearly represent the mean of the corresponding bin of scores. (y-axis) Prevalence of chronic kidney disease (n = 2,780 cases and n = 89,409 total, defined by verbal questionnaire and hospital in-patient record ICD code data) within each quantile bin of the polygenic risk score. Error bars represent the standard error around each measurement, and individuals evaluated are held-out European-ancestry individuals in UK Biobank. (b) ROC curve with AUC for chronic kidney disease, comparing the snpnet-derived polygenic score to a multi-PRS model trained across biomarkers as well. Individuals evaluated are held -out European-ancestry individuals in UK Biobank. (c) AUC-ROC estimates for prediction of 10 disease outcomes in a held-out test set of the UK Biobank. Diabetes was run using both a strict definition (excluding from control individuals with HbA1c < 39) and the complete sample (Methods). (d) Hazard ratios for the incidence of type 2 diabetes (n = 17,519), chronic kidney disease (n=3,058), myocardial infarction (n=7,913), heart failure (n = 13,965), gout (n = 1,936), gallstones (n = 11,629), and cirrhosis (n=845) in FinnGen using the standard single-disease PRS trained on UK Biobank using snpnet versus the multi-PRS including both biomarker PRSs and the trait PRS. The strict definition of type 2 diabetes is shown. Error bars represent 95% confidence intervals and points represent mean hazard ratio estimates.

Encouraged by these findings, we evaluated the potential of these improved polygenic scores in identifying disease cases by applying both trait-specific PRS and combined PRS in an independent replication cohort, FinnGen (R3, n = 135,500, Supplementary Tables 24–26). Here, we found evidence that the combination of PRS increased the effect size in chronic kidney disease (hazard ratio (HR) = 0.99, p = 0.46 for snpnet PRS and hazard ratio = 1.12, p = 2.09 × 10−10 for multi-PRS, Figure 5d, Supplementary Figure 12, Supplementary Table 24), type 2 diabetes (HR = 1.37, p < 2 × 10−16 for snpnet PRS and HR = 1.49 for multi-PRS), gout (HR = 1.39, p < 2 × 10−16 for snpnet PRS and HR = 1.58 for multi-PRS), heart failure (HR = 1.01, p = 0.38 for snpnet PRS and HR = 1.08, p < 2 × 10−16 for multi-PRS), and alcoholic cirrhosis (HR = 0.97, p = 0.35 for snpnet PRS and HR = 1.18, p = 1.04 × 10−6 for multi-PRS, Supplementary Table 24). Similar results to UK Biobank were found in models including existing polygenic scores (Supplementary Table 26, Supplementary Figure 13), with the integrated type 2 diabetes model, including both pre-existing PRS and biomarker PRS, resulting in 1.67 HR change per standard deviation. This suggests that multiple regression of polygenic risk for biomarkers might capture multiple underlying disease states and/or underlying causes, and that these multiple states are predictive of disease.

Discussion

Using data from 35 biomarkers in ~363,000 UKB samples, we provide an assessment of genetic associations with biomarker levels, the relevance of these associations in disease phenotypes, and their utility in risk stratification.

Protein-altering variants that modify biomarker levels and disease risk can provide in vivo validation of therapeutic targets[75,76]. Here, we found multiple protein-altering variants that directly implicate genes associated with the studied biomarkers, and we hypothesize that some of these genes may provide potential therapeutic targets.

To assess the translatability of our findings, we built predictive models aggregating trait PRS with those of the biomarkers, improving the predictive accuracy of multiple disease outcomes both overall and especially at the extremes of genetic risk. Given that biomarker values are already routinely collected in structured data formats, we anticipate that the multi-PRS methods could inform clinical practice in the coming years, as a larger fraction of the population is genotyped and sequenced.

In addition to the discovery of multiple individual loci and candidate causal variants, we can also draw some general conclusions across the traits evaluated with our multi-PRS models. Traits and diseases were predicted best when they had individually predictive biomarkers and a complex etiology (e.g. chronic kidney disease), but underpowered genetic studies. We believe that a large number of disease cases is typically most useful in developing well-powered models, as it helps both with the baseline polygenic score and fitting of the multi-PRS components. Further exploration of the conditions where multi-PRS models perform particularly well is an area of future study.

Numerous limitations to this work are present. We assigned individuals to ancestry groups based on self-reported ethnicity categories and the top two principal components of the genotype matrix. We included many technical covariates in order to reduce bias in the measurements of the biomarkers, but doing so has the potential to reduce power. We fine-map based on imputed genotypes and summary statistics, and both reduce power to detect true causal variants. In addition, the large and complex linkage present at some loci, including notably the LPA locus, might result in spurious fine-mapped and rare coding variant associations, though conditional analyses (e.g. of a rare coding variant in SLC22A2) are inconclusive (Supplementary Figure 14). Similarly, causal inference using individual-level data[77] can increase power and reduce bias, and we recommend it for future studies. Lastly, we anticipate that including other genetic risk scores will fit well into the multi-PRS framework to further improve prediction of common complex disease.

The genome-wide resource made available with this study, including the association summary statistics, fine-mapped regions, and polygenic prediction models (Data availability)[78], provides a starting point for causal mapping of genetic variants affecting the 35 biomarkers and their relevance to medical phenotypes. These results highlight the benefits of direct measurements of biomarkers in population cohorts for interpreting the genetic basis of biomarkers and improved prediction of multiple common diseases.

Online Methods

Genotype and phenotype data in UK Biobank

We used genotype datasets from the UK Biobank (release version 2 for the directly genotyped variants and the imputed HLA allelotype datasets and release version 3 for the imputed genotype dataset), the copy number variation dataset[34], and the hg19 human genome reference for all analyses in the study[10]. To minimize the variability due to population structure in our dataset, we restricted our analyses to unrelated individuals based on the following four criteria reported by the UK Biobank in the sample QC file, “ukb_sqc_v2.txt”: 1) used to compute principal components (“used_in_pca_calculation” column); 2) not marked as outliers for heterozygosity and missing rates (“het_missing_outliers” column); 3) do not show putative sex chromosome aneuploidy (“putative_sex_chromo- some_aneuploidy” column); and 4) have at most 10 putative third-degree relatives (“excess_relatives” column).

Additionally, we used the “in_white_British_ancestry_subset” column in the sample QC file as a part of the population definition as shown below.

We used a combination of self-reported ethnicity (UK Biobank field ID: 21000) and principal component analysis and analyzed 5 subpopulations in the study: self-identified White British (n = 337,151 individuals), African (6,498), East Asian (1,772), South Asian (7,962), and self-identified non-British White (24,909). We first used the genotype principal components (PCs) of the genotyped variants from the UK Biobank and defined thresholds on PC1 and PC2 and further refine the population definition (described in Supplementary Note). We subsequently focused on a subset of individuals with non-missing values for covariates and biomarkers as described below.

Variant annotation and quality control

Detailed information on variant annotation and quality control is described in Supplementary Note.

Biomarker phenotype definition

Phenotype and covariate quality control excluded rheumatoid factor and estradiol from further analyses, and fasting glucose (available for 17,439 self-reported fasting individuals) was used as a phenotype-level quality control for the glucose measurements -- throughout the text, “glucose” refers to glucose levels adjusted for fasting time rather than the GWAS among only fasting individuals (self-report of more than 7 and less than 24 hours of fasting, n = 17,439) unless otherwise noted. We focused on 32 biomarkers for genetic analysis and also defined three derived phenotypes, estimated glomerular filtration rate (eGFR), non-albumin protein, and AST to ALT ratio, for a total of 35 biomarkers (Supplementary Table 4a). The eGFR measure is an indicator of renal function and is defined by the CKD-EPI equation[79]. We defined non-albumin protein as the difference between the total protein and albumin. Then after applying covariate correction (see Covariate correction below, Supplementary Table 4b), we additionally defined the AST to ALT ratio as the difference of the (log-transformed) estimates for aspartate aminotransferase and alanine aminotransferase.

Statin identification and LDL adjustment

Statin identification and LDL adjustment is described in Supplementary Note.

Covariate correction

Covariate adjustment is described in Supplementary Note.

Definition of type 2 diabetes

We used the definition of type 2 diabetes from our previous paper, including removal of type 1 diabetes cases from both cases and controls[80]. We use the terminology from Eastwood et al. throughout this description[81]. Type 2 diabetes was assigned case status for “probable type 2 diabetes” and “possible type 2 diabetes” and control status for “type 2 diabetes unlikely”; in addition, individuals with “probable type 1 diabetes,” “possible type 1 diabetes,” or “ probable gestational diabetes” were excluded. Finally, for the “strict” type 2 diabetes definition, we removed controls with HbA1c ≥ 39 mmol/mol.

Genome-wide association analysis

We performed genome-wide association analyses using the following four datasets: 1) The directly genotyped variants on the array (for protein-truncating and protein-altering variants); 2) The imputed variants (version 3); 3) The imputed HLA alleles; and 4) The copy number variations (CNVs) and gene-level aggregated CNV burden[34]. All the p-values from the association analyses are from two-sided tests. Detailed description of the association analysis is described in Supplementary Note.

Meta-analysis

Using the GWAS summary statistics for four analyzed populations (White British, Non-British White, South Asian, and African; East Asian GWAS were excluded) in UK Biobank, we performed inverse-variance weighted meta-analysis using METAL (version 2011-03-25) and included heterogeneity of effects analysis.

For the summary statistics from the meta-analysis, we checked whether the A1 and A2 alleles match with the alternate and reference allele in GRCh37/hg19 reference genome (fasta file) using bedtools getfasta subcommand[82], and canonicalized our association summary statistics such that we always report the effect size with respect to the alternate allele in the reference genome.

Derivation of independent loci

Once we ran the GWAS, full summary statistics were clumped to r2 > 0.1 using the following clump command:

plink1.9 --bfile <1000G Phase 3 European plink file> --clump