GPS derivation and selection

To calculate a GPS, we used a recently developed computational algorithm, LDPred [11], to reweight each variant and obtain the average effects for each of 103,155 genetic variants on fracture and estimated BMD from the most extensive GWAS study of osteoporosis published to date [4]. This algorithm can reweight each variant according to the given effect size of the prior distribution and incorporate more variants observed in the prior GWAS. Moreover, with the input of a comprehensive reference panel, this method can incorporate the degree of correlation between a variant and others nearby and a tuning parameter that denotes the proportion of variants with non-zero effect size [11]. Vilhjálmsson et al. [11] recommended testing a range of different tuning parameters to capture the non-zero effect size. In our study, we tested the seven candidate tuning parameter values in order to identify the best GPS.

We tested the seven candidate-GPSs with measured hip and spine BMD in a validation dataset of 3739 women from the Women’s Health Initiative (WHI) Genomics and Randomized Trials Network (GARNET) study. The WHI GARNET study aimed to identify genetic factors affecting myocardial infarction, stroke, venous thrombosis, and diabetes phenotypes through genome-wide analysis using a nested case-control study design [32]. Of the 3739 women in the GARNET study, 2458 participants had genotype and BMD measurements available. All 2458 participants were Non-Hispanic White population. We compared the maximal correlation coefficient with each BMD to select the best GPS. Each candidate-GPS was strongly correlated with measured BMD at both spine and hip (all p < 0.001). The correlation coefficients ranged from − 0.22 to 0.21 for hip BMD (Additional file 1: Figure S1) and from 0.22 to 0.21 for spine BMD (Additional file 1: Figure S2), respectively. Similar results were obtained for each GPS after the adjustment of genetic background, which was evaluated and quantified by the top ten principal components of ancestry.

The highest correlation with BMD in both spine (− 0.22) and hip (− 0.22) was achieved with the GPS (\(\rho =0.03\)); hence we used the GPS (\(\rho =0.03\)) as the best one to move forward in the analysis of the testing dataset. The best-performing GPS (ρ = 0.03) contains all 103,155 variants. The GPS (ρ = 0.03) is normally distributed with the empirical risk of fracture (Fig. 3). The median GPS (ρ = 0.03) percentile score was \(-25\) for individuals without the fracture vs. 42 for individuals with the fracture. The testing dataset (n = 15,776) is independent of the validation dataset (n = 2458) studied above. Additional details of GPS derivation and validation are shown in Fig. 1.

Fig. 1

Study Overview. Derivation, Validation, and Testing of a Genome-wide Polygenic Score (GPS) for Osteoporosis. A genome-wide polygenic score for osteoporosis was used, including 103,155 variants from the previously published [4] GWAS study. The LDPred computational algorithm [11] was used to calculate GPS vary for the tuning parameter, \(\rho \), the fraction of causal variants, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, and 1. A conventional genome-wide significant variants score with conditionally independent 1,103 variants from the previously published GWAS study [4]

Our GPS of 103,155 variants showed significantly higher predictive power than the conventional GRS, which comprises the conditionally independent 1103 SNPs [4] from a comprehensive GWAS study. Within the validation dataset (\(n=2458)\), correlation with each BMD for GRS derived from 1103 SNPs [4] ranged from − 0.146 to − 0.159. This lower strength of association with BMD using only the conditionally independent, genome-wide significant 1103 SNPs was consistent with the previous study [33], ranging from − 0.154 to − 0.188. Having derived and validated a new polygenic score that considerably outperformed the conventional GRS calculated from the 1103 SNPs, we further studied the predictive power of the new GPS on major osteoporotic fracture, weight, and osteoporosis in an independent testing dataset (n = 15,776).

Polygenic susceptibility to major osteoporotic fracture (MOF) and osteoporosis

We examined the extent to which the GPS predicted MOF and osteoporosis in an independent testing dataset (total n = 15,776) of WHI three sub-studies, including Integrative genomics and risk of coronary heart disease and related phenotypes, Population Architecture using Genomics and Epidemiology (PAGE), and Women’s Health Initiative Memory Study (WHIMS). This testing dataset is independent of the validation dataset (n = 2458) used earlier in this study. We also investigated the transferability of GPS in the independent testing dataset of 15,776 participants stratified by race/ethnicity: the Non-Hispanic White (n = 3427), Black or African-American (n = 9742), Hispanic/Latino (n = 2929) and others (n = 429, including American Indian or Alaskan Native, and Asian or Pacific Islander).

