Genetically adjusted PSA levels for prostate cancer screening

Abstract

Prostate-specific antigen (PSA) screening for prostate cancer remains controversial because it increases overdiagnosis and overtreatment of clinically insignificant tumors. Accounting for genetic determinants of constitutive, non-cancer-related PSA variation has potential to improve screening utility. In this study, we discovered 128 genome-wide significant associations (P < 5 × 10⁻⁸) in a multi-ancestry meta-analysis of 95,768 men and developed a PSA polygenic score (PGS_PSA) that explains 9.61% of constitutive PSA variation. We found that, in men of European ancestry, using PGS-adjusted PSA would avoid up to 31% of negative prostate biopsies but also result in 12% fewer biopsies in patients with prostate cancer, mostly with Gleason score <7 tumors. Genetically adjusted PSA was more predictive of aggressive prostate cancer (odds ratio (OR) = 3.44, P = 6.2 × 10⁻¹⁴, area under the curve (AUC) = 0.755) than unadjusted PSA (OR = 3.31, P = 1.1 × 10⁻¹², AUC = 0.738) in 106 cases and 23,667 controls. Compared to a prostate cancer PGS alone (AUC = 0.712), including genetically adjusted PSA improved detection of aggressive disease (AUC = 0.786, P = 7.2 × 10⁻⁴). Our findings highlight the potential utility of incorporating PGS for personalized biomarkers in prostate cancer screening.

Main

Prostate-specific antigen (PSA) is an enzyme produced by the prostate gland that degrades gel-forming seminal proteins to release motile sperm and is encoded by the KLK3 (kallikrein 3) gene^1,2,3. As prostate epithelial tissue becomes disrupted by a tumor, greater PSA concentrations are released into circulation^2,3. PSA levels can also rise due to prostatic inflammation, infection, benign prostatic hyperplasia, older age and increased prostate volume^3,4,5. Increased body mass index is associated with lower PSA levels, but the underlying mechanisms remain unclear^6,7. Low PSA levels, thus, do not rule out prostate cancer, and PSA elevation is not sufficient for a conclusive diagnosis⁸. Although PSA testing reduces deaths from prostate cancer⁹, between 20% and 60% of cancers detected using PSA testing are estimated to be overdiagnoses^10,11,12. In addition, the long-term risk of lethal prostate cancer remains low, especially in men with PSA below the age-specific median^13,14. As a result, clinical guidelines in the United States and globally advise against population-level PSA screening and promote a shared decision-making model^15,16.

One avenue for refining PSA screening is to account for variability in PSA due to genetic factors. PSA is highly heritable, with 40 independent loci identified in the largest previous genome-wide association study (GWAS)^17,18. The goal of genetically correcting PSA levels is to increase the relative variation in PSA attributable to prostate cancer, thereby improving their predictive value for disease detection. The first study to genetically correct PSA using just four variants reclassified 3% of participants to warranting biopsy and 3% to avoiding biopsy¹⁹. Incorporating additional genetic predictors has the potential to personalize PSA testing, reduce overdiagnosis-related morbidity and improve detection of lethal disease. To maximize the utility of this approach, it is critical to distinguish genetic variants that influence constitutive PSA levels from those affecting prostate tumor development. PSA and prostate cancer share many genetic loci^{17,19,20,21,22}, but the extent to which this overlap reflects screening bias remains unclear, as GWASs of prostate cancer may capture signals for disease susceptibility and incidental detection due to benign PSA elevation.

Our study explores the genetic architecture of PSA levels in men without prostate cancer, with a view toward assessing whether genetic adjustment of PSA improves clinical decision-making related to prostate cancer diagnosis. It also provides a novel framework for the clinical translation of polygenic scores (PGSs) for non-causal cancer biomarkers.

Results

The study design of the Precision PSA study is illustrated in Fig. 1. Using data from five studies (Methods), we conducted genome-wide analyses of PSA levels ≤10 ng ml⁻¹ in cis-gender men never diagnosed with prostate cancer. GWAS results were meta-analyzed within ancestry groups and then combined across populations for a total sample size of 95,768 individuals.

Fig. 1: Overview of the Precision PSA study design.

Genome-wide association analyses were conducted in men without prostate cancer and meta-analyzed within each population: European ancestry (EUR), African ancestry (AFR), East Asian ancestry (EAS) and Hispanic/Latino ancestry (HIS/LAT). Ancestry-stratified results were used to develop a genome-wide PGS_PSA comprised of approximately 1.1 million variants and were also combined into a multi-ancestry meta-analysis of 95,768 men. PGS_PSA was validated in the PCPT and the SELECT and was used to compute PSA^G values. We examined how using PSA^G values affects eligibility for prostate biopsy and evaluated associations with incident prostate cancer.

Genetic architecture of PSA variation

The heritability (h²) of PSA levels was investigated using several methods to assess sensitivity to underlying modeling assumptions (Methods). Across 26,491 men of European ancestry in the UK Biobank (UKB) with linked clinical records, the median PSA value was 2.35 ng ml⁻¹ (Supplementary Fig. 1). Using individual-level data for variants with minor allele frequency (MAF) ≥ 0.01 and imputation INFO > 0.80, PSA heritability was h² = 0.41 (95% confidence interval (CI): 0.36–0.46) based on GCTA²³ and h² = 0.30 (95% CI: 0.26–0.33) based on LDAK²⁴ (Supplementary Table 1 and Extended Data Fig. 1). Applying LDAK to GWAS summary statistics generated from the same individuals produced similar estimates (h² = 0.35, 95% CI: 0.28–0.43), whereas other methods^25,26 were biased downward. In the European ancestry GWAS meta-analysis (n_EUR = 85,824), LDAK estimated h² = 0.30 (95% CI: 0.29–0.31). Sample sizes for other ancestries were too small for reliable heritability estimates.

The multi-ancestry meta-analysis of 95,768 men from five studies identified 128 independent index variants (P < 5.0 × 10⁻⁸, linkage disequilibrium (LD) r² < 0.01 within ±10-Mb windows) across 90 chromosomal cytoband regions (Fig. 2). The strongest associations were in known PSA loci^17,19,21,22, such as KLK3 (rs17632542, P = 3.2 × 10⁻⁶³⁸), 10q26.12 (rs10886902, P = 8.2 × 10⁻¹¹⁸), MSMB (rs10993994, P = 7.3 × 10^–87), NKX3-1 (rs1160267, P = 6.3 × 10⁻⁸³), CLPTM1L (rs401681, P = 7.0 × 10⁻⁵⁴) and HNF1B (rs10908278, P = 2.1 × 10⁻⁴⁶). Eighty-two index variants were independent of previously detected associations in the Genetic Epidemiology Research on Adult Health and Aging (GERA) cohort¹⁷; they mapped to 56 cytobands where PSA signals have not previously been reported. Associations initially detected in the UKB (Extended Data Fig. 1b) strengthened in the meta-analysis: TEX11 in Xq13.1 (rs62608084, P = 1.7 × 10⁻²⁴); THADA in 2p21 (rs11899863, P = 1.7 × 10⁻¹³); OTX1 in 2p15 (rs58235267, P = 4.9 × 10⁻¹³); SALL3 in 18q23 (rs71279357, P = 1.8 × 10⁻¹²); and ST6GAL1 in 3q27.3 (rs12629450, P = 2.6 × 10⁻¹⁰). Additional novel findings included CDK5RAP1 (rs291671, P = 1.2 × 10⁻¹⁸), LDAH (rs10193919, P = 1.5 × 10⁻¹⁵), ABCC4 (rs61965887, P = 3.7 × 10⁻¹⁴), INKA2 (rs2076591, P = 2.6 × 10⁻¹³), SUDS3 (rs1045542, P = 1.2 × 10⁻¹³), FAF1 (rs12569177, P = 3.2 × 10⁻¹³), JARID2 (rs926309, P = 1.6 × 10⁻¹²), GPC3 (rs4829762, P = 5.9 × 10⁻¹²), EDA (rs2520386, P = 4.2 × 10⁻¹¹) and ODF3 (rs7103852, P = 1.2 × 10⁻⁹) (Supplementary Tables 2 and 3).

Fig. 2: Multi-ancestry GWAS of PSA levels.

a, Manhattan plot depicting the results of the GWAS meta-analysis of PSA levels in 95,768 men without prostate cancer. The genome-wide significance threshold of P < 5 × 10⁻⁸ is indicated by the dotted black line. Index variants within known PSA-associated loci are annotated with the corresponding cytoband. Novel findings are highlighted in yellow. b, Circular dendogram shows the nearest gene(s) for novel PSA-associated variants. Genome-wide significant (P < 5 × 10⁻⁸) index variants were selected using LD-based clumping (LD r² < 0.01 within ±10-Mb windows). All GWAS P values are two-sided and derived from a fixed-effects inverse-variance-weighted meta-analysis using METAL.

