Introduction

Osteoporosis (OP) is a chronic disease characterized by low bone mineral density (BMD) and an increased risk of fracture. Due to the increasing aging population, OP becomes a major public health threat worldwide. In the US, OP affected 53.6 million adults in 2010 and was estimated to affect 70.8 million people in 20301. Osteoporotic fracture (OF) can occur with low trauma or other activities, such as minor falls, twisting and bending. As the most severe complication of OP, OF causes a long-term disability, a high mortality to afflicted individuals, and substantial economic burden on health care systems2. According to a recent investigation in U.S., the annual cost of OF was $5.1 billion, which is much higher than the annual cost of stroke and myocardial infarction3.

BMD is the density of bone mineral in bone tissue, which is highly heritable. To date, over 200 loci associated with variations in BMD have been identified through genome-wide association studies (GWASs) and their meta-analyses4,5,6,7,8,9. BMD is negatively associated with the risk of OF within ethnic groups. BMD is considered to be a key and major risk factor of OF and is widely used to construct fracture risk assessment models5,10,11,12. However, BMD is not the only risk factor of OF, since not all people with a low BMD would develop OF during their lives, and about two-thirds of individuals who ever suffered a fracture were not diagnosed with osteoporosis before13.

Besides BMD, several clinical risk factors (CRFs), such as female gender, advancing age, and low body weight, family fracture history, smoking, and alcohol addiction are confirmed to be significantly associated with OF14,15. Some of the CRFs (such as ethnicity, weight, and height) are highly heritable7,8,16,17,18,19. These CRFs are also commonly used as prediction factors in fracture risk assessment models20,21. As a phenotype largely affected by BMD and CRFs, fracture susceptibility is also considered to be partially genetically determined, and several studies based on genome-wide association study (GWAS) showed that genetic risk factors (GRFs) can help improve the accuracy of fracture prediction. For example, Ho-Le et al.12 showed that genetic profiling using results derived from populations of European descent can improve the accuracy of fracture risk assessment, especially non-hip fractures in Australians and Lee et al.22 showed that genetic factors identified from populations of European descent and East-Asian can improve prediction of nonvertebral fracture in postmenopausal women in South Korea.

Most of the risk factors and subjects used in previous studies were from people of European descent or European populations. However, genetic factors for BMD or bone disease vary greatly among different ethnicities23,24,25,26. Whether genetic factors identified from people of European descent can help improve the prediction of BMD and the accuracy of fracture in Chinese populations remains uncertain. Moreover, in practice, only BMD-associated genetic factors were identified in most of the osteoporosis genetics studies, but whether BMD-associated genetic factors in people of European descent can improve the fracture prediction power of osteoporotic fractures in Chinese populations is not clear. The aim of this study is to investigate (1) whether fracture-associated genetic factors identified from people of European descent studies can improve the risk prediction of OF in Chinese, (2) whether BMD-associated genetic factors identified from people of European descent can improve the prediction of BMD in Chinese, and (3) whether BMD-associated genetic factors identified from people of European descent can improve the fracture prediction power in Chinese.

Materials and Methods

Ethics Statement

This study was approved by the institutional review board of Xi’an Jiaotong University. This study incorporated samples from two researches6,9. All participants provided written informed consent. All methods in this study were carried out in accordance with relevant guidelines and regulations.

Study population

This study includes two Chinese Han ancestry samples. The inclusion and exclusion criteria for these two samples were detailed in our earlier publication6,9. Briefly, the first sample consists of 700 unrelated individuals, of which 350 are patients with osteoporotic (low trauma) hip fractures (including 226 females and 124 males) and 350 are elderly controls (including 177 females and 173 males) living in the city of Xi’an and its neighboring areas in China. This sample, with a case-control design, was initially used in a GWAS discovery stage for SNPs of the potential significance of OF in Chinese Han6. Individuals with low trauma hip fractures were recruited from affiliated hospitals and their associated clinics of Xi’an Jiaotong University. Healthy control subjects were obtained using local advertisements. Controls were geographic- and age-matched to the cases6. The second sample was derived from the China osteoporosis study (COS) comprised of 1620 unrelated individuals, whose BMD were measured at lumbar spine (LS), hip and/or femoral neck (FN) by dual-energy X-ray absorptiometry (DXA)9. Table 1 describes basic characteristics of the samples.