The baseline characteristics of the participants in an independent testing dataset were stratified by MOF status (Table 1). About 6% (n = 941) of participants had at least one major osteoporotic fracture (MOF) during an average of 12 years of follow-up. The participants with MOF have significantly higher GPS (p < 0.001), are older, have a lower weight, and decreased BMD in both hip and spine. 50% of participants with MOF are Not-Hispanic White, and 59% without MOF are Black or African-Americans. Rheumatoid arthritis, previous fragility fracture, previous osteoporosis, and parental fracture history significantly differ between MOF and non-MOF participants. In particular, the participants with MOF have higher previous osteoporosis and parental fracture history. 46.9% of participants in the testing dataset had normal BMD (T-score \(\ge -1.0\)), 37.2% of the participants were osteopenia (\(-2.5<\) T-score \(< -1.0\)), 13.6% of the participants were osteoporosis (T-score \(\le -2.5\)), and 2.3% met criteria for severe osteoporosis (T-score \(\le -2.5\) and presence of one or more fragility fractures).

Table 1 Baseline descriptive statistics of 15,776 women in an independent testing dataset stratified by Major osteoporotic fracture (MOF) status

The baseline characteristics of the participants in the independent testing dataset were also stratified by race/ethnicity: Non-Hispanic White (n = 3427), Black or African-American (n = 8991), Hispanic/Latino (n = 2,929), and others (including American Indian or Alaskan Native, and Asian or Pacific Islander) (n = 429). (Additional file 1: Table S1). Among the participants who had at least one MOF in the independent testing dataset (n = 941), the Non-Hispanic White participants had a significantly higher percentage of MOF (49.9%), followed by Black or African-American (26.5%), Hispanic/Latino (17.9%), and others (5.7%). The participants in the Non-Hispanic White group are older and have higher GPS, lower weight, and lower BMD in both hip and spine compared with other racial and ethnic groups.

The GPS approximated a normal distribution in the independent testing dataset (Fig. 3A). The correlation between the GPS and observed BMD (hip and spine) ranged from − 0.225 to − 0.218. These results are similar to our observation in the validation dataset. We then stratified the participants in the testing dataset according to GPS decile and found a remarkable gradient with respect to BMD, MOF, and body weight (Fig. 2A–C). For each decile group, we calculated the mean for the continuous phenotype variable and frequency for the categorical variable. For example, the mean hip BMD was 0.862 g/cm2 for those in the top decile of the GPS and 0.931 g/cm2 for those in the bottom decile, a difference of 0.069 g/cm2 (p < 0.001). Similarly, the average body weight was 76.5 kg for those in the top decile of the GPS and 82.2 kg for those in the bottom decile, a difference of 5.8 kg (p < 0.001). MOF was present in 139 of 941 (14.7%) in the top decile of the GPS versus 27 of 941 (2.8%) in the bottom decile, corresponding to a fivefold gradient in fracture risk (p < 0.001).

Fig. 2

Relationship of a Genome-wide Polygenic Score (GPS, LDPred with \(\rho =0.03\)) distribution in the testing dataset (\(n=\mathrm)\) with A Hip bone mineral density, B Major osteoporotic fracture (MOF), and C Body weight. Significant differences in osteoporosis categories were observed D when participants were stratified into the bottom decile, deciles 2–9, and top decile

Despite the strong associations observed in this study, polygenic susceptibility of the GPS to osteoporosis is not deterministic. Among those in the top decile of the GPS, 57.5% of participants in the testing dataset were osteoporosis and severe osteoporosis (Fig. 2D). In contrast, among those in the bottom decile of the GPS, 29.2 % of participants were osteoporosis and severe osteoporosis. However, among those in the top decile of the GPS, 15.8% had a normal range (Fig. 2D). These results were very similar after adjusting the top ten principal components.

A high polygenic score is common among those with severe osteoporosis

Conventional analyses of rare genetic mutations are conducted by comparing heterozygous mutation carriers with non-carriers. Individuals carrying the variants within or close to the LRP5, SOST, OPN, and TNFRSF11A genes had a significantly higher fracture risk, with odds ratios ranging from 1.13 and 1.43 per allele [34]. We tried to mimic this method using the new GPS by labeling the top decile of the GPS distribution as “carriers” and those in the remainder of the distribution as non-carriers (Fig. 3A). The magnitude of risk conferred by a high GPS increased at more extreme levels of observed disease risk. The proportion of high-GPS carriers was 15.8% among the 7410 individuals with normal BMD, 26.7% among the 5872 individuals with osteopenia, 29.8% among the 2140 individuals with osteoporosis, and 27.7% among the 354 individuals with severe osteoporosis. Compared with the remainder of the GPS distribution, the top GPS decile was associated with a 15.25-, 3.62-, and 1.89-fold increased risk of severe osteoporosis, osteoporosis, and osteopenia, respectively (Fig. 3B). Using the same method, we calculated the odds ratio of various fracture types for the top GPS decile versus the 90% remainder of the distribution with adjustment of clinical factors including age, height, weight, parental fracture, previous fracture, smoking, glucocorticoid use, rheumatoid arthritis, BMD (hip or spine), and previous osteoporosis. The results show that compared with the remainder of the GPS distribution, the top GPS decile was associated with a 3.59-, 2.48-, 1.92-, and 1.58-fold increased risk of any-fracture, MOF, hip fracture, and spine fracture, respectively (Fig. 3B).