Of the 128 index variants, 96 reached genome-wide significance in the European ancestry meta-analysis, as did three in the East Asian ancestry meta-analysis (KLK3: rs2735837 and rs374546878; MSMB: rs10993994; n_EAS = 3,337), two in the Hispanic/Latino meta-analysis (KLK3: rs17632542 and rs2735837; n_HIS/LAT = 3,098) and one in the African ancestry meta-analysis (FGFR2: rs10749415; n_AFR = 3,509) (Supplementary Table 4). Effect sizes from the European ancestry GWAS were modestly correlated with estimates from other ancestries (Spearman’s ρ_HIS/LAT = 0.48, P = 1.1 × 10⁻⁸; ρ_AFR = 0.27, P = 2.0 × 10⁻³; ρ_EAS = 0.16, P = 0.068) (Supplementary Fig. 2). However, cross-population comparisons of correlations should be interpreted with caution as they are confounded by higher sampling error in groups with smaller sample sizes.

There was heterogeneity (Cochran’s Q P_Q < 0.05) across ancestry-specific fixed-effects meta-analyses for 12 of 128 index variants, four of which had effects in different directions: rs58235267 (OTX1), rs1054713 (KLK1), rs10250340 (EIF4HP1) and rs7020681 (SLC35D2) (Supplementary Table 5). An alternative meta-analysis approach, MR-MEGA²⁷, which partitions effect size heterogeneity into components correlated with ancestry and residual variation, identified one additional signal in 5q15 (rs291812, P_MR-MEGA = 1.0 × 10⁻⁸) that was driven by the East Asian ancestry results (P_EAS = 1.2 × 10⁻⁶) (Supplementary Table 6).

Predicted functional consequences of the 128 index variants were explored using CADD²⁸. Scores >13 (corresponding to the 5% most deleterious substitutions genome-wide) were observed for 16 of the 128 index variants detected in the original fixed effects meta-analysis, including ten new signals: rs10193919 (LDAH); rs7732515 in 5q14.3; rs11899863 (THADA); rs58235267 (OTX1); rs926309 (JARID2); rs4829762 (GPC3) and rs13268, a missense variant in FBLN1; rs78378222 in TP53 and rs3760230 in SMG6; and rs712329 in SLC25A21 (Supplementary Table 7). Sixty-one variants had significant (false discovery rate (FDR) < 0.05) effects on gene expression, including 15 prostate tissue expression quantitative trait loci (eQTLs) for 17 eGenes, 55 blood eQTLs for 185 eGenes and nine eQTLs with effects in both tissues. Notable eGenes included RUVBL1, a chromatin-remodeling factor that modulates pro-inflammatory NF-κB signaling and transcription of Myc and β-catenin²⁹; ODF3, which maintains elastic structures in the sperm tail³⁰; and LDAH, which promotes cholesterol mobilization in macrophages³¹. Several PSA-associated variants were eQTLs for genes involved in immune response (IFITM2, IFITM3 and HS1BP3).

Impact of PSA-related selection bias on prostate cancer GWAS

Because prostate cancer detection often hinges on PSA elevation, genetic factors resulting in higher constitutive PSA levels may appear to increase prostate cancer risk because of more frequent screening. Of the 128 lead PSA variants, 52 (41%) were associated with prostate cancer at the Bonferroni-corrected threshold (P < 0.05/128) in the PRACTICAL consortium’s European ancestry GWAS³² (Supplementary Table 8). Using the method by Dudbridge et al.³³, we investigated whether index event bias could partly explain these shared signals^33,34 (Methods, Fig. 3 and Supplementary Table 9). Applying the estimated bias correction factor (b = 1.144) decreased the number of variants associated with prostate cancer from 52 to 34 (Extended Data Fig. 2). When we corrected 209 European ancestry prostate cancer risk variants (P < 5.0 × 10⁻⁸, LD r² < 0.01) for screening bias, 93 (45%) remained genome-wide significant. Notably, rs76765083 (KLK3) remained genome-wide significant but reversed direction. Sensitivity analyses using SlopeHunter³⁵ resulted in 150 (72%) variants with P < 5 × 10⁻⁸ (Supplementary Table 10).

Fig. 3: Influence of PSA-related index event bias on prostate cancer GWAS.

a, Conceptual diagram depicts how selection on PSA levels induces an association between genetic variant G and U, a composite confounder that captures polygenic and non-genetic factors. This selection induces an association with prostate cancer (PrCa) via path G–U→PrCa, in addition to the direct G→PrCa effects. Bi-directional dotted lines show that PSA is not only a disease biomarker but also influences screening behavior and the likelihood of prostate cancer detection. b,c, The impact of bias correction is shown for 209 prostate cancer risk variants. Independent risk variants were selected from the PRACTICAL GWAS meta-analysis (85,554 cases and 91,972 controls of European ancestry) by Conti et al.³² using LD clumping (LD r² < 0.01, P < 5 × 10⁻⁸). For each variant, associations with PSA (β_PSA) are based on an inverse-variance-weighted fixed-effects meta-analysis in men of European ancestry (n = 85,824). b, GWAS effect sizes for prostate cancer (β_PrCa) are aligned to the risk-increasing allele. Bias-adjusted effect sizes (β_adj) are denoted by triangles. c, Two-sided GWAS P values for prostate cancer (P_PrCa) were derived from an inverse-variance-weighted fixed-effects meta-analysis. Two-sided bias-adjusted P values (P_adj), denoted by triangles, were calculated from a chi-squared test statistic based on β_adj and corresponding standard errors. Genome-wide significance threshold (P < 5 × 10⁻⁸) is indicated by the horizontal dotted line.

Development and validation of PGS_PSA

We considered two approaches for constructing a PGS for PSA: clumping genome-wide significant associations from the multi-ancestry meta-analysis (PGS₁₂₈) and a genome-wide score generated using the Bayesian PRS-CSx algorithm (PGS_CSx) (ref.³⁶) (Methods). Each score was validated in the Prostate Cancer Prevention Trial (PCPT) and the Selenium and Vitamin E Cancer Prevention Trial (SELECT), which were excluded from the discovery GWAS. Most of the men in both cohorts were of European ancestry, although SELECT offered larger sample sizes for other ancestry groups (Extended Data Fig. 3). PGS_CSx was ultimately selected, as it was more predictive of baseline PSA than PGS₁₂₈ in multi-ancestry analyses and most ancestry subgroups (Supplementary Table 11).

In the PCPT, PGS_CSx accounted for 8.13% of variation in baseline PSA levels (β per s.d. increase = 0.186, P = 3.3 × 10⁻¹¹²) in the pooled multi-ancestry sample of 5,883 men (Fig. 4a–c and Supplementary Table 11). PGS_CSx was associated with PSA across age groups, although effects attenuated in participants aged ≥70 years (Extended Data Fig. 4). PGS_CSx was validated in 5,725 participants of European ancestry (EUR ≥ 0.80) (PGS_CSx: β = 0.194, P = 1.7 × 10⁻¹¹⁵), but neither PGS₁₂₈ nor PGS_CSx reached nominal significance in the admixed European and African ancestry (0.20 < AFR/EUR < 0.80, n = 103) or East Asian ancestry (EAS ≥ 0.80, n = 55) populations.

Fig. 4: Validation of the PGS_PSA in two cancer prevention trials.

a–f, Performance of PGS_PSA was evaluated in the PCPT and the SELECT. a,b, Violin plots show the distribution of baseline log(PSA) within quantiles of PGS_PSA, comprising 1,058,173 and 1,071,278 variants in the PCPT (a) and the SELECT (b). Box plots extend from the 25th to the 75th percentiles, with a trend line connecting the median value within each age stratum. Two-sided P values were derived from linear regression models for the effect of a quantile increase in PGS_PSA on log(PSA). c,d, Crossbar plots show the effect estimates (β) and corresponding 95% CIs per s.d. increase in the standardized PGS_PSA on baseline log(PSA) in the PCPT (c) and the SELECT (d). Ancestry-stratified and pooled multi-ancestry estimates are presented. Two-sided P values based on linear regression models are annotated. e,f, Comparison of distributions for PSA and PSA^G, with the horizontal line at 4 ng ml⁻¹, a commonly used threshold for further diagnostic testing. Box plots show the median value, with lower and upper hinges corresponding to the 25th and 75th percentiles or first and third quartiles. Whiskers extend as a multiple of the interquartile range (IQR × 1.5). Outlying values beyond the end of the whiskers are plotted individually.