Table 1 Baseline characteristics of the two study samples.
Full size table

Data pre-processing

The first sample was derived from an in-house study performed by Guo et al.6. IMPUTE program27 was utilized to impute the genotypes of all SNPs located in the genomic regions of interest based on Asian HapMap data. SNPTEST27 was used to test the association between the imputed SNPs and OF. The second sample was derived from a GWAS discovery sample used by Zhang et al.9 and the SNPs were imputed by the 1000 genomes project (1KG) sequence variants (as of August 2010). Reference haplotypes representing 193 individuals with Asian ancestry were downloaded from the MACH website28. Prior to imputation, a consistency test of allele frequency between GWAS and reference samples was performed with chi-squared test9. To correct for potential mis-strandedness, GWAS SNPs that failed a consistency test (P < 1.0 × 10−6) were transformed into reverse strands. SNPs that further failed consistency were removed from the GWAS sample. Each GWAS sample was imputed by the respective reference panel with the closest ancestry9,29.

Selection of genes and clinical risk factors

In total, all the SNPs in 17 fracture-associated genes and 107 BMD-associated genes (including 74 LS-BMD and 82 FN-BMD associated genes) were analyzed in our study. These genes were derived from predominantly samples of European descent in GEFOS25 which was a meta-analysis of GWASs with multiple studies of populations across North America, Europe, East Asia and Australia for FN-BMD and LS-BMD performed by the Genetic Factors for Osteoporosis (GEFOS) Consortium. The 17 fracture-associated genes were directly selected from a recent meta-analysis of GWAS5. The 107 BMD-associated genes were obtained by two methods. Firstly, similar to the way of selecting fracture-associated genes, 62 genes (47 LS-BMD and 48 FN-BMD associated genes with 33 common genes) were directly selected from 56 BMD-associated loci identified in the meta-analysis of GWAS5. Secondly, we used gene-based analysis to identify and select BMD-associated genes from the dataset of the GEFOS2 study5. Briefly, we performed gene-based analysis in a Knowledge-based mining system for Genome-wide Genetic studies (KGG) (http://bioinfo.hku.hk/kggweb/) to identify the BMD-associated genes harboring multiple association SNP signals, using the Extended Simes procedure method (GATES)30. Appendix detailed the gene-based analysis. We used the threshold of P-value < 2 × 10−6 31 to obtain 62 BMD-associated genes, among which 45 were associated with FN, 40 were associated with LS, and 23 were overlapped between them. Here, 13 LS-BMD associated genes and 11 FN-BMD associated genes were overlapped between the two methods. CRFs available and used here consisted of gender, age, height, and weight. Tables 1 and 2 presented the basic characteristics of CRFs of the two samples and the selected genes.

Table 2 A list of genes selected in the analysis.
Full size table

Statistics

We used logistic regression and multiple linear regression models to assess the contribution of genetic factors to prediction of fracture risk (in the first sample) and BMD (in the second sample) respectively. 10-fold cross-validation was employed in the two samples (that is, the sample was randomly partitioned into 10 equal size subsamples; of the 10 subsamples, a single subsample was retained as the validation data for testing the model, and the remaining 9 subsamples are used as training data). We obtained a genetic risk score (GRS) using following steps:

  1. (1)

    Each SNP was coded as the number of “risk” alleles (0, 1, 2) in the sense that the allele was associated with the risk increment of fracture or lower BMD.

  2. (2)

    Lasso regression was used in the training data sets to capture the significant SNPs and their coefficients. GRS was constructed by summing all weighted scores which were calculated by multiplying the number of risk alleles with regression coefficients derived in the Lasso regression.

  3. (3)

    Logistic regression or multiple linear regression was used in training sets to estimate regression coefficients.

  4. (4)

    Predictive analysis was done in the testing data sets.

In order to investigate whether genetic factors can help improve the prediction of fracture risk and BMD and whether BMD-associated genetic factors can improve fracture prediction, we constructed three models in the first sample (Model I-I, Model I-II, Model I-III) and in the second sample (Model II-I, Model II-II, and Model II-III), respectively.

In the first sample:

Model I-I: \({\rm{Fracture}}={\beta }_{0}+\beta \cdot {\rm{CRFs}}\), which includes CRFs data only.

Model I-II: \({\rm{Fracture}}={\beta }_{0}+{\beta }_{1}\cdot {\rm{CRFs}}+{\beta }_{2}\cdot {{\rm{GRS}}}_{{\rm{F}}}\), which includes CRFs data and GRSF data

Model I-III: \({\rm{Fracture}}={\beta }_{0}+{\beta }_{1}\cdot {\rm{CRFs}}+{\beta }_{2}\cdot {{\rm{GRS}}}_{{\rm{B}}}\), which includes CRFs data and GRSB data of FN.

In the second sample:

Model II-I: \({\rm{BMD}}={\beta }_{0}+\beta \cdot {\rm{CRFs}}\), which includes CRFs data only.

Model II-II: \({\rm{BMD}}={\beta }_{0}+\beta \cdot {{\rm{GRS}}}_{B}\), which includes GRSB data (LS or FN) only.

Model II-III: \({\rm{BMD}}={\beta }_{0}+{\beta }_{1}\cdot \mathrm{CRFs}+{\beta }_{2}\cdot {{\rm{GRS}}}_{{\rm{B}}}\), which includes CRFs data and GRSB data (LS or FN).

where GRSB and GRSF are GRS for BMD and fracture, respectively.

The area under the receiver operating characteristic (ROC) curve (AUC)32 and the net reclassification improvement (NRI) index were used to assess the ability of the GRS to predict fracture with logistic regression in the first sample. We used Mean Squared Error (MSE) to measure the testing error and performed t-test to test the difference between pairwise models under multiple linear regression models for BMD in the second sample. In the first sample, we performed the reclassification analysis33. The predicted risk fracture was estimated for each individual using Model I-I, Model I-II and Model I-III in the reclassification analysis, and then similar to Lee et al.22, individuals were classified into three risk groups in terms of their estimated fracture probability (i.e., higher risk group (≥15%), middle risk group (10–15%), lower risk group (<10%)). We calculated the proportion of individuals who would be reclassified into the three risk groups using trained assessment models. If genetic factors are useful for fracture prediction, the probability of fracture estimated by the model with GRS (Model I-II and Model I-III) would be expected to increase for the fracture group, and to decrease for the non-fracture group in the model without GRS (Model I-I). NRI index could be used to measure the predictive improvement and is computed as:

$$\begin{array}{ccc}{\rm{N}}{\rm{R}}{\rm{I}} & = & Pr({\rm{u}}{\rm{p}}|{\rm{f}}{\rm{r}}{\rm{a}}{\rm{c}}{\rm{t}}{\rm{u}}{\rm{r}}{\rm{e}})-\,Pr({\rm{d}}{\rm{o}}{\rm{w}}{\rm{n}}|{\rm{f}}{\rm{r}}{\rm{a}}{\rm{c}}{\rm{t}}{\rm{u}}{\rm{r}}{\rm{e}})+Pr({\rm{d}}{\rm{o}}{\rm{w}}{\rm{n}}|{\rm{n}}{\rm{o}}{\rm{n}}{\textstyle \text{-}}{\rm{f}}{\rm{r}}{\rm{a}}{\rm{c}}{\rm{t}}{\rm{u}}{\rm{r}}{\rm{e}})\\ & & -\,Pr({\rm{u}}{\rm{p}}|{\rm{n}}{\rm{o}}{\rm{n}}{\textstyle \text{-}}{\rm{f}}{\rm{r}}{\rm{a}}{\rm{c}}{\rm{t}}{\rm{u}}{\rm{r}}{\rm{e}}).\end{array}$$

where (Pr(up|fracture) − Pr(down|fracture)) represents the probability increment (gain) on the fracture group of two comparison models, and (Pr(down|non-fracture) - Pr(up|non-fracture)) represents the probability decrement (gain) on the non-fracture group of two comparison models.

The statistical significance of NRI was estimated by a Z-test, a simple asymptotic test for the null hypothesis of NRI = 033. All the analyses were performed using MATLAB statistical software, and P ≤ 0.05 was considered statistically significant.

Results

The baseline characteristics of the two samples are shown in Table 1. All the data are presented as mean ± SD, or as numbers and percentages. In the first sample, there are 700 individuals whose clinical and FN-fracture data are available for analyses. The average age is 70.52, and there are 350 individuals with fracture in femoral neck. In the second sample, there are 1620 individuals whose clinical and LS-BMD and FN-BMD data are available. The average age is 34.46. The average LS-BMD and FN-BMD are 0.948 g/cm2 and 0.920 g/cm2, respectively.

Table 2 listed all the genes in the analysis. 62 BMD-associated genes (40 associated with LS- BMD, 45 associated with FN- BMD, and 23 overlapped between them) were selected by using the gene-based analysis30. The other 62 BMD-associated genes (47 associated with LS-BMD, 48 associated with FN- BMD, and 33 overlapped between them) were directly selected from meta-analysis of GWAS5. There are 17 common genes in the two methods, including 13 common LS- BMD associated genes and 11 common FN-BMD associated genes. Other genes from samples of European descent in GEFOS25 are not associated with the trait in gene-based method. At last, a total of 107 genes (including 74 LS-BMD associated genes and 82 FN-BMD associated genes) were selected for our analysis. Among the 17 fracture-associated genes, SPTBN1 is associated with LS-BMD, FUBP3 and C18orf19 are associated with FN-BMD, and other genes are associated with both LS-BMD and FN-BMD5.

Table 3 showed the AUCs of the models for FN-fracture in the first sample. The AUC value of the model with CRF (Model I-I) is 0.653. When the GRS of fracture was added to the model with CRF (Model I-II), the AUC value decreases to 0.587. The AUC value of the model with CRF and GRS of FN-BMD (Model I-III) is 0.588, which is 0.065 smaller than that of the model with CRF only and similar to that of the model with CRF and GRSF. We also presented the odds ratio and 95% confidence interval for each CRF and GRS (Table 4). It can be seen that adding genetic factors of fracture or genetic factors of BMD did not approve the prediction of fracture risk. The relative importance analysis suggested 62.3% and 61.8% of the explained variance in Model I-II and Model I-III, respectively, were slightly lower than 65.5% in Model I-I. Genetic risk factor of fracture or BMD accounted for only 1.4% and 1.3% of the observed variance in fracture risk, respectively. These data showed that GRF was not a good predictive factor and declined the model predictive performance; GRFs of fracture could not help improve the ability of FN-fracture prediction, meanwhile BMD-associated genetic factors could not improve the fracture prediction power either in our Chinese sample.

Table 3 AUC of the Models.
Full size table
Table 4 The odds ratio and 95% confidence interval for each CRF and GRS.
Full size table

Table 5 showed the results of the reclassification analysis which compares the predictability of model I-II (GRSF + model I-I) with that of model I-I (using CRF only). It could be seen that, compared to model I-I, the addition of GRS of fracture (Model I-II) reclassified none into the higher risk group and reclassified 4 of the 350 individuals with FN-fracture (1.14%) into the lower risk group, a gain of 1.14% (0% minus 1.14%). Model I-II reclassified 7 of the 350 individuals with no FN-fracture (2.00%) into the lower risk group and 1 with no FN-fracture (0.30%) into the higher risk group, a gain of 1.70% (2.00% minus 0.30%). In total, there is a net gain of 0.56% with model I-II compared with model I-I, but it is not statistically significant (P = 0.57).

Table 5 Reclassification analysis to determine the effectiveness of adding the genetic risk score of fracture to the model with clinical risk factors.
Full size table

Table 6 showed the results of the reclassification analysis which compares the predictability of model I-III (GRSB + model I-I) with that of model I-I. While compared to model I-I, model I-III reclassified 4 of the 350 individuals with FN-fracture (1.14%) into the higher risk group and reclassified 4 of the 350 individuals (1.14%) into the lower risk group, which means no gain (1.14% minus 1.14%). For the non- fracture group, Model I-III reclassified 5 of the 350 individuals with no FN-fracture (1.43%) into the lower risk group and 1 with no FN-fracture (0.30%) into the higher risk group, a gain of 1.13% (1.43% minus 0.30%). Thus, there is a net gain of 1.13% with model I-III compared with model I-I, but similar to model I-II, it was not statistically significant either (P = 0.29).

Table 6 Reclassification analysis to determine the effectiveness of adding the genetic risk score of FN-BMD to the model with clinical risk factors.
Full size table

We performed t-tests for MSE of the testing error (here, the MSE for the testing error is the MSE between the observed and predicted values) between pairwise models under the multiple linear regression models in the second sample. MSE in Model II-I is significantly lower than that in Model II-II (P < 0.0001 for LS-BMD and FN-BMD), while MSE in Model II-II is significantly greater than that in Model II-III (P = 0.002 for LS-BMD and P = 0.003 for FN-BMD). There is no significant difference between Model II-I and Model II-III (P = 0.15 for LS-BMD and P = 0.10 for FN-BMD). These results indicated that adding genetic risk factors into the model with CRF could not help improve the ability of BMD prediction in this population of Chinese Han ancestry.

Discussion and Conclusion

BMD, as the strongest predictor of fracture risk, has high genetic determinations. To date, over 200 loci associated with the variations in BMD have been identified through GWAS and meta-analyses2,5,6,7,8,9. However, most of these factors were from European descent populations and studies. Our study is to investigate whether genetic factors identified predominantly from people of European descent can help improve the prediction of BMD and the prediction accuracy of osteoporotic fracture in Chinese. Here, we conducted assessments for FN-fracture prediction in a Chinese Han ancestry population using CRFs, GRFs associated with fracture and/or BMD from European descent studies. We found that genetic factors from people of European descent could not contribute to FN-fracture prediction and BMD prediction in our Chinese sample. Comparing with GRFs of fracture, BMD-associated genetic factors could not help improve the accuracy of the fracture prediction either. We further analyzed assessment models for BMD prediction using BMD-associated GRFs in another Chinese Han ancestry population, but there is still no significant difference between the model with CRFs only and the model using CRFs and BMD-associated GRFs. These results suggested that BMD-associated genetic factors from European descent studies may not help improve BMD prediction for Chinese.

As we introduced before, there were limited studies12,22,23 supported that fracture-associated GRFs play an important role in fracture predicting. For instance, the genetic factors identified from populations of European descent and East-Asian can improve prediction of nonvertebral fracture in postmenopausal women in South Korea22 and the genetic profiling using results derived from populations of European descent can improve the accuracy of non-hip fractures assessment in Australian12. Ho-Le et al.12 showed that a significant association between GRS and fracture was observed for the vertebral and wrist fractures, but not for hip fracture. With GRS, the correct reclassification of fracture versus nonfracture ranged from 12% for hip fracture to 23% for wrist fracture. Although our study here showed that fracture-associated GRFs identified from largely European descent studies could not increase the predictive accuracy of FN-fracture, it is difficult to compare our finding with other studies because of several reasons. First, and perhaps foremost, BMD-associated or fracture-associated genes used in our analysis for Chinese population were not originally derived from Chinese populations, but initially identified largely among European descent populations. The diversity between populations of different ancestries might carry different effects of genetic profiles. For example, a recent study showed that there are differences in the underlying linkage disequilibrium (LD) structure and in the frequency of sequence variants between European and Asian populations23. The genetic diversity in different populations might influence the results of association studies24 and thus might further influence the results of fracture prediction. Genes or SNPs selected to construct GRS in our analysis are also different from those in other studies. For example, 62 SNPs in 29 genes from previous GWAS7 were adopted for fracture prediction in Dubbo in Ho-Le et al.12, and in Lee et al.22, 39 SNPs in 30 genes were selected for fracture prediction in postmenopausal women in Koreans. In our analysis, 17 genes associated with fracture were directly selected from a recent meta-analysis of genome wide association study5 and 107 genes (including 74 LS- BMD associated genes and 82 FN-BMD associated genes) were selected using the same methods as that of fracture-associated genes and the gene-based analysis. There are only eight common genes (AKAP11, ESR1, F2, RPS6KA5, SP7, SPTBN1, ZBTB40, and WNT4) between Ho-Le’s study12 and our study, and there are only 17 common genes (CTNNB1, MEF2C, MEPE, SOST, SOX6, TNFRSF11B, ZBTB40, JAG1, LRP5, MARK3, SPTBN1, TNFRSF11A, ARHGAP1, C6orf97, DCDC5, ESR1, SP7) between Lee’s study22 and our study. Second, most of the previous studies directly selected the associated SNPs for prediction, while our study used the LASSO algorithm to select all the SNPs in the associated genes in the training model and then performed testing the prediction model. In fact, we found that directly using the associated SNPs for prediction achieved lower accuracy than that using the LASSO algorithm (data not shown). The reason is that we do not know which SNPs are actually associated with the BMD value and fracture in our Chinese samples. The LASSO algorithm can automatically capture the significant SNPs in our selected genes and thus can improve the accuracy of the prediction. Third, the fracture-associated genes were initially identified to be associated with fracture at any body site5, but in this study, the fracture prediction power we analyzed is for a homogeneous fracture type at skeletal site femoral neck only. Fourth, recent studies have shown that the utility of models for risk prediction depends on the sample size in training data sets34. The sample size in our study is round 1,000, which might to some extent limit the capability of the assessment models. Actually, all the subjects in our study were drawn from a homogenous population with the homogeneous phenotype of hip fractures. However, it took us several years to recruit such a sample because of high mortality rates after the hip fractures. We are now keeping the recruitment of hip OF subjects. OF Research with a larger sample awaits to be done in the future.

Because BMD is the strongest predictor for fracture risk, a natural idea could be that variants associated with BMD should contribute to fracture prediction. Indeed, in the Dubbo Osteoporosis Epidemiology Study, Ho-Le et al.12 provided evidence that BMD-associated genetic variants could improve the accuracy of fracture prediction, comparing with that using clinical risk factors alone. However, our study in the first sample of Chinese Han ancestry does not support this hypothesis. We found out that BMD-associated genetic factors do not contribute more to fracture prediction than clinical risk factors alone. Our analyses in the second sample of Chinese Han ancestry might partially explain these results. There is no significant difference between the prediction performance of the model using CRFs only and that of the model using CRF and GRFs of LS-BMD or FN-BMD, which indicates that genetic factors associated with BMD from people of European descent may not contribute to BMD prediction in our Chinese sample, and thus do not contribute more to fracture prediction with clinical risk factors in subjects from our Chinese sample. It should be noted that the second population here had an average age of 34. Although the BMD-associated genes were derived from populations with all ages5, our study was set up reasonably because populations where BMD-associated genes came from included the young. We also noted that the CRFs used in other studies for fracture prediction included BMD, fracture history, falls history etc.12,13,14. Although the CRFs in our study are just basic anthropomorphic variables, it did not affect our analysis because our aim was to assess whether genetic factors can help improve the prediction of BMD and the prediction accuracy of fracture. With additional useful CRFs included, the relatively performance of models with both CRFs and GRFs may be further decreased relative to that of using CRFs alone.

Indeed in the general studies of how to improve risk prediction, many efforts have repeatedly shown only moderate success in improving predictive accuracy by adding genetic information35,36,37,38. Several studies in assessing the incremental value of a limited number of SNPs compared to traditional risk factors showed that AUC has just achieved very small or minimum improvements39,40,41,42. For example, in the study of prediction of coronary heart disease (CHD), de Vries et al.43 found that the AUC improved from 0.716 to 0.718 and 0.719 using 49 SNPs significant at p < 5 × 10−8 and 152 SNPs significant at approximately p < 10−3, respectively. The AUC even was observed some decrease in the study of Morris et al.36 (as also seen in our study here). The reason of achieving small improvement or decrease of predictive ability, as Dudbridge et al.35 point out, is that nonpolygenic risk factors are themselves heritable, and may therefore mediate some or all of the genetic risk and/or in addition as we postulate, is that heterogeneity in the study populations and/or minuscule effects of the GWAS findings may further decreased the prediction power.

In conclusion, this study unfortunately does not support that fracture-associated genetic variants from people of European descent could help improve the FN-fracture prediction accuracy of clinical risk factors. Meanwhile BMD-associated genetic variants from people of European descent could not help improve the prediction power of FN-fracture and BMD in the two Chinese populations. More work in Chinese populations awaits to be done for uncovering genetic risk factors that are useful for disease prediction in Chinese. For example, large GWASs should be designed to discover genetic risk variants associated with BMD and fracture in Chinese, and more factors such as the family history of fracture and high-sensitivity C-reactive protein44,45 may be considered in the further study of fracture prediction.