Fig. 3

Association of High Genome-wide Polygenic Score (GPS, LDPred with \(\rho =0.03\)) with osteoporosis category and fracture status in the testing dataset (\(n=\mathrm)\). A The top 10% of the GPS distribution was considered high GPS “Carriers,” shaded in red, compared to the remaining 90%. B The relationship of high GPS to the osteoporosis category and fracture status was measured using a multiple logistic regression. OR, odds ratio, CI, confidence interval, MOF, major osteoporotic fractures. 1) Osteopenia: T-score between − 1 and − 2.5. 2) Osteoporosis: T-score \(\le -2.5\). 3) Severe Osteoporosis: Reserved for patients with a fragility fracture(s) and a T-score \(\le -2.5\)

We further estimated the odds ratio and 95% CI for individuals in the top 30%, 20%, 10%, and 5% of the GPS compared with the remaining individuals (Table 2). The odds ratio per standard deviation increment of MOF risk in the top 5% and 10% of GPS distribution were 3.12 (95% CI 2.25–5.42, p < 0.001) and 2.48 (95% CI 1.86–2.76, p < 0.001), compared with the remaining 95% and 90% of the individuals, respectively. In contrast, the odds ratio per standard deviation increment of MOF risk in the top 5% and 10% of GRS distribution were 1.25 (95% CI 0.98–1.75, p < 0.13) and 0.97 (95% CI 0.67–1.68), respectively.

Table 2 The Odds Ratio (OR) estimate of Major osteoporotic fractures (MOF) is derived from two methods in an independent testing dataset (\(n=\mathrm\))

We also examined the transferability of GPS by comparing the odds ratio and its 95% CI stratified by four populations: Non-Hispanic White (n = 3427), Black or African-American (n = 8991), Hispanic/Latino (n = 2929) and others (including American Indian or Alaskan Native, and Asian or Pacific Islander) (n = 429). (Additional file 1: Table S2). The odds ratio per standard deviation increment of MOF risk in the top 5% of GPS distribution in the Non-Hispanic White, Black or African American, and Hispanic/Latino were 2.26 (95% CI 1.56–2.63, p < 0.001), 1.53 (95% CI 1.13–1.84, p < 0.001) and 1.19 (95% CI 1.01–1.47, p < 0.001), compared to the remaining 95% of GPS distribution, respectively.

Postmenopausal women’s risk of developing severe osteoporosis varies according to polygenic score

Although only a small percentage of postmenopausal women experienced severe osteoporosis in their middle age (around 50) at the baseline, the prevalence of severe osteoporosis increases substantially over subsequent decades at the 20-year follow-up. We hypothesized that the GPS might significantly predict who will develop severe osteoporosis during the transition from middle age to the elderly. Among individuals in the top decile of the GPS, 215 of 1565 (13.7%) developed severe osteoporosis compared with 6.2% of those in deciles 2–9 (Fig. 4). By contrast, among those in the lowest decile, only 35 of 12,522 (< 1%) individuals developed severe osteoporosis.

Fig. 4

Association of High Genome-wide Polygenic Score (GPS, LDPred with \(\rho =0.03\)) with incident severe osteoporosis. Among 15,776 individuals in an independent testing dataset, GPS was stratified into top decile, deciles 2–9, and bottom decile. Incident severe osteoporosis is plotted with three GPS group stratification

With respect to the discriminative capacity, we first evaluated a baseline model of age, height, and weight, yielding a C statistic of 0.681 (95% CI 0.674–0.689). Each of the nine additional risk factors was then added (individually) to this baseline model for the outcome of major osteoporotic fractures (Fig. 5). GPS had a higher discriminative capacity with a C statistic of 0.723 (95% CI 0.715–0.729). By contrast, the addition of GRS has a C statistic of 0.708 (95% CI 0.703–0.714).

Fig. 5

Discriminative capacity of the genome-wide polygenic score and clinical risk factors in the testing dataset (\(n=\mathrm)\). The C statistic estimate of major osteoporotic fractures is first obtained with a baseline model of age, height, and weight using a multiple logistic regression model. Next, the C statistic was calculated after the additional inclusion of individual clinical risk factors; parental fracture, previous fracture, smoking, glucocorticoid use, rheumatoid arthritis, hip bone mineral density, spine bone mineral density, genetic risk score (1103 SNPs from GWAS), and genome-wide polygenic score (LDPred with ρ = 0.03)