In the SELECT, PGS_CSx was associated with baseline PSA levels in the pooled sample of 25,917 men (β = 0.258, P = 1.3 × 10⁻⁶¹⁹) and among men of European ancestry (n = 22,253, β_PGS = 0.283, P = 5.5 × 10⁻⁶¹⁰), accounting for 9.61% to 10.94% of variation, respectively (Fig. 4b–d and Supplementary Table 11). PGS_CSx also validated in the East Asian (n = 257, β = 0.258, P = 5.9 × 10⁻⁷) and admixed EAS/EUR (n = 321, β = 0.315, P = 5.2 × 10⁻¹²) ancestry groups. In men with admixed AFR/EUR ancestry (n = 1,763), PGS_CSx explained 4.22% of PSA variation (β = 0.157, P = 4.8 × 10⁻¹⁹). PGS₁₂₈ was more predictive than PGS_CSx (β = 0.163, P = 8.2 × 10⁻¹¹ versus β = 0.098, P = 8.0 × 10⁻⁶) in men of African ancestry (AFR ≥ 0.80, n = 1,173) and the pooled AFR and admixed (0.20 <EUR/AFR < 0.80) group (n = 2,936).

We also examined associations with temporal trends in PSA: velocity, calculated using log(PSA) values at two timepoints, and doubling time in months (Methods and Supplementary Table 12). In men with a PSA increase (SELECT pooled sample: n = 14,908), PGS_CSx was associated with less rapid velocity (PGS_CSx: β = −4.06 × 10⁻⁴, P = 3.7 × 10⁻⁵) and longer doubling time (PGS_CSx: β = 10.41, P = 1.9 × 10⁻⁸). In men with a PSA decrease between the first and last timepoint (SELECT pooled sample: n = 6,970), PGS_CSx was only suggestively associated with slowing PSA decline (β = 5.02 × 10⁻⁴, P = 0.068). The same pattern was observed in the PCPT, with higher PGS_CSx values conferring less rapid changes in PSA.

PGS_CSx, referred to as PGS_PSA from here onward, was used to genetically adjust baseline or earliest pre-randomization PSA values (PSA^G) for each individual, relative to the population mean (Methods and equations 1 and 2). PSA^G and unadjusted PSA were strongly correlated in the PCPT (Pearson’s r = 0.841, 0.833–0.848) and the SELECT (r = 0.854, 0.851–0.857). The number of participants with PSA^G > 4 ng ml⁻¹, a commonly used threshold for diagnostic testing, increased from 0 to 24 in the PCPT and from 5 to 413 in the SELECT (Fig. 4e,f), reflecting the preferential trial selection of men with low PSA^8,37.

Impact of PSA-related bias on PGS associations

In men of European ancestry in the UKB excluded from the PSA GWAS, there was a strong positive relationship between the 269-variant prostate cancer PGS (PGS₂₆₉)³² and PGS_PSA in cases (n = 11,568, β = 0.190, P = 2.3 × 10⁻⁹⁶) and controls (n = 152,884, β = 0.236, P < 10⁻⁷⁰⁰) (Extended Data Fig. 5 and Supplementary Table 13). Re-fitting PGS₂₆₉ using weights corrected for index event bias (PGS₂₆₉^adj) substantially attenuated associations in cases (β_adj = 0.029, P = 2.7 × 10⁻³) and controls (β_adj = 0.052, P = 2.2 × 10⁻⁸⁹).

To further characterize the impact of this bias, we examined PGS₂₆₉ associations with prostate cancer status in 3,673 cases and 2,363 biopsy-confirmed, European ancestry controls from GERA. PGS₂₆₉^adj had a larger magnitude of association with prostate cancer (OR for top decile = 3.63, 95% CI: 3.01–4.37) than PGS₂₆₉ (odds ratio (OR) = 2.71, 95% CI: 2.28–3.21) and higher area under the curve (AUC: 0.685 versus 0.677, P = 3.91 × 10⁻³) (Supplementary Table 14). The impact of bias correction was most pronounced for Gleason ≥7 tumors (PGS₂₆₉^adj AUC = 0.692 versus PGS₂₆₉ AUC = 0.678, P = 1.91 × 10⁻³), although these AUC estimates are inflated due to overlap with the GWAS used to develop PGS₂₆₉ (ref. ³²). In case-only analyses, PGS_PSA and PGS₂₆₉ were inversely associated with Gleason score, illustrating how screening bias decreases the likelihood of identifying high-grade disease (Supplementary Table 15). Compared to Gleason ≤6 tumors, an s.d. increase in PGS_PSA was inversely associated with Gleason 7 disease (OR = 0.79, 95% CI: 0.76–0.83) and Gleason ≥8 disease (OR = 0.71, 95% CI: 0.64–0.81). Patients in the top decile of PGS₂₆₉ were approximately 30% less likely to have Gleason ≥8 tumors (OR = 0.72, 95% CI: 0.54–0.96) than Gleason ≤6 tumors, but this association was attenuated after bias correction (PGS₂₆₉^adj: OR = 0.94, 95% CI: 0.75–1.17).

Impact of genetic adjustment of PSA on biopsy eligibility

Among GERA participants who underwent prostate biopsy, we examined how adjustment using PGS_PSA reclassified individuals for biopsy recommendation at age-specific thresholds used by Kaiser Permanente: 40–49 years old = 2.5 ng ml⁻¹; 50–59 years old = 3.5 ng ml⁻¹; 60–69 years old = 4.5 ng ml⁻¹; and 70–79 years old = 6.5 ng ml⁻¹ (Methods). For men of European ancestry, mean PSA levels in men with a negative biopsy (n = 2,363, 7.2 ng ml⁻¹) were higher than in men without prostate cancer who did not have a biopsy (n = 24,811, 1.5 ng ml⁻¹) (Supplementary Table 16). Relative to all controls, where standardized $\overline {PGS} _{PSA}$ = 0, biopsied men were enriched for PSA-increasing alleles (cases: $\overline {PGS} _{PSA}$ = 0.278; controls: $\overline {PGS} _{PSA}$ = 0.934). After genetic adjustment, 31.7% of biopsy-negative men were reclassified below the PSA level for recommending biopsy, and 2.5% became biopsy eligible, resulting in a net reclassification of 29.3% (27.5% to 31.21%) (Fig. 5a). Among 3,673 cases, PSA^G values below the biopsy referral threshold were more prevalent than upward adjustment, resulting in a net reclassification of −8.6% (−9.48% to −7.67%) (Fig. 5a). Of the patients who became ineligible, most had Gleason <7 tumors (n = 300, 72%; Supplementary Table 16). In men of African ancestry, there were few changes in biopsy eligibility among patients (n = 392), with 3.1% reclassified upward and 4.6% downward (Fig. 5b and Supplementary Table 16). Of 108 biopsy-negative controls, 75 (69.4%) were reclassified below the referral threshold based on PSA^G, reflecting high enrichment for predisposition to PSA elevation ($\overline {PGS} _{PSA}$ = 1.710). The overall net reclassification was positive, suggesting that PSA^G has some clinical utility in both populations.

Fig. 5: Genetically adjusted PSA influences biopsy eligibility.

a–b, Flow diagrams illustrate changes in PSA values after genetic adjustment for participants in the GERA cohort and subsequent reclassification at PSA thresholds used to recommend prostate biopsy. Genetic adjustment was applied to the last pre-biopsy PSA value to obtain PSA^G. Analyses were performed separately in men of European (a) and African (b) ancestry. Size of the nodes and flows are proportional to the number of individuals in each category. Patients with prostate cancer (cases) were stratified by Gleason score categories, where Gleason <7 represents potentially indolent disease. Gleason score is not applicable to men with a negative prostate biopsy (controls).

PSA genetic adjustment improves prostate cancer detection

The utility of PSA^G, alone and in combination with PGS₂₆₉, was first assessed in the PCPT, where end-of-study biopsies were performed in all participants, effectively eliminating potential misclassification of prostate cancer status. Among 335 cases and 5,548 controls, PGS_PSA was not associated with prostate cancer incidence (pooled: OR per s.d. = 1.01, P = 0.83), confirming that it captures genetic determinants of non-cancer PSA variation. The magnitude of association for genetically adjusted baseline PSA^G with prostate cancer (OR per unit increase in log(PSA ng ml⁻¹) = 1.90, 95% CI: 1.56–2.31) was slightly larger than for PSA (OR = 1.88, 95% CI: 1.55–2.29) in the European ancestry group (Supplementary Table 17). The magnitude of association with prostate cancer was larger for PGS₂₆₉^adj (pooled and European: OR per s.d. = 1.57, 95% CI: 1.40–1.76) than for PGS₂₆₉ without bias correction (pooled: OR = 1.52, 95% CI: 1.36–1.70; European: OR = 1.53, 95% CI: 1.36–1.72) (Supplementary Table 17). The model with PGS₂₆₉^adj and PSA^G achieved the best classification in the pooled (AUC = 0.686) and European ancestry (AUC = 0.688) populations and outperformed PGS₂₆₉^adj alone (pooled: AUC = 0.656, P_AUC = 7.5 × 10⁻⁴; European: AUC = 0.658, P_AUC = 1.4 × 10⁻³).

The benefit of genetically adjusting PSA was most evident for detection of aggressive prostate cancer, defined as Gleason ≥7, PSA ≥ 10 ng ml⁻¹, T3–T4 stage and/or distant or nodal metastases. In the PCPT, PSA^G conferred an approximately threefold risk increase (pooled: OR = 2.87, 95% CI: 1.98–4.65, AUC = 0.706; European: OR = 2.99, 95% CI: 1.95–4.59, AUC = 0.711) compared to PGS₂₆₉^adj (pooled: OR = 1.55, 95% CI: 1.23–1.95, AUC = 0.651; European: OR = 1.55, 95% CI: 1.22–1.96, AUC = 0.657) (Fig. 6a and Supplementary Table 18). The model with PSA^G and PGS₂₆₉^adj achieved AUC = 0.726 (European: AUC = 0.734) for aggressive tumors but had lower discrimination for non-aggressive disease (pooled and European: AUC = 0.681) (Supplementary Table 19). Among patients with prostate cancer, PSA^G (pooled: OR = 2.06, 95% CI: 1.23–3.45) and baseline PSA (pooled: OR = 1.81, 85% CI: 1.12–3.10) were associated with higher likelihood of aggressive compared to non-aggressive tumors, whereas PGS₂₆₉ (pooled: OR = 0.91, P = 0.54) and PGS₂₆₉^adj (OR = 0.97, P = 0.85) were not (Supplementary Table 20).

Fig. 6: Genetic associations with aggressive prostate cancer.

a–b, Comparison of models for aggressive disease, defined as Gleason score ≥7, PSA ≥ 10 ng ml⁻¹, T3–T4 stage and/or distant or nodal metastases in the PCPT (a) and the SELECT (b). The pooled study population includes all ancestry groups. Logistic regression models were adjusted for baseline age, randomization arm, the top ten population-specific genetic ancestry principal components and proportions of African and East Asian genetic ancestry. ORs and 95% CIs were estimated per 1-unit increase in log(PSA ng ml⁻¹) and log(PSA^G ng ml⁻¹) and per s.d. increase in the prostate cancer genetic risk score (PGS₂₆₉) from Conti et al.³², which was standardized to achieve s.d. equal to 1. All P values are two-sided. AUC is based on the full model with all covariates.

In the SELECT, associations with risk of prostate cancer overall (Supplementary Table 21), aggressive disease (Fig. 6b and Supplementary Table 22) and non-aggressive disease (Supplementary Table 23) in the pooled and European ancestry analyses were similar to the PCPT. In men of East Asian ancestry, associations for PSA^G (OR = 2.15, 95% CI: 0.82–5.62) were attenuated compared to PSA (OR = 2.60, 95% CI: 1.03–6.54). This was also observed in men of African ancestry, although the effect size for PSA^G derived using PGS₁₂₈ (OR = 3.37, 95% CI: 2.38–4.78) was larger than for PSA^G based on PGS_CSx (OR = 2.68, 95% CI: 1.94–3.69), consistent with the larger proportion of variation in PSA explained by PGS₁₂₈ than PGS_CSx in this population. Models for prostate cancer including PSA^G were calibrated in the pooled and European ancestry individuals, whereas, in the African ancestry subgroup, PSA^G inaccurately estimated risk in upper deciles (Supplementary Figs. 3–6).

The largest improvement in discrimination from PSA^G (OR = 3.81, 95% CI: 2.62–5.54, AUC = 0.777) relative to PSA (OR = 3.40, 95% CI: 2.34–4.93, AUC = 0.742, P_AUC = 0.026) and to PGS₂₆₉ (OR = 1.76, 95% CI: 1.41–2.21, AUC = 0.726, P_AUC = 0.057) was for aggressive tumors in men of European ancestry (106 cases, 23,667 controls). In the pooled African ancestry population (18 cases, 2,733 controls), PSA^G based on PGS₁₂₈ (OR = 2.96, 95% CI: 1.43–6.12), but not PGS_CSx (OR = 2.48, 95% CI: 1.24–4.97), was more predictive than unadjusted PSA (OR = 2.82, 95% CI: 1.33–5.99) (Supplementary Table 22). The best model for aggressive disease included PSA^G and PGS₂₆₉^adj for pooled (AUC = 0.788, 95% CI: 0.744–0.831) and European ancestry (AUC = 0.804, 95% CI: 0.757–0.851) populations, but, for African ancestry individuals, unadjusted PSA and PGS₂₆₉ without bias correction achieved the highest AUC of 0.828 (95% CI: 0.739–0.916). PSA^G was better calibrated than PSA in pooled and European ancestry groups but not in African ancestry participants (Supplementary Figs. 7 and 8).

Discussion

Serum PSA is the most widely used biomarker for prostate cancer detection, although concerns with specificity and, to a lesser degree, sensitivity have limited adoption of PSA testing for population-level screening. Leveraging PGS to personalize diagnostic biomarkers, such as PSA, provides a new avenue for translating GWAS discoveries into clinical practice. This concept, termed ‘de-Mendelization’, is essentially Mendelian randomization in reverse—subtracting the genetically predicted component of trait variance instead of using it to estimate causal effects. De-Mendelization of non-causal predictive biomarkers can maximize disease-related signal and improve disease detection^38,39. Although previous work on PSA genetics¹⁹ and other biomarkers^38,40 has alluded to the potential of genetic adjustment to produce clinically meaningful shifts in the PSA distribution, the value of this approach for reducing overdiagnosis and detecting aggressive disease has not been previously shown.

Risk-stratified, personalized screening for prostate cancer will require parallel efforts to elucidate the genetic architecture of prostate cancer susceptibility and PSA variation in individuals without disease. Our GWAS advances these efforts by discovering 82 novel PSA-associated variants. The strongest novel signals map to genes involved in reproductive processes, potentially reflecting non-cancer function of PSA in liquefying seminal fluid. TEX11 on Xq13.1, for example, is preferentially expressed in male germ cells and early spermatocytes. TEX11 mutations cause meiotic arrest and azoospermia, and this gene regulates homologous chromosome synapsis and double-strand DNA break repair⁴¹. ODF3 encodes a component of sperm flagella fibers and has been linked to regulation of platelet count and volume⁴². Other novel loci contained genes involved in embryonic development, epigenetic regulation and chromatin organization, including DNMT3A, OTX1, CHD3, JARID2, HMGA1, HMGA2 and SUDS3. DNMT3A is a methyltransferase that regulates imprinting and X-chromosome inactivation and has been studied extensively in the context of height⁴³, clonal hematopoiesis and hematologic cancers⁴⁴. CHD3 is involved in chromatin remodeling during development and suppresses herpes simplex virus infection⁴⁵. Multiple PSA-associated variants were in genes related to infection and immunity, including HLA-A; ST6GAL1, involved in IgG N-glycosylation⁴⁶; KLRG1, which regulates natural killer (NK) cell function and IFN-γ production⁴⁷; and FUT2, which affects ABO precursor H antigen presentation and confers susceptibility to viral and bacterial infections⁴⁸.

Although our GWAS was restricted to men without prostate cancer, several cancer susceptibility genes were among the PSA-associated loci, including a pan-cancer risk variant in TP53 (rs78378222) (ref.⁴⁹) and signals in TP63, GPC3 and THADA. Although we cannot rule out undiagnosed prostate cancer in our participants, its prevalence is unlikely to be high enough to produce appreciable bias. Pervasive pleiotropy and omnigenic architecture⁵⁰ may explain the diverse functions of PSA loci implicated in inflammation, epigenetic regulation and growth factor signaling. Even established tumor suppressor genes, such as TP53, GPC3 and THADA, have pleiotropic effects on obesity via dysregulation of cell growth and metabolism^51,52,53. Furthermore, distinct p63 isoforms regulate epithelial and craniofacial development as well as apoptosis of male germ cells and spermatogenesis^54,55. Mutations in GPC3 cause Simson–Golabi–Behmel syndrome, which is characterized by visceral and skeletal abnormalities and excess risk of embryonic tumors⁵⁶.

Distinguishing variants that influence prostate cancer detection via PSA screening from genetic signals for prostate carcinogenesis has implications for deciphering biological mechanisms and developing risk prediction models. Prostate cancer detection depends on PSA testing, whereas PSA screening is influenced by genetic factors affecting constitutive PSA levels. The bias arising from this complex relationship may be substantial. Our findings suggest that bias-corrected effect sizes more accurately capture the contribution of GWAS-identified variants to prostate cancer risk, without conflating it with detection. Correction for PSA-related bias and subsequent improvement in PGS₂₆₉ performance for detecting aggressive disease is an extension of de-Mendelization. Adjusting risk allele weights may be a more effective strategy than filtering out variants based on associations with PSA. Generally, the improvements in PSA^G and PGS₂₆₉ are proportional to the extent of their de-noising of signals for PSA elevation unrelated to prostate cancer. The impact of bias correction was most pronounced in populations selected for high PSA, such as men who underwent prostate biopsy in GERA, but it was also observed in the PCPT and the SELECT, which enrolled men with low PSA.

Our investigation of index event bias has several limitations. The Dudbridge method assumes that direct genetic effects on PSA levels and prostate cancer susceptibility are uncorrelated, and violations of this assumption over-attribute shared genetic signals to selection bias³³. Although SlopeHunter relaxes this assumption³⁵, analyses of PGS₂₆₉ suggest that it under-corrects selection bias. SlopeHunter relies on clustering to distinguish PSA-specific from pleiotropic variants³⁵, with small or poorly separated clusters resulting in unstable bias estimates. Disentangling genetic associations between PSA and prostate cancer with greater certainty will require experiments such as CRISPR screens and massively parallel reporter assays.

Another limitation is that the reported magnitude of biopsy reclassification may be specific to GERA and Kaiser Permanente clinical guidelines and biased because GERA controls comprised 30% of the PSA discovery GWAS. Because it was unlikely for men with low PSA to be biopsied, and most patients with prostate cancer already had PSA values at or above the biopsy referral cutoff, there were limited opportunities to increase biopsy eligibility in this population. Despite these limitations, our findings indicate that genetically adjusted PSA may reduce overdiagnosis and overtreatment, albeit accompanied by some undesirable loss of sensitivity. Although reclassifying cases to not receive biopsy is concerning, most such reclassifications occurred among patients with non-aggressive disease, a group susceptible to overdiagnosis⁵⁷.

Our PGS-based approach updates the first application of PSA genetic correction by Gudmundsson et al.¹⁹ while retaining straightforward calculation of the genetic correction factor. Increasing the specificity of an established, clinically useful biomarker is efficient and would have low adoption barriers. However, analytic choices, such as selecting an optimal PGS algorithm and reference population for obtaining mean PGS_PSA, are not trivial. The choice of reference population affects the magnitude of correction and clinical decisions based on absolute PSA values. Furthermore, any new biomarker would require validation in real-world settings to identify populations who would benefit most and characterize barriers to implementation, such as physician familiarity with PGS and patient education about genetic testing. Genetically adjusted PSA should also be evaluated in conjunction with other procedures used for prostate cancer detection, such as targeted magnetic resonance imaging, and explored as a criterion for refining selection of participants into screening trials.

Our study highlights the importance and challenge of developing a PGS that adequately performs across the spectrum of ancestry. Compared to PGS₁₂₈, PGS_CSx did not improve performance in men of African ancestry. This may reflect the ‘meta’ estimation procedure, which does not require a separate dataset for hyperparameter tuning but is less accurate³⁶. GWAS efforts in larger and more diverse cohorts are underway and will expand the catalog of PSA-associated variants and increase their utility. Genetic adjustment using a PGS_PSA that does not explain a sufficiently high proportion of trait variation risks decreasing the accuracy of PSA screening.

Future research should assess whether genetically adjusted PSA levels improve prediction of prostate cancer mortality and investigate PSA-related biomarkers, such as the ratio of free to total PSA and pro-PSA (a precursor PSA isoform), which may have higher specificity for prostate cancer detection^58,59. Although PGS_PSA was associated with PSA doubling time and velocity, these metrics assess change between two timepoints and may not capture PSA trajectories that are meaningful for disease detection⁶⁰. Clinical guidelines for PSA kinetics are also lacking in the context of prostate cancer screening. Regardless, we think that genetic adjustment may improve the accuracy of any heritable PSA biomarker and may be a valuable addition to multi-omic biomarkers.

In summary, by detecting genetic variants associated with non-prostate cancer PSA variation, we developed a PGS_PSA that captures the contribution of common genetic variants to a man’s inherent PSA level. We showed that a straightforward calculation of genetically adjusted, personalized PSA levels using PGS_PSA provides clinically meaningful improvements in prostate cancer diagnostic characteristics. Moreover, genetic determinants of PSA provide an avenue for mitigating selection bias due to PSA screening in prostate cancer GWASs and improving disease prediction. These results illustrate a proof of concept for incorporating genetic factors into PSA screening for prostate cancer and expanding this approach to other diagnostic biomarkers.

Methods

Informed consent was obtained from all study participants. The UKB received ethics approval from the Research Ethics Committee (reference: 11/NW/0382) in accordance with the UKB Ethics and Governance Framework. The research was conducted with approved access to UKB data under application number 14105. We used previously published PSA GWAS results from the GERA cohort by Hoffmann et al.¹⁷. The original study was approved by the Kaiser Permanente Northern California institutional review board and the University of California, San Francisco Human Research Protection Program Committee on Human Research. The Prostate, Lung, Colorectal and Ovarian (PLCO) Cancer Screening Trial was approved by the institutional review board at each participating center and the National Cancer Institute. The informed consent document signed by PLCO study participants allows use of these data by investigators for discovery and hypothesis generation in the investigation of the genetic contributions to cancer and other adult diseases. Our study includes publicly posted genomic summary results from the PLCO Atlas⁶¹. No institutional review board review is required for PLCO summary data use. The Vanderbilt University Medical Center institutional review board approved the BioVU study. The Malmö Diet and Cancer Study (MDCS) was approved by the local ethics committee.

Study populations and phenotyping

Genome-wide association analyses of PSA levels were conducted using germline genetic data derived from DNA extracted from non-prostatic tissues (for example, blood and buccal swabs). Analyses were restricted to cis-gender men, defined as individuals of biological male sex and self-reported male gender identity who had never been diagnosed with prostate cancer. Men with a history of surgical resections of the prostate were excluded in studies for which this information was available. To reduce potential for reverse causation, analyses were limited to PSA values ≤10 ng ml⁻¹, which corresponds to low-risk prostate cancer based on the D’Amico prostate cancer risk classification system⁶², and PSA > 0.01 ng ml⁻¹, to ensure that individuals had a functional prostate not impacted by surgery or radiation.

The UKB is a population-based prospective cohort of over 500,000 individuals aged 40–69 years at enrollment in 2006–2010 with genetic and phenotypic data⁶³. Health-related outcomes were ascertained via individual record linkage to national cancer and mortality registries and hospital inpatient encounters. PSA values were abstracted from primary care records for a subset of participants with genetic data. Field code mappings used to identify PSA values included any serum PSA measure except for free PSA or ratio of free to total PSA (Supplementary Table 25).

The Kaiser Permanente GERA cohort used in this analysis was previously described in Hoffmann et al.¹⁷. In brief, prostate cancer status was ascertained from the Kaiser Permanente Northern California Cancer Registry, the Kaiser Permanente Southern California Cancer Registry or through review of clinical electronic health records. PSA levels were abstracted from Kaiser Permanente electronic health records from 1981 through 2015.

The PLCO Cancer Screening Trial is a completed randomized trial that enrolled approximately 155,000 participants between November 1993 and July 2001. The PLCO Cancer Screening Trial was designed to determine the effects of screening on cancer-related mortality and secondary endpoints in men and women aged 55–74 years⁶⁴. Men randomized to the screening arm of the trial underwent annual screening with PSA for 6 years and digital rectal exam (DRE) for 4 years⁶⁴. These analyses were limited to men with a baseline PSA measurement who were randomized to the screening arm of the trial (n = 29,524). Men taking finasteride at the time of PSA measurement were excluded from analysis.

The Vanderbilt University Medical Center BioVU resource is a synthetic derivative biobank linked to de-identified electronic health records⁶⁵. Analyses were based on PSA levels that were measured as part of routine clinical care.

The MDCS is a population-based prospective cohort study that recruited men and women aged between 44 years and 74 years of age who were living in Malmö, Sweden between 1991 and 1996 to investigate the impact of diet on cancer risk and mortality⁶⁶. These analyses included men from the MDCS who were not diagnosed with prostate cancer as of December 2014 and had available genotyping and baseline PSA measurements⁶⁶.

The PCPT is a completed phase 3 randomized, double-blind, placebo-controlled trial of finasteride for prostate cancer prevention that began in 1993 (ref. ⁸). The PCPT randomly assigned 18,880 men aged 55 years or older who had a normal DRE and PSA level ≤3 ng ml⁻¹ to either finasteride or placebo. For men with multiple pre-randomization PSA values, the earliest value was selected. Cases included all histologically confirmed prostate cancers detected during the 7-year treatment period and tumors that were detected by the end-of-study prostate biopsy. Our analyses included the subset of PCPT participants that was genotyped on the Illumina Infinium Global Screening Array 24 v2.0.

The SELECT is a completed phase 3 randomized, placebo-controlled trial of selenium (200 µg per day from l-selenomethionine) and/or vitamin E (400 IU per day of all rac-α-tocopheryl acetate) supplementation for prostate cancer prevention³⁷. Between 2001 and 2004, 34,888 eligible participants were randomized. The minimum enrollment age was 50 years for African American men and 55 years for all other men³⁷. Additional eligibility requirements included no prior prostate cancer diagnosis, ≤4 ng ml⁻¹ of PSA in serum and a DRE not suspicious for cancer. For men who had multiple pre-randomization PSA values, the earliest value was selected. Our analyses included a subset of SELECT participants genotyped on the Illumina Infinium Global Screening Array 24 v2.0.

Quality control and genome-wide association analyses

Standard genotyping and quality control (QC) procedures were implemented in each participating study. Before meta-analysis, we applied variant-level QC filters that included low imputation quality (INFO < 0.30), MAF < 0.005 and deviations from Hardy–Weinberg equilibrium (P_HWE < 1 × 10⁻⁵). Sample-level filtering was performed to remove samples with discordant genetic sex and self-reported gender and call rate < 0.97. One sample from each pair of first-degree relatives was also excluded. GWAS phenotypes and adjustment covariates are reported in Supplementary Table 26. Genome-wide association analyses performed linear regression of log(PSA) as the outcome, using age and genetic ancestry principal components (PCs) as the minimum set of covariates.

UKB

Genotyping and imputation for the UKB cohort were previously described⁶³. In brief, participants were genotyped on the UKB Affymetrix Axiom array (89%) or the UK BiLEVE array (11%) with imputation performed using the Haplotype Reference Consortium (HRC) and the merged UK10K and 1000 Genomes phase 3 reference panels. Genetic ancestry PCs were computed using fastPCA based on a set of 407,219 unrelated samples and 147,604 genetic markers⁶³. Association analyses in the UKB were restricted to individuals of European ancestry based on self-report (‘White’) and after excluding samples with either of the first two genetic ancestry PCs outside of 5 s.d. of the population mean, as previously descibed⁴⁹. We removed samples with discordant self-reported and genetic sex as well as one sample from each pair of first-degree relatives identified using KING⁶⁷. Using a subset of genotyped autosomal variants with MAF ≥ 0.01 and call rate ≥ 97%, we filtered samples with heterozygosity >5 s.d. from the mean. For participants with multiple PSA measurements, the median value PSA was used. Sensitivity analyses were conducted comparing this approach to a GWAS of individual-specific random effects derived from fitting a linear mixed model to repeated log(PSA) values.

GERA

Genotyping, imputation and QC of the GERA cohort were previously described^17,68,69. In brief, all men were genotyped for over 650,000 single-nucleotide polymorphisms (SNPs) on four race/ethnicity-specific Affymetrix Axiom arrays that were optimized for individuals who self-identified as non-Hispanic white, Latino, East Asian and African American, respectively^68,69. Genotype QC procedures and imputation for the original GERA cohort were performed on an array-wise basis, as previously described^17,70. Pre-phasing was done by SHAPEIT version 2.5 (ref. ⁷¹) and imputation with IMPUTE2 version 2.3.1 (ref. ⁷²) using the 1000 Genomes phase 3 release with 2,504 samples. The top ten genetic ancestry PCs from EIGENSTRAT version 4.2 were included in the linear model as ancestry covariates⁷³. Analyses were conducted according to self-identified race/ethnicity groups. Residuals were computed from linear mixed models that were fit to repeated log(PSA) measures. This approach was nearly identical to a long-term average, except that it used the median instead of the mean to handle any potential outlier PSA level values.

PLCO Atlas

Our study used GWAS summary statistics from the PLCO Atlas Project, a resource for multi-trait GWAS. Genotyping, QC and imputation procedures for this resource are described by Machiela et al.⁶¹. The Atlas Project combined genotyping data previously generated by high-density arrays for 25,831 participants (OncoArray, Omni2.5M and OmniExpress) with a new round of genotyping using the Illumina Global Screening Array (GSA). For participants genotyped on multiple genotyping arrays (n = 1,192), data from only one array were retained, with the following prioritization: GSA > OncoArray > Omni2.5M > OmniExpress. Extensive QC filtering was performed for subsequent imputation and association analyses. Iterative 80% and 95% sample-level and variant-level call rate filters were applied to remove poorly genotyped samples and variants. Samples with > 20% estimated contamination based on VerifyIDintensity⁷⁴ were also removed. Samples with discordant self-reported gender and genetically inferred sex were identified based on X-chromosome method-of-moments F coefficient from PLINK, using 0.5 as the threshold (F coefficients are close to 0.0 for males and 1.0 for females). Heterozygosity outliers were detected using absolute values from PLINK method-of-moments F coefficients > 0.2.

Genetic ancestry was determined using GRAF⁷⁵ on a set of 10,000 pre-selected fingerprinting variants. Participants were assigned to nine ancestral groups: ‘African’, ‘African American’, ‘East Asian’, ‘European’, ‘Hispanic1’, ‘Hispanic2’, ‘Other’, ‘Other Asian’ and ‘South Asian’. Hispanic1 included individuals of Dominican or Puerto Rican ancestry, whereas Hispanic2 included individuals of Mexican or Latin American ancestry. For parsimony, we merged ‘African’ and ‘African American’ into an ‘African American (Combined)’ and ‘East Asian’ and ‘Other Asian’ into an ‘East Asian (Combined)’. Imputation was performed using the TOPMed 5b reference panel, which is accessible via the TOPMed Imputation Server hosted on the Michigan Imputation Server. Before imputation, variants with MAF ≤ 0.01, missingness ≥ 0.05 and Hardy–Weinberg deviations (P_HWE ≤1 × 10⁻⁶) were removed. Genotyped data were aligned to reference datasets using a community-recommended script (HRC-1000G-check-bim.pl from https://www.well.ox.ac.uk/~wrayner/tools/) that was modified to support the TOPMed 5b reference panel using a pre-existing test imputation with 1000 Genomes subjects. Pre-phasing using phased reference data from TOPMed release 5b was conducted using Eagle 2.4 (ref. ⁷⁶). Imputation was conducted against the same reference panel using minimac4. GWAS was based on the first PSA value for each PLCO participant.

BioVU

Participants were identified using Vanderbilt University Medical Center’s BioVU resource, a DNA biobank comprising ~270,000 individuals and linked to a de-identified electronic health record⁶⁵. All participants (n = 8,074) were genotyped on Illumina’s Expanded Multi-Ethnic Genotyping Array (MEGA^EX) platform. Genetic ancestries were assigned by running principal component analysis using SNPRelate⁷⁷ on a set of pruned SNPs (Rsq < 0.5, MAF ≥ 0.1). Participants were classified as European ancestry if their first two PCs were within 4 s.d. of the median for the participants reporting ‘White’ as their race. Participants were classified as African ancestry if their first two PCs were within 4 s.d. of the median for participants reporting their race as ‘Black’. All QC procedures were performed using PLINK version 1.90. We removed one randomly selected sample out of each pair of related individuals (pi-hat ≥ 0.2) identified using identity-by-descent. We excluded participants with SNP missingness > 3% or heterozygosity >5 s.d. from the mean. Before imputation, data were pre-processed using the HRC-1000G-check-bim.pl (from http://www.well.ox.ac.uk/~wrayner/tools/) and pre-phased using Eagle version 2.4 (ref. ⁷⁶). Genetic data were imputed on the Michigan Imputation Server using 1000 Genomes phase 3 version 5 as the reference panel. For men with multiple PSA measurements, the median PSA was used.

MDCS

Data from multiple batches of genotyping of 4,069 MDCS participants using different Illumina Omni arrays were merged. For variants that appeared more than once under different names on the same Illumina array, those with the higher genotyping rate were retained. Indels, ambiguous palindromic (for example, A/T or C/G alleles) and multi-allelic variants were removed. Only SNPs that we could unambiguously map to the 1000 Genomes phase 1 dataset were kept. Individuals with > 10% missingness were removed. Next, SNPs with a missingness rate  > 10% or deviation from Hardy–Weinberg equilibrium (P_HWE < 0.001) were removed. At this stage, the PCs of ancestry were computed. Individuals for whom the inferred sex based on X-chromosome heterozygosity was not male, or for whom there were more than two genetic mismatches with 40 SNPs that we had previously genotyped in these samples with targeted genotyping⁶⁶, were excluded.

To assess genetic ancestry, MDCS data were combined with data from HapMap phase 3 for variants present in all genotyping batches. These SNPs were further filtered to have < 0.01% missingness and LD pruned (–indep-parwise 50 5 0.05). SMARTPCA in EIGENSOFT (https://github.com/chrchang/eigensoft) was run on the resulting 18,299 SNPs to generate the top ten genetic ancestry PCs. Analyses were restricted to individuals of European ancestry based on clustering with HapMap reference populations and exclusion of outliers with a z-score on PC1 and PC2 > 5. Imputation was performed using the TOPMed 5b reference panel, which is accessible via the TOPMed Imputation Server hosted on the Michigan Imputation Server. Before imputation, the input file was aligned to the build37 reference genome on the basis of chromosome, position and alleles. A total of 847,133 SNPs that passed pre-imputation QC were uploaded to the imputation server. From the resulting imputed files, analyses were restricted to individuals without a prostate cancer diagnosis by 31 December 2014, with individual missingness < 3% and a z-score < 5.0 for heterozygosity. Log(PSA) values were analyzed using robust linear regression with Tukey biweights. GWAS was performed using linear regression on the residuals extracted from the fitted models.

PCPT and SELECT

Participants from PCPT and SELECT were genotyped on the Illumina Infinium Global Screening Array 24 v2.0 and underwent the same QC and imputation procedures. Genotyping calling and QC were performed at the Center for Inherited Disease Research at Johns Hopkins. After removal of samples that failed to produce valid output during initial processing and clustering, the completion rate was 0.9951 and 0.9959 in PCPT and SELECT, respectively. A two-stage filter by completion rate threshold of 0.8 for samples and 0.8 for variants, followed by 0.95 for samples and 0.95 for variants, was performed. Samples with discordant self-reported gender and genetically inferred sex were identified based on X-chromosome method-of-moments F coefficient from PLINK, using 0.5 as the threshold (F coefficients are close to 0.0 for males and 1.0 for females). Identity-by-descent for all subject pairs was determined using PLINK, with close (first and second degree) relatives identified based on a threshold of 0.20. One randomly selected sample from each pair of relatives was retained.

Ancestry was estimated using a set of LD-pruned markers and running SNPWEIGHTS⁷⁸ with the reference panel provided containing the following populations: European, West African and East Asian, with a threshold of 0.8 used for imputed ancestry designation. Participants were assigned to a single ancestry group if the ancestry score was ≥0.80 for just one group. Participants were assigned to an admixed cluster if their ancestry score was > 0.20 and <0.80 for only one group (for example, ADMIXED_AFR where AFR = 0.75, EUR = 0.17, EAS = 8). Intermediate ancestry clusters included individuals with ancestry scores matching those criteria in multiple groups: 0.20 < AFR_EUR < 0.80 (for example, AFR = 0.65, EUR = 0.33) and 0.20 < EAS_EUR < 0.80 (for example, EUR = 0.55, EAS = 0.43). Autosomal heterozygosity was assessed using the method-of-moments F coefficient calculated within each ancestry cluster. Heterozygosity outliers were identified and excluded using a threshold of 0.10. Principal component analysis was performed with SMARTPCA in EIGENSOFT (https://github.com/chrchang/eigensoft) on a set of LD-pruned markers after splitting by ancestry cluster, to resolve more detailed population substructure. Genetic ancestry PCs were not computed for small clusters (n < 50) or individuals who failed other QC filters. For validation of PGS_PSA in PCPT and SELECT, we combined ADMIXED_AFR and AFR_EUR and treated this as a single group with admixed AFR and EUR ancestry proportions (AFR/EUR). ADMIXED_EAS and EAS_EUR were also combined into a single cluster with admixed EAS and EUR ancestry (EAS/EUR).

To prepare genotype data for imputation with the TOPMed 5b reference panel, variants with MAF < 0.001, call rate < 98% or evidence of deviation from Hardy–Weinberg equilibrium (P_HWE < 10⁻⁶) were removed. After these QC steps, a total of 474,046 variants remained for PCPT, and 491,015 variants were retained for SELECT. Before submitting the data to the TOPMed Imputation Server, files were pre-processed using the check-bim.pl script (http://www.well.ox.ac.uk/~wrayner/tools/). Next, chromosomal positions were lifted over from GRCh37/hg19 to GRCh38 and aligned against the TOPMed reference SNP list based on chromosome, position and alleles to ensure that reference and alternate alleles were correct in the resulting VCF files.

Heritability of PSA levels attributed to common variants

Heritability of PSA levels was estimated using individual-level data and GWAS summary statistics. UKB participants with available PSA and genetic data were analyzed using LDAK version 5.1 (ref. ²⁴) and GCTA version 1.93 (ref. ²³), following the approach previously implemented in the GERA cohort¹⁷. Genetic relationship matrices were filtered to ensure that no pairwise relationships with kinship estimates >0.05 remained. Heritability was estimated using common (MAF ≥ 0.01) LD-pruned (r² < 0.80) variants with imputation INFO > 0.80. We implemented the LDAK-Thin model using the recommended genetic relatedness matrix (GRM) settings (INFO > 0.95, LD r² < 0.98 within 100 kb) and the same parameters as GCTA for comparison (LD r² < 0.80, INFO > 0.80). For both methods, sensitivity analyses were conducted using more stringent GRM settings (kinship = 0.025, genotyped variants).

Summary statistics from GWAS results based on the same set of UKB participants (n = 26,491) and from a European ancestry GWAS meta-analysis (n = 85,824) were analyzed using LDAK, LD score regression (LDSR)²⁵ and an extension of LDSR using a high-definition likelihood (HDL) approach²⁶. For LDSR, we used the default panel comprising variants available in HapMap3 with weights computed in 1000 Genomes version 3 EUR individuals and in-house LD scores computed in UKB European ancestry participants⁴⁹. The baseline linkage disequilibrium (BLD)-LDAK model was fit using pre-computed tagging files calculated in UKB GBR (white British) individuals for HapMap3 variants from the LDSR default panel. HDL analyses were conducted using the UKB-derived panel restricted to high-quality imputed HapMap3 variants²⁶. All GWAS summary statistics had sufficient overlap with the reference panels, not exceeding the 1% missingness threshold for HDL and the 5% missingness threshold for LDAK and LDSR.

Genome-wide meta-analysis

Each ancestral population was analyzed separately, and GWAS summary statistics were combined via meta-analysis (Fig. 1). We first used METAL⁷⁹ to conduct an inverse-variance-weighted fixed-effects meta-analysis in each ancestry group and then meta-analyzed the ancestry-stratified results. Multi-ancestry meta-analysis results were processed using clumping to identify independent association signals by grouping variants based on LD within specific windows. Clumps were formed around index variants with the lowest genome-wide significant (P < 5 × 10⁻⁸) meta-analysis P value. All other variants with LD r² > 0.01 within a ±10-Mb window were considered non-independent and assigned to that lead variant. Since over 90% of the meta-analysis consisted of individuals of European ancestry, clumping was performed using 1000 Genomes phase 3 EUR and UKB reference panels, which yielded concordant results. We confirmed that LD among the resulting lead variants did not exceed r² = 0.05 using a merged 1000 Genomes ALL reference panel.

We first examined heterogeneity in the multi-ancestry fixed-effects meta-analysis results using Cochran’s Q statistic. To assess heterogeneity specifically due to ancestry, we applied MR-MEGA²⁷, a meta-regression approach for aggregating GWAS results across diverse populations. Summary statistics from each GWAS were meta-analyzed using MR-MEGA without combining by ancestry first. The MR-MEGA analysis was performed across four axes of genetic variation derived from pairwise allele frequency differences, based on the recommendation for separating major global ancestry groups. Index variants from the MR-MEGA analysis were selected using the same clumping parameters as described above (LD r² < 0.01 within a ±10-Mb window), based on the merged 1000 Genomes ALL reference panel. For each variant, we report two heterogeneity P values: one that is correlated with ancestry and accounted for in the meta-regression (P_Het-Anc) and the residual heterogeneity that is not due to population genetic differences (P_Het-Res).

PGS_PSA development and validation

We implemented two strategies for generating a genetic score for PSA levels. In the first approach, we selected 128 variants that were genome-wide significant (P < 5 × 10⁻⁸) in the multi-ancestry meta-analysis and were independent (LD r² < 0.01 within a ±10-Mb window) in 1000 Genomes EUR and (LD r² < 0.05) 1000 Genomes ALL populations (PGS₁₂₈). Each variant in PGS₁₂₈ was weighted by the meta-analysis effect size estimated using METAL. As an alternative strategy to clumping and thresholding, we fit a genome-wide score using the PRS-CSx algorithm³⁶, which takes GWAS summary statistics from each ancestry group as inputs and estimates posterior SNP effect sizes under coupled continuous shrinkage priors across populations (PGS_CSx). Analyses were conducted using pre-computed population-specific LD reference panels from the UKB, which included 1,287,078 HapMap3 variants that are available in both the UKB and 1000 Genomes phase 3.

We calculated a single trans-ancestry PGS that can be applied to all participants in the target cohort, rather than optimizing a PGS within each ancestry group. This approach is more robust to differences in genetic ancestry assignments across studies and does not require separate testing and validation datasets for parameter tuning each ancestry group³⁶. To facilitate this type of analysis, PRS-CSx provides a –meta option that integrates population-specific posterior SNP effects using an inverse-variance-weighted meta-analysis in the Gibbs sampler³⁶. The global shrinkage parameter was set to φ = 0.0001. PRS-CSx was run on the intersection of variants that were in the LD reference panel and had imputation quality (INFO > 0.90), resulting in 1,058,163 variants in PCPT and 1,071,268 variants in SELECT. Because PRS-CSx considers only autosomes, chrX variants that were included in PGS₁₂₈ were added to PGS_CSx separately, when output files from each chromosome produced by the PLINK–score command were concatenated.

The predictive performance of PGS_CSx and PGS₁₂₈ was evaluated in two independent cancer prevention trials that were not included in the meta-analysis: PCPT and SELECT. Analyses were conducted in the pooled sample for each cohort, which included individuals of all ancestries who passed QC filters (Supplementary Note). Ancestry-stratified analyses were conducted for clusters with n > 50 with available genetic ancestry PCs. Ancestry scores were computed with SNPWEIGHTS⁷⁸. Individuals with ancestry scores ≥0.80 for a single group were assigned to clusters for predominantly European (EUR), West African (AFR) and East Asian (EAS) ancestry. Admixed individuals with intermediate ancestry scores for at least one group were assigned to separate clusters: 0.20 < EUR/AFR < 0.80 or 0.20 < EUR/EAS < 0.80. Pooled analyses were adjusted for ten within-cluster PCs and global ancestry proportions (AFR and EAS).

Index event bias analysis

Index event bias occurs when individuals are selected based on the occurrence of an event or specific criterion. This is analogous to the direct dependence of one phenotype on another, as in the commonly used example of cancer survival³⁴. Due to unmeasured confounding, this dependence can induce correlations between previously independent risk factors among those selected^33,34. Genetic effects on prostate cancer can be viewed as conditional on PSA levels, because elevated PSA typically triggers diagnostic investigation. Genetic factors resulting in higher constitutive PSA levels may also increase the likelihood of prostate cancer detection due to more frequent testing (Fig. 4). This selection mechanism could bias prostate cancer GWAS associations by capturing both direct genetic effects on disease risk and selection-induced PSA signals. In the GWAS setting, methods using summary statistics have been developed to estimate and correct for this bias^33,35. Although typically derived assuming a binary selection trait, these methods are still applicable to selection or adjustment based on quantitative phenotypes³³. In this study, we conceptualized PSA variation as the selection trait and prostate cancer incidence as the outcome trait (Fig. 4).

We applied the method described in Dudbridge et al.³³, which tests for index event bias and estimates the corresponding correction factor (b) by regressing genetic effects on the selection trait (PSA) against their effects on the subsequent trait (prostate cancer), with inverse variance weights: w = 1/(SE_PrCa)². Summary statistics for prostate cancer were obtained from the most recent prostate cancer GWAS from the PRACTICAL consortium³². Sensitivity analyses were performed using SlopeHunter³⁵, an extension of the Dudbridge approach that allows for direct genetic effects on the index trait and subsequent trait to be correlated. For both methods, analyses were conducted using relevant summary statistics and 127,906 variants pruned at the recommended threshold³³ (LD r² < 0.10 in 250-kb windows) with MAF ≥ 0.05 in the 1000 Genomes EUR reference panel. After merging the pruned 1000 Genomes variants with each set of summary statistics, variants with large effects, (|β| > 0.20) on either log(PSA) or prostate cancer, were excluded. The resulting estimate (b), adjusted regression dilution using the SIMEX algorithm, was used as a correction factor to recover unbiased genetic effects for each variant:$\beta _{PrCa}^\prime$ = β_PrCa−b×β_PSA, where β_PSA is the per-allele effect on log(PSA), andβ_PrCa is the log(OR) for prostate cancer.

The impact of the bias correction was assessed in three ways. First, genome-wide significant prostate cancer index variants were selected from the European ancestry PRACTICAL GWAS meta-analysis (85,554 cases and 91,972 controls) using clumping (LD r² < 0.01 within 10 Mb) (ref. ³²). We tabulated the number of variants that remained associated at P < 5 × 10⁻⁸ after bias correction. Next, we fit genetic scores for PSA and prostate cancer in men of European ancestry in the UKB who were not included in the PSA or prostate cancer GWAS (11,568 prostate cancer cases and 152,884 controls). We compared the correlation between the PGS for PSA (PGS_PSA), comprising 128 lead variants, and the 269-variant prostate cancer risk score fit with original risk allele weights (PGS₂₆₉) and with weights corrected for index event bias (PGS₂₆₉^adj). To allow adjustment for genetic ancestry PCs and genotyping array, associations between the two scores were estimated using linear regression models. Next, we examined associations for each genetic score (PGS₂₆₉, PGS₂₆₉^adj, PGS₂₆₉^adj-S) with prostate cancer in a subset of GERA participants who underwent a biopsy. Because GERA controls were included in the PSA GWAS meta-analysis, AUC estimates and corresponding bootstrapped 95% CIs were obtained using tenfold cross-validation. We also examined PGS associations with Gleason score, a marker of disease aggressiveness, which was not available in the UKB. Multinomial logistic regression models with Gleason score ≤6 (reference), 7 and ≥8 as the outcome were fit for each score in 4,584 cases from the GERA cohort.

Application of genetically adjusted PSA for biopsy referral and prostate cancer detection

Genetically corrected PSA values were calculated for individual i as follows^17,19:

$$PSA_i^G = \frac{{PSA_i}}{{a_i}}$$

(1)

where a_i is a personalized adjustment factor derived from PGS_PSA. Because genetic effects were estimated for log(PSA), a_i for correcting PSA in ng ml⁻¹ was derived as:

$$a_i = \frac{{\exp \left( {PGS_i} \right)}}{{\exp ( {\overline {PGS} } )}}$$

(2)

$\overline {PGS}$ can be estimated in controls without prostate cancer or obtained from an external control population^17,19. We see that a_i > 1 when an individual has a higher multiplicative increase in PSA than the sample average due to their genetic profile, resulting in a lower genetically adjusted PSA compared to the observed value ($PSA_i^G < PSA_i$).

We evaluated the potential utility of PGS_PSA in two clinical contexts. First, we quantified the impact of using $PSA_i^G$ on biopsy referrals by examining reclassification at age-specific PSA thresholds used in the Kaiser Permanente health system. Analyses were conducted in GERA participants with information on biopsy date and outcome, comprising prostate cancer cases not included in the PSA GWAS and controls that were part of the PSA GWAS. To use the same normalization factor for both cases and controls while mitigating bias due to control overlap with the PSA discovery GWAS, a_i for GERA participants was calculated by substituting $\overline {PGS}$ from out-of-sample UKB controls (n = 152,884). Upward classification resulting in biopsy eligibility occurred when $PSA_i^G > PSA_i \cap PSA_i^G > ref$, where ref is the biopsy referral threshold. Downward classification resulting in biopsy ineligibility was defined as: $PSA_i^G < PSA_i \cap PSA_i^G < ref$. Net reclassification (NR) was summarized separately for cases and controls:

$$NR_{{{{\mathrm{case}}}}} = P\left( {up{{{\mathrm{|case}}}}} \right) - P\left( {down{{{\mathrm{|case}}}}} \right)$$

$$NR_{{{{\mathrm{control}}}}} = P\left( {down{{{\mathrm{|control}}}}} \right) - P\left( {up{{{\mathrm{|control}}}}} \right)$$

This is equivalent to tabulating the proportion of individuals in each biopsy eligibility category:

$$NR_{{{{\mathrm{case}}}}} = \left( {\frac{{n_{{{{\mathrm{eligible}}}}}}}{{n_{{{{\mathrm{case}}}}}}}} \right) - \left( {\frac{{n_{{{{\mathrm{ineligible}}}}}}}{{n_{{{{\mathrm{case}}}}}}}} \right)$$

$$NR_{{{{\mathrm{control}}}}} = \left( {\frac{{n_{{{{\mathrm{ineligible}}}}}}}{{n_{{{{\mathrm{control}}}}}}}} \right) - \left( {\frac{{n_{{{{\mathrm{eligible}}}}}}}{{n_{{{{\mathrm{control}}}}}}}} \right)$$

For each NR proportion, 95% CIs were obtained using the normal approximation:

$$NR \pm 1.96 \times \sqrt {\frac{{\left| {NR} \right| \times \left( {1 - \left| {NR} \right|} \right)}}{n}}$$

Next, we assessed the performance of risk prediction models for prostate cancer overall, aggressive prostate cancer and non-aggressive prostate cancer in the PCPT and the SELECT.

Because both studies were excluded from the PSA GWAS meta-analysis, a_i and $PSA_i^G$ for then PCPT and the SELECT were calculated using $\overline {PGS}$ observed in each respective study. Consistent with the PGS_PSA validation analysis, pooled analyses included individuals of all ancestries who passed QC filters. To facilitate ancestry-stratified analyses in SELECT, especially for aggressive disease, we combined AFR and AFR/EUR clusters into a single group (AFR pooled) and similarly pooled EAS and EAS/EUR (EAS pooled). Aggressive prostate cancer was defined as Gleason score ≥7, PSA ≥ 10 ng ml⁻¹, T3–T4 stage and/or distant or nodal metastases. We compared AUC estimates for logistic regression models using the following predictors, alone and in combination: baseline PSA, genetically adjusted baseline PSA (PSA^G) PGS_PSA, prostate cancer risk score with original weights (PGS₂₆₉) (ref. ³²) and weights corrected for index event bias (PGS₂₆₉^adj).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.