Abstract
Objectives
Osteoporosis is a chronic disease. The aim of this study was to identify key genes in osteoporosis.
Methods
Microarray data sets GSE56815 and GSE56814, comprising 67 osteoporosis blood samples and 62 control blood samples, were obtained from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were identified in osteoporosis using Limma package (3.2.1) and Meta-MA packages. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed to identify biological functions. Furthermore, the transcriptional regulatory network was established between the top 20 DEGs and transcriptional factors using the UCSC ENCODE Genome Browser. Receiver operating characteristic (ROC) analysis was applied to investigate the diagnostic value of several DEGs.
Results
A total of 1320 DEGs were obtained, of which 855 were up-regulated and 465 were down-regulated. These differentially expressed genes were enriched in Gene Ontology terms and Kyoto Encyclopedia of Genes and Genomes pathways, mainly associated with gene expression and osteoclast differentiation. In the transcriptional regulatory network, there were 6038 interactions pairs involving 88 transcriptional factors. In addition, the quantitative reverse transcriptase-polymerase chain reaction result validated the expression of several genes (VPS35, FCGR2A, TBCA, HIRA, TYROBP, and JUND). Finally, ROC analyses showed that VPS35, HIRA, PHF20 and NFKB2 had a significant diagnostic value for osteoporosis.
Conclusion
Genes such as VPS35, FCGR2A, TBCA, HIRA, TYROBP, JUND, PHF20, NFKB2, RPL35A and BICD2 may be considered to be potential pathogenic genes of osteoporosis and may be useful for further study of the mechanisms underlying osteoporosis.
Cite this article: B. Xia, Y. Li, J. Zhou, B. Tian, L. Feng. Identification of potential pathogenic genes associated with osteoporosis. Bone Joint Res 2017;6:640–648. DOI: 10.1302/2046-3758.612.BJR-2017-0102.R1.
Article focus
-
The objective of this study was to identify key genes in osteoporosis.
Key messages
-
Genes including VPS35, FCGR2A, TBCA, HIRA, TYROBP, JUND, PHF20, NFKB2, RPL35A and BICD2 were considered to play potential roles in the pathology of osteoporosis.
Strengths and limitations
-
Bioinformatics analysis was used to identify and study the biological function of differentially expressed genes.
-
The validation sample for quantitative reverse transcriptase-polymerase chain reaction was small and it is necessary to collect larger samples for further validation.
Introduction
Osteoporosis, characterised by the impairment of bone microarchitecture and the loss of bone mass and strength, has become an important clinical problem in ageing populations.1,2 The spine is the most common site of osteoporotic fractures, followed by the hip, forearm and proximal humerus.3 Osteoporosis is characterised by deterioration of the microstructure of bone, specifically at trabecular sites, which leads to pain, deformity, disability and possibly death.4,5 A variety of factors contribute to the development of osteoporosis, such as genetic variants, gender, steroid production, age, lifestyle and environment.6-8 Additionally, low calcium intake, cigarette-smoking and intake of excessive alcohol may be secondary causes.4,9 Generally, osteoporosis is considered a silent disease because it is asymptomatic until a fracture occurs. Recently, the treatment of osteoporosis has been mainly pharmaceutical, but treatment may not be satisfactory due to its time-consuming nature and high cost, as well as the side effects of the drugs.
As Sims et al10 have reported, a number of studies have explored the mechanism of osteoporosis. The members of the Wnt signalling pathway, such as Wnt3a, low-density lipoprotein receptor-related protein 5, secreted frizzled-related protein 1 and sclerostin, have been suggested to be related to variation in bone mineral density (BMD).10 It is known that osteoporosis is defined clinically by measuring the BMD with heritability estimates of 0.5 to 0.9.11 This further demonstrated that BMD is an important clinical marker in osteoporosis. The underlying aetiology of osteoporosis is still not fully understood and the identification of novel therapeutic target genes for osteoporosis is needed.
It is worth mentioning that microarray data analysis is available to identify vital genes and gene regulatory networks associated with disease.12 In this study, we downloaded the data sets GSE56815 and GSE56814 from the Gene Expression Omnibus database, and identified the differentially expressed genes (DEGs) in blood samples obtained from osteoporosis patients, followed by Gene Ontology (GO)13 and Kyoto Encyclopedia of Genes and Genomes (KEGG)14 enrichment analyses and interaction network construction between transcription factors (TFs) and DEGs. We aimed to find involvement of key genes in osteoporosis which may be potential modulators in the pathology of osteoporosis.
Methods
Microarray data
The two microarray data sets (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56815 and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56814) were downloaded using the GEOquery package with R version (The R Foundation for Statistical Computing, Vienna, Austria) from the Gene Expression Omnibus (GEO) database. These two microarray data sets were from the blood samples of 67 osteoporosis patients (based on World Health Organization (WHO) criteria) and 62 control groups. It is reported that the WHO recommends the use of BMD of the spine and proximal femur measured by double energy X-ray absorption as the benchmark to diagnose osteoporosis and its severity.15 Therefore, diagnosis of osteoporosis depended on the BMD in these two datasets; these osteoporosis patients all had low BMD. In addition, the patients with osteoporosis did not receive any anti-osteoporosis medication or other medication that may affect bone metabolism. Detailed information of two data sets is shown in Table I.
Table I.
Two Gene Expression Omnibus (GEO) datasets of osteoporosis
| Gene type | GEO ID | Platform | Sample count (control group: patients with osteoporosis) | Notes |
|---|---|---|---|---|
| mRNA | GSE56815 | GPL96 [HG-U133A] Affymetrix Human Genome U133A Array | 80 (40:40) | Liu YZ, USA, 201616 |
| mRNA | GSE56814 | GPL5175 [HuEx-1_0-st] Affymetrix Human Exon 1.0 ST Array [transcript (gene) version] | 49 (22:27) | Liu YZ, USA, 201616 |
-
mRNA, messenger ribonucleic acid
Screening of DEGs
Normalised microarray data were downloaded and gene expression value was calculated as the mean value of its corresponding probe values. Limma package17 and Meta-MA18 package were used to identify the DEGs between osteoporosis and control group; p-values and false discovery rates (FDR) were further calculated. FDR < 0.01 was selected as the threshold for screening DEGs.
Functional analyses of DEGs
GO and KEGG enrichment analyses were carried out for the identified DEGs using GeneCodis3 (http://www.genecodis.cnb.csic.es/analysis).19-21 Significant GO terms and KEGG pathways were identified according to the threshold of p < 0.05.
Network construction of DEGs and TFs
To gain deeper insight into the molecular functions of DEGs, the gene regulatory relationships between DEGs and TFs were selected based on human TF binding sites data and genetic coordinate position information, which were available at the UCSC ENCODE Genome Browser.22 The identified TFs were considered to be associated with DEGs and the regulatory network between DEGs and TFs was visualised using Cytoscape software (Cytoscape Corporation, San Diego, California).
Quantitative reverse transcriptase-polymerase chain reaction (qRT-PCR) in vitro
Four women diagnosed with osteoporosis with low BMD were enrolled in this study, as were three women with high BMD but no diagnosis of osteoporosis. The clinical information of osteoporosis patients is shown in Table II. All blood samples were collected for further qRT-PCR experimentation. All participating individuals provided informed consent with the approval of the ethics committee.
Table II.
The clinical information of osteoporosis patients
| Patient number | Gender | Age (yrs) | Weight (kg) | Height (cm) | Smoking history | Drinking history | Family history | Menopause status | BMD | Diagnostic method | Symptom |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Female | 82 | 50 | 160 | No | No | No | Yes | -2.57 | DXA | Elderly osteoporosis |
| 2 | Female | 81 | 45 | 155 | No | No | No | Yes | -2.66 | DXA | Elderly osteoporosis |
| 3 | Female | 72 | 55 | 163 | No | No | No | Yes | -2.93 | DXA | Elderly osteoporosis |
| 4 | Female | 66 | 60 | 165 | No | No | No | Yes | -2.82 | DXA | Elderly osteoporosis |
| 5 | Female | 69 | 53 | 158 | No | No | No | Yes | -0.49 | DXA | Normal |
| 6 | Female | 79 | 46 | 162 | No | No | No | Yes | -0.58 | DXA | Normal |
| 7 | Female | 73 | 57 | 157 | No | No | No | Yes | -0.43 | DXA | Normal |
-
BMD, bone mineral density; DXA, double energy X-ray absorption
Total ribonucleic acid (RNA) of the blood samples was extracted using TRIzol Reagent (Invitrogen, Carlsbad, California), in accordance with the manufacturer’s protocols. Then 1 ug RNA was applied to synthesise DNA by SuperScript III Reverse Transcriptase (Invitrogen) and qRT-PCR was performed in an ABI 7500 Real-time PCR system (Invitrogen) with SYBR Green PCR Master Mix (Invitrogen). Glyceraldehyde 3-phosphate dehydrogenase was used as internal control and all reactions were performed in triplicate. Relative gene expressions were analysed by the 2-△△Ct method.
Receiver operating characteristic (ROC) analyses
By using the pROC package23 in R language, we performed the ROC analyses to assess the diagnostic value of DEGs. The area under the curve (AUC) under the binomial exact confidence interval was calculated and the ROC curve was generated.
Results
DEGs identification
In this study, we identified 1320 DEGs with significantly altered expression in osteoporosis blood samples, of which 855 were up-regulated and 465 were down-regulated compared with the control groups. The top 20 (ten up-regulated and ten down-regulated) DEGs abbreviations are defined in Table III and presented in Table IV. The results indicated that the expression pattern of these DEGs could observably distinguish the osteoporosis samples from control groups.
Table III.
Differentially expressed gene abbreviations
| Abbreviation | Full name |
|---|---|
| ALKBH1 | alkB homolog 1, histone H2A dioxygenase |
| BICD2 | BICD cargo adaptor 2 |
| CDKN2D | cyclin dependent kinase inhibitor 2D |
| DPP8 | dipeptidyl peptidase 8 |
| FCGR2A | Fc fragment of IgG receptor IIa |
| HIRA | histone cell cycle regulation |
| IER2 | immediate early response 2 |
| JUND | JunD proto-oncogene, AP-1 transcription factor subunit |
| METTL4 | methyltransferase like 4 |
| NBEAL2 | neurobeachin like 2 |
| NFKB2 | nuclear factor kappa B subunit 2 |
| NIF3L1 | NGG1 interacting factor 3 like 1 |
| NIT2 | nitrilase family member 2 |
| PAF1 | PAF1 homolog, Paf1/RNA polymerase II complex component |
| PHF20 | PHD finger protein 20 |
| POGLUT1 | protein O-glucosyltransferase 1 |
| RPL35A | ribosomal protein L35a |
| SAP130 | Sin3A associated protein 130 |
| SH3GLB2 | SH3 domain containing GRB2 like, endophilin B2 |
| TBCA | tubulin folding cofactor A |
| TYROBP | TYRO protein tyrosine kinase binding protein |
| VPS35 | VPS35, retromer complex component |
Table IV.
The top 20 differentially expressed genes (DEGs) in osteoporosis
| ID | Symbol | Combined ES | p-value | FDR | Regulation |
|---|---|---|---|---|---|
| 55737 | VPS35 | 1.41E+00 | 5.58E-13 | 4.37E-09 | Up |
| 56983 | POGLUT1 | 1.40E+00 | 7.55E-13 | 4.37E-09 | Up |
| 54878 | DPP8 | 1.39E+00 | 1.32E-12 | 5.09E-09 | Up |
| 60491 | NIF3L1 | 1.31E+00 | 1.32E-11 | 3.81E-08 | Up |
| 51230 | PHF20 | 1.25E+00 | 9.79E-11 | 2.03E-07 | Up |
| 56954 | NIT2 | 1.26E+00 | 1.05E-10 | 2.03E-07 | Up |
| 64863 | METTL4 | 1.21E+00 | 1.77E-10 | 2.92E-07 | Up |
| 79595 | SAP130 | 1.21E+00 | 3.95E-10 | 5.71E-07 | Up |
| 2212 | FCGR2A | 1.16E+00 | 1.32E-09 | 1.70E-06 | Up |
| 8846 | ALKBH1 | 1.17E+00 | 2.41E-09 | 2.79E-06 | Up |
| 4791 | NFKB2 | -1.06E+00 | 1.97E-08 | 6.92E-06 | Down |
| 6165 | RPL35A | -1.07E+00 | 2.51E-08 | 7.98E-06 | Down |
| 23299 | BICD2 | -1.05E+00 | 2.55E-08 | 7.98E-06 | Down |
| 6902 | TBCA | -1.07E+00 | 2.88E-08 | 8.35E-06 | Down |
| 56904 | SH3GLB2 | -1.03E+00 | 3.04E-08 | 8.58E-06 | Down |
| 7290 | HIRA | -1.02E+00 | 3.50E-08 | 9.22E-06 | Down |
| 54623 | PAF1 | -1.03E+00 | 3.76E-08 | 9.68E-06 | Down |
| 23218 | NBEAL2 | -1.05E+00 | 4.25E-08 | 1.07E-05 | Down |
| 9592 | IER2 | -1.02E+00 | 5.34E-08 | 1.19E-05 | Down |
| 1032 | CDKN2D | -1.07E+00 | 6.00E-08 | 1.24E-05 | Down |
-
ES, effect size; FDR, false discovery rate
GO and KEGG analyses of DEGs
Among 1320 DEGs, 1251 genes were recognised and significantly involved in different GO terms and KEGG pathways. The top 15 enriched GO terms, such as biological process, molecular function and cellular component of the identified DEGs are shown in Table V. The results showed that these DEGs were mainly involved in the GO terms associated with apoptosis processes, gene expression and signal transduction. On the other hand, the results of KEGG analysis revealed that a total of 15 pathways were enriched; for example, Ubiquitin mediated proteolysis and osteoclast differentiation (Table VI).
Table V.
Top 15 enriched Gene Oncology (GO) terms in osteoporosis
| GO ID | GO term | Genes (n) | FDR |
|---|---|---|---|
| Biological process | |||
| GO:0006915 | Apoptotic process | 75 | 2.26E-17 |
| GO:0010467 | Gene expression | 56 | 2.31E-14 |
| GO:0007165 | Signal transduction | 105 | 2.99E-14 |
| GO:0015031 | Protein transport | 50 | 9.06E-11 |
| GO:0006355 | Regulation of transcription, DNA-dependent | 119 | 1.16E-10 |
| GO:0006810 | Transport | 61 | 4.45E-10 |
| GO:0048011 | Nerve growth factor receptor signalling pathway | 33 | 1.19E-09 |
| GO:0016070 | RNA metabolic process | 36 | 1.73E-09 |
| GO:0016032 | Viral reproduction | 41 | 2.23E-09 |
| GO:0044419 | Interspecies interaction between organisms | 41 | 2.24E-09 |
| GO:0008380 | RNA splicing | 36 | 3.86E-09 |
| GO:0042981 | Regulation of apoptotic process | 31 | 4.62E-09 |
| GO:0000278 | Mitotic cell cycle | 38 | 8.32E-09 |
| GO:0016071 | mRNA metabolic process | 32 | 8.35E-09 |
| GO:0051437 | Positive regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle | 18 | 9.82E-09 |
| Molecular function | |||
| GO:0005515 | Protein binding | 446 | 2.10E-93 |
| GO:0000166 | Nucleotide binding | 180 | 7.43E-24 |
| GO:0003677 | DNA binding | 142 | 4.86E-16 |
| GO:0005524 | ATP binding | 126 | 5.21E-16 |
| GO:0046872 | Metal ion binding | 193 | 1.71E-14 |
| GO:0003723 | RNA binding | 69 | 5.81E-14 |
| GO:0016787 | Hydrolase activity | 90 | 1.13E-13 |
| GO:0016740 | Transferase activity | 59 | 2.82E-09 |
| GO:0008270 | Zinc ion binding | 129 | 3.94E-09 |
| GO:0003676 | Nucleic acid binding | 68 | 2.21E-08 |
| GO:0003700 | Sequence-specific DNA binding transcription factor activity | 71 | 1.04E-07 |
| GO:0016874 | Ligase activity | 40 | 1.10E-07 |
| GO:0003824 | Catalytic activity | 37 | 3.27E-06 |
| GO:0008233 | Peptidase activity | 44 | 4.47E-06 |
| GO:0042803 | Protein homodimerisation activity | 45 | 5.19E-06 |
| Cellular component | |||
| GO:0005737 | Cytoplasm | 464 | 3.22E-78 |
| GO:0005634 | Nucleus | 457 | 4.04E-71 |
| GO:0005829 | Cytosol | 246 | 6.44E-58 |
| GO:0005739 | Mitochondrion | 155 | 1.62E-33 |
| GO:0005654 | Nucleoplasm | 116 | 6.90E-31 |
| GO:0016020 | Membrane | 295 | 3.54E-30 |
| GO:0005730 | Nucleolus | 141 | 3.18E-24 |
| GO:0005622 | Intracellular | 156 | 2.02E-18 |
| GO:0005794 | Golgi apparatus | 97 | 7.96E-18 |
| GO:0005783 | Endoplasmic reticulum | 99 | 1.43E-17 |
| GO:0005789 | Endoplasmic reticulum membrane | 73 | 4.91E-16 |
| GO:0016021 | Integral to membrane | 259 | 8.22E-14 |
| GO:0005743 | Mitochondrial inner membrane | 43 | 2.07E-12 |
| GO:0005625 | Soluble fraction | 44 | 1.45E-09 |
| GO:0000139 | Golgi membrane | 46 | 1.55E-09 |
-
FDR, false discovery rate; RNA, ribonucleic acid; mRNA, messenger RNA; ATP, adenosine triphosphate
Table VI.
Top 15 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) terms in osteoporosis
| KEGG ID | KEGG term | Count | FDR | Genes |
|---|---|---|---|---|
| Hsa04120 | Ubiquitin mediated proteolysis | 24 | 2.83E-08 | UBE2A,MDM2,HERC1,UBE2Z,UBE4B,FZR1,DET1,CDC23,FBXO4,ANAPC2,ANAPC13,CDC27,UBE2W,UBE2G1,NEDD4,SOCS3,UBE2D3,UBE2E3,UBE4A,VHL,CUL1,UBE2D4,UBE2N,HERC2 |
| Hsa04380 | Osteoclast differentiation | 22 | 1.18E-07 | TYROBP,NFKB2,JUND,IL1R1,TYK2,NFKBIA,CREB1,TNFRSF1A,PPP3CC,LILRB3,SOCS3,FOSB,CYBB,FCGR2A,NFKB1,LILRA2,JUN,MITF,MAPK14,JUNB,PIK3CG,IFNAR1 |
| Hsa05016 | Huntington's disease | 26 | 1.32E-07 | POLR2J,NDUFA8,ATP5G1,COX6A2,POLR2H,POLR2B,NDUFC1,CREB1,COX5B,NRF1,GRIN2B,PLCB3,GNAQ,TBP,UQCRC2,NDUFA7,UQCR10,NDUFA13,BAX,NDUFB1,SDHD,PPARGC1A,NDUFB5,NDUFS3,SDHA,ATP5B |
| Hsa00190 | Oxidative phosphorylation | 22 | 1.45E-07 | NDUFA8,ATP5G1,COX6A2,NDUFC1,COX10,COX5B,ATP5L,ATP5I,ATP6V1E1,ATP6V0E2,UQCRC2,NDUFA7,UQCR10,NDUFA13,NDUFB1,SDHD,NDUFB5,ATP6V1B2,NDUFS3,SDHA,LHPP,ATP5B |
| Hsa04141 | Protein processing in endoplasmic reticulum | 24 | 1.76E-07 | SEC63,SSR2,UBE4B,CANX,ATF6,RPN2,SSR4,PPP1R15A,STT3A,UBE2G1,UBE2D3,UBE2E3,RRBP1,PDIA6,CUL1,BAX,UBE2D4,HERPUD1,SEL1L,DNAJC10,EIF2S1,LMAN2,MBTPS2,SIL1 |
| Hsa04142 | Lysosome | 20 | 5.87E-07 | NPC2,SCARB2,CTSC,GLA,LAPTM4A,MANBA,CLN5,FUCA1,AGA,LIPA,ASAH1,NPC1,PPT1,ARSB,HEXB,DNASE2,CTSO,IDS,IGF2R,LAPTM5 |
| Hsa05010 | Alzheimer's disease | 23 | 9.66E-07 | NDUFA8,ATP5G1,COX6A2,NDUFC1,IDE,ATF6,TNFRSF1A,COX5B,PPP3CC,GRIN2B,ITPR2,PLCB3,GNAQ,UQCRC2,NDUFA7,UQCR10,NDUFA13,NDUFB1,SDHD,NDUFB5,NDUFS3,SDHA,ATP5B |
| Hsa00020 | Citrate cycle (TCA cycle) | 10 | 1.50E-06 | CS,DLAT,IDH3A,SDHD,MDH1,DLST,ACO2,SDHA,PDHB,SUCLG2 |
| Hsa03040 | Spliceosome | 19 | 1.54E-06 | SF3A3,WBP11,PRPF6,RBM17,DHX38,PPIE,TXNL4A,PRPF38B,EIF4A3,PRPF3,RBM25,SART1,SRSF3,CDC5L,CRNKL1,ACIN1,LSM3,SNRNP40,SNRPB2 |
| Hsa03050 | Proteasome | 11 | 2.77E-06 | SHFM1,PSMC2,PSMD11,PSMC4,PSME3,PSMB3,PSMD1,PSMD2,PSME1,PSMD8,PSMA4 |
| Hsa00970 | Aminoacyl-tRNA biosynthesis | 10 | 2.27E-05 | IARS,DARS2,RARS,WARS,DARS,GARS,NARS2,NARS,WARS2,EPRS |
| Hsa05200 | Pathways in cancer | 31 | 2.43E-05 | VEGFB,MDM2,SOS1,NFKB2,ARNT,HRAS,CCNE1,NFKBIA,E2F1,CCND1,PTEN,TPM3,RALB,BCR,RALBP1,CTBP1,CTBP2,FZD2,NFKB1,TFG,JUN,VHL,FLT3,MAP2K2,MITF,BAX,VEGFC,SOS2,PIK3CG,MLH1,DVL1 |
| Hsa03013 | RNA transport | 19 | 2.53E-05 | EIF2S3,RPP30,EIF3H,NUP98,PRMT5,POM121,EIF3D,NUP43,NUP155,NXT2,EIF3G,RPP38,EIF4E2,EIF4A3,SNUPN,POP4,EIF2S1,ACIN1,EIF3E |
| Hsa05220 | Chronic myeloid leukemia | 13 | 3.34E-05 | MDM2,SOS1,HRAS,NFKBIA,E2F1,CCND1,BCR,CTBP1,CTBP2,NFKB1,MAP2K2,SOS2,PIK3CG |
| Hsa04144 | Endocytosis | 22 | 3.68E-05 | MDM2,CXCR4,SH3GLB2,ADRB2,RAB11FIP5,HRAS,FOLR2,ADRBK1,RAB11FIP3,GRK6,TSG101,ARAP1,NEDD4,HGS,PARD6A,SRC,RAB5C,RAB5A,ARFGAP1,CHMP2A,PDCD6IP,RNF41 |
-
FDR, false discovery rate; TCA, tricarboxylic acid; tRNA, trnaser ribonucleic acid; RNA, ribonucleic acid
Transcriptional regulatory relationships between DEGs and TFs
Regulatory relationships were predicted between top 20 DEGs (ten up-regulated and ten down-regulated) and TFs, and the regulatory network was established and visualised by Cytoscape software (Cytoscape Corporation) (Fig. 1). In this network, there were 6038 interactions pairs involving 88 TFs. The top seven TFs covering the most downstream genes were FOXD3, Nkx2-5, Pax-4, Oct-1, HNF-4, Pax-6 and COMP1.
Fig. 1
Transcription factors (TFs)-differentially expressed genes interaction networks. Rectangles and ellipses represent the TFs and target genes, respectively. The red and green colors represent up-regulation and down-regulation, respectively.
qRT-PCR
Among the identified DEGs, VPS35, FCGR2A, TBCA, HIRA, TYROBP and JUND were selected to verify the integrated result. The qRT-PCR results showed that TYROBP and JUND were up-regulated, while VPS35, FCGR2A, TBCA and HIRA were down-regulated. The expression of TBCA, HIRA and TYROBP were consistent with integrated analyses except VPS35, FCGR2A and JUND. The qRT-PCT results are shown in Figure 2.
Fig. 2
Validation of differentially expressed genes (DEGs) in the osteoporosis blood by quantitative reverse transcriptase-polymerase chain reaction. *p < 0.05, †p < 0.01.
ROC curve analysis
We performed ROC curve analyses and calculated the AUC to assess the discriminatory ability of DEGs (VPS35, HIRA, PHF20 and NFKB2) in data set GSE56815. The AUC of four DEGs including VPS35 (0.789), HIRA (0.77), PHF20 (0.851) and NFKB2 (0.741) was > 0.7 (Fig. 3). PHF20 had the largest AUC among these four DEGs. For the diagnosis of osteoporosis, the sensitivity (proportion of true positives) and 1-specificity (proportion of false positives) of VPS35 was 52.5% and 92.5%, respectively; the sensitivity and 1-specificity of HIRA was 60% and 87.5%, respectively; the sensitivity and 1-specificity of PHF20 was 77.5% and 85%, respectively; and the sensitivity and 1-specificity of NFKB2 was 67.5% and 75%, respectively.
Fig. 3
The receiver operating characteristic (ROC) curves of a) VPS35, b) HIRA, c) PHF20 and d) NFKB2 between osteoporosis patients and healthy controls. The ROC curves were used to show the diagnostic ability of these selected differentially expressed genes (DEGs) with 1-specificity (x-axis; the proportion of false positives) and sensitivity (y-axis; the proportion of true positives).
Discussion
Osteoporosis is a complex disease that is characterised by reduced bone mass and the deterioration of bone-tissue microarchitecture, ultimately leading to increased risk of fractures.24 In the present study, we identified 855 up-regulated and 465 down-regulated genes in blood samples of osteoporosis patients compared with control groups. It is well known that osteoclast formation and activation are critical events in maintaining the normal bone mass and structure. Herein, GO and KEGG analyses indicated that these DEGs were significantly involved in osteoclast differentiation, which demonstrated the important role in osteoclast development of osteoporosis. Additionally, the TF-DEGs regulatory network was constructed involving top 20 (ten up-regulated and ten down-regulated) genes and 88 TFs, which further illustrated the role of these DEGs under the regulation of TFs in osteoporosis. Finally, qRT-PCR in vitro validated the expression patterns of several genes including VPS35, FCGR2A, TBCA, HIRA, TYROBP and JUND. Some results have been inconsistent with the microarray analyses, probably because of the heterogeneity of the studies, including different inclusion criteria and the small number of patients in the validation set. In summary, ten genes including VPS35, FCGR2A, TBCA, HIRA, TYROBP, JUND, PHF20, NFKB2, RPL35A and BICD2 were considered to play a potential role in the pathology of osteoporosis.
It has been reported that VPS35 is highly expressed in osteoclasts as well as osteoblasts and loss of function will increase hyper-resorptive osteoclast formation.25 Specific VPS35 knockdown in the osteoblast lineage resulted in mildly lowered bone mass in the primary spongiosa.26 Our study found up-regulated expression of VPS35, which indicated that it may function in the regulation of osteoclast and osteoblast activity in osteoporosis. In addition, VPS35 had significant diagnostic value for osteoporosis, which may serve as a diagnostic biomarker of osteoporosis. FCGR2A has been reported to be related to rheumatoid arthritis and is a target gene of drugs in rheumatoid arthritis treatment.27,28 It has been indicated that FCGR2A is also associated with ankylosing spondylitis and axial spondyloarthritis.29 In the current study, we found an increased expression of FCGR2A in osteoporosis, suggesting that FCGR2A may be associated with the pathology of osteoporosis.
In human cell lines, TBCA knockdown will decrease a mass of a- and b-tubulin levels, subtle changes in the microtubule cytoskeleton and cell death.30 Herein, it is down-regulated in osteoporosis blood samples compared with controls. Therefore, we suggest that TBCA may participate in the formation process of cytoskeleton in osteoporosis. HIRA has been recorded associated with H3.3-containing nucleosomes in transcriptional active genomic loci in bone tumours.31 In this study, we found that HIRA was down-regulated in osteoporosis blood samples, which suggested that HIRA may play a significant role in histone modification in bone development in osteoporosis. In addition, HIRA has diagnostic value for osteoporosis, suggesting that HIRA may be associated with the pathology of osteoporosis and may serve as a diagnostic biomarker of osteoporosis.
TYROBP is a protein involved in osteoclast differentiation and function, such as the generation of the actin cytoskeleton, which is important for bone resorption.32 In the current study, the expression of TYROBP was up-regulated, indicating that it may play a vital role in the regulation of osteoclast differentiation in osteoporosis. JUND has proven to control bone formation, osteoblast proliferation, and differentiation.33-36 Furthermore, it is a DEG in mesenchymal stem cells in osteoporosis.37 Our results showed that JUND was expressed differentially in osteoporosis, which demonstrated the essential role of JUND in bone formation and osteoporosis.
It has been reported that PHF20-null mice show delay in bone formation, defects in skeletal composition and haematopoiesis.38 We found that PHF20 was up-regulated, suggesting that PHF20 may play a role in bone development through different mechanisms or signal pathways in human osteoporosis. It is noted that PHF20 had a significant diagnostic value for osteoporosis, which may be considered as a biomarker in osteoporosis diagnosis. NFKB2 is reported up-regulated under estrogen treatment in needle bone biopsies of osteoporosis.39 Additionally, it is a key gene in the TRAIL pathway for osteoporosis fractures.40 It is noted that bone resorption marker TRAP5b was undetectable in Nfkb2+/– and Nfkb2–/– mice, and was slightly but significantly increased by TNF injection in Nfkb2+/– mice,41 which suggested the relationship between NFKB2 and bone resorption. Moreover, our study revealed that NFKB2 was significantly enriched in osteoblast differentiation, which further indicated that it was an essential molecule in bone metabolism. Additionally, we found that NFKB2 was remarkably associated with diagnosis and may play a valuable role in the clinical and laboratory diagnosis of osteoporosis.
RPL35A participates in cytoplasmic ribosomal protein pathways in osteoarthritis chondrocytes.42 In addition, it is linked to inherited bone marrow failure syndromes.43 Herein, we found that RPL35A was down-regulated in osteoporosis blood samples compared with normal controls, which provided another pathogenic role in bone disease. BICD2 mutation is involved in spinal muscular atrophy.44-47 Moreover, it has been reported that mutations in BICD2 will cause early onset non-length dependent lower-limb predominant weakness and contractures.48 In the present study, we found decreased expression of BICD2 in osteoporosis, suggesting the vital role in bone formation and metabolism.
There are limitations to our study. First, the sample size in the qRT-PCR data set was small and larger numbers of blood samples of osteoporosis patients are needed for further research. Secondly, the deregulated DEGs in osteoporosis were identified and the definite biological function was not investigated in our study. In vivo and in vitro experiments are essential for elucidation of the biological roles of DEGs in osteoporosis in future work.
In summary, we identified several key genes including VPS35, FCGR2A, TBCA, HIRA, TYROBP, JUND, PHF20, NFKB2, RPL35A and BICD2 involved in the regulation of bone formation and metabolism under the regulation of TFs (FOXD3, Nkx2-5, Pax-4, Oct-1, HNF-4, Pax-6 and COMP1) in osteoporosis. It is noted that VPS35, HIRA, PHF20 and NFKB2 had significant diagnostic value for osteoporosis and may serve as diagnostic biomarkers of osteoporosis. Our results may provide important information for studying the pathogenic mechanisms and consequences of osteoporosis.
Clinical information on the patients in both datasets, a figure showing a heat map image displaying the top 100 differentially expressed genes in osteoporosis compared with relative control groups and a table showing the top seven transcription factors covering the most downstream genes in osteoporosis are available alongside this article online at www.bjr.boneandjoint.org.uk
None declared
None declared
References
1 Nguyen TV , Eisman JA . Genetic profiling and individualized assessment of fracture risk. Nat Rev Endocrinol2013;9:153-161.CrossrefPubMed Google Scholar
2 Rachner TD , Khosla S , Hofbauer LC . Osteoporosis: now and the future. Lancet2011;377:1276-1287.CrossrefPubMed Google Scholar
3 Holroyd C , Cooper C , Dennison E . Epidemiology of osteoporosis. Best Pract Res Clin Endocrinol Metab2008;22:671-685.CrossrefPubMed Google Scholar
4 Chen HL , Deng LL , Li JF . Prevalence of osteoporosis and its associated factors among older men with type 2 diabetes. Int J Endocrinol2013;2013:285729.CrossrefPubMed Google Scholar
5 Wongdee K , Charoenphandhu N . Osteoporosis in diabetes mellitus: possible cellular and molecular mechanisms. World J Diabetes2011;2:41-48.CrossrefPubMed Google Scholar
6 Ralston SH , Uitterlinden AG . Genetics of osteoporosis. Endocr Rev2010;31:629-662.CrossrefPubMed Google Scholar
7 Pietschmann P , Rauner M , Sipos W , Kerschan-Schindl K . Osteoporosis: an age-related and gender-specific disease–a mini-review. Gerontology2009;55:3-12. Google Scholar
8 Seeman E . Bone quality: the material and structural basis of bone strength. J Bone Miner Metab2008;26:1-8.CrossrefPubMed Google Scholar
9 Lampropoulos CE , Papaioannou I , D’Cruz DP . Osteoporosis–a risk factor for cardiovascular disease?Nat Rev Rheumatol2012;8:587-598. Google Scholar
10 Sims AM , Shephard N , Carter K , et al. Genetic analyses in a sample of individuals with high or low BMD shows association with multiple Wnt pathway genes. J Bone Miner Res2008;23:499-506.CrossrefPubMed Google Scholar
11 Huang QY , Recker RR , Deng HW . Searching for osteoporosis genes in the post-genome era: progress and challenges. Osteoporos Int2003;14:701-715.CrossrefPubMed Google Scholar
12 Segal E , Shapira M , Regev A , et al. Module networks: identifying regulatory modules and their condition-specific regulators from gene expression data. Nat Genet2003;34:166-176.CrossrefPubMed Google Scholar
13 No authors listed. The Gene Ontology Consortium. Gene Ontology Consortium: going forward. Nucl Acids Res2015; 43 Database issue D1049–D1056. Google Scholar
14 Kanehisa M . Post-genome Informatics, Oxford University Press, 2000. Google Scholar
15 Kanis JA , McCloskey EV , Johansson H , et al. A reference standard for the description of osteoporosis. Bone2008;42:467-475.CrossrefPubMed Google Scholar
16 Zeng Y , Zhang L , Zhu W. , et al. Quantitative proteomics and integrative network analysis identified novel genes and pathways related to osteoporosis. J Proteomics2016;142:45-52.CrossrefPubMed Google Scholar
17 Parnell LD , Lindenbaum P , Shameer K , et al. 2011 BioStar: An Online Question & Answer Resource for the Bioinformatics Community. PLoS Comput Biol 2011;7:e1002216.CrossrefPubMed Google Scholar
18 Marot G , Foulley JL , Mayer CD , Jaffrezic F . Moderated effect size and P-value combinations for microarray meta-analyses. Bioinformatics2009;25:2692-2699.CrossrefPubMed Google Scholar
19 Tabas-Madrid D , Nogales-Cadenas R , Pascual-Montano A . GeneCodis3: a non-redundant and modular enrichment analysis tool for functional genomics. Nucleic Acids Research2012;40:W478-W483CrossrefPubMed Google Scholar
20 Nogales-Cadenas R , Carmona-Saez P , Vazquez M , et al. GeneCodis: interpreting gene lists through enrichment analysis and integration of diverse biological information. Nucleic Acids Research2009;37:W317-W322.CrossrefPubMed Google Scholar
21 Carmona-Saez P , Chagoyen M , Tirado F , Carazo JM , Pascual-Montano A . GENECODIS: A web-based tool for finding significant concurrent annotations in gene lists. Genome Biology2007;8:R3.CrossrefPubMed Google Scholar
22 ENCODE Project Consortium. The ENCODE (ENCyclopedia Of DNA Elements) Project. Science2004;306:636-640.CrossrefPubMed Google Scholar
23 Robin X , Turck N , Hainard A , et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics2011;12:77.CrossrefPubMed Google Scholar
24 Rivadeneira F , Styrkársdottir U , Estrada K , et al. Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. Nat Genet2009;41:1199-1206.CrossrefPubMed Google Scholar
25 Xia WF , Tang FL , Xiong L , et al. Vps35 loss promotes hyperresorptive osteoclastogenesis and osteoporosis via sustained RANKL signaling. J Cell Biol2013;200:821-837.CrossrefPubMed Google Scholar
26 Xiong L , Xia WF , Tang FL , et al. Retromer in osteoblasts interacts with protein phosphatase 1 regulator subunit 14c, terminates parathyroid hormone’s signaling, and promotes its catabolic response. EBioMedicine. 2016;9:45-60. Google Scholar
27 Raychaudhuri S , Thomson BP , Remmers EF , et al. Genetic variants at CD28, PRDM1 and CD2/CD58 are associated with rheumatoid arthritis risk. Nat Genet2009;41:1313-1318.CrossrefPubMed Google Scholar
28 Tutuncu Z , Kavanaugh A , Zvaifler N , et al. Fcgamma receptor type IIIA polymorphisms influence treatment outcomes in patients with inflammatory arthritis treated with tumor necrosis factor alpha-blocking agents. Arthritis Rheum2005;52:2693-2696.CrossrefPubMed Google Scholar
29 Reveille JD , . Biomarkers for diagnosis, monitoring of progression, and treatment responses in ankylosing spondylitis and axial spondyloarthritis. Clin Rheumatol2015;34:1009-1018.CrossrefPubMed Google Scholar
30 Nolasco S , Bellido J , Gonçalves J , Zabala JC , Soares H . Tubulin cofactor A gene silencing in mammalian cells induces changes in microtubule cytoskeleton, cell cycle arrest and cell death. FEBS Lett2005;579:3515-3524.CrossrefPubMed Google Scholar
31 Kato S , Ishii T , Kouzmenko A . Point mutations in an epigenetic factor lead to multiple types of bone tumors: role of H3.3 histone variant in bone development and disease. Bonekey Rep2015;4:715.CrossrefPubMed Google Scholar
32 Denninger KC , Litman T , Marstrand T , et al. Kinetics of gene expression and bone remodelling in the clinical phase of collagen-induced arthritis. Arthritis Res Ther2015;17:43.CrossrefPubMed Google Scholar
33 McCabe LR , Kockx M , Lian J , Stein J , Stein G . Selective expression of fos- and jun-related genes during osteoblast proliferation and differentiation. Exp Cell Res1995;218:255-262.CrossrefPubMed Google Scholar
34 Wagner EF . Functions of AP1 (Fos/Jun) in bone development. Ann Rheum Dis2002;61(suppl 2):ii40-ii42.CrossrefPubMed Google Scholar
35 McCabe LR , Banerjee C , Kundu R , et al. Developmental expression and activities of specific fos and jun proteins are functionally related to osteoblast maturation: role of Fra-2 and Jun D during differentiation. Endocrinology1996;137:4398-4408.CrossrefPubMed Google Scholar
36 Yeo H , Beck LH , McDonald JM , Zayzafoon M . Cyclosporin A elicits dose-dependent biphasic effects on osteoblast differentiation and bone formation. Bone2007;40:1502-1516.CrossrefPubMed Google Scholar
37 Liu L , Zhu Q , Wang J , et al. Gene expression changes in human mesenchymal stem cells from patients with osteoporosis. Mol Med Rep2015;12:981-987.CrossrefPubMed Google Scholar
38 Badeaux AI , Yang Y , Cardenas K , et al. Loss of the methyl lysine effector protein PHF20 impacts the expression of genes regulated by the lysine acetyltransferase MOF. J Biol Chem2012;287:429-437.CrossrefPubMed Google Scholar
39 Fujita K , Roforth MM , Demaray S , et al. Effects of estrogen on bone mRNA levels of sclerostin and other genes relevant to bone metabolism in postmenopausal women. J Clin Endocrinol Metab2014;99:E81-E88.CrossrefPubMed Google Scholar
40 Zhang YP , Liu YZ , Guo Y , et al. Pathway-based association analyses identified TRAIL pathway for osteoporotic fractures. PLoS One2011;6:e21835.CrossrefPubMed Google Scholar
41 Yao Z , Xing L , Boyce BF . NF-kappaB p100 limits TNF-induced bone resorption in mice by a TRAF3-dependent mechanism. J Clin Invest2009;119:3024-3034.CrossrefPubMed Google Scholar
42 Tsolis KC , Bei ES , Papathanasiou I , et al. Comparative proteomic analysis of hypertrophic chondrocytes in osteoarthritis. Clin Proteomics2015;12:12.CrossrefPubMed Google Scholar
43 Chung NG , Kim M . Current insights into inherited bone marrow failure syndromes. Korean J Pediatr2014;57:337-344.CrossrefPubMed Google Scholar
44 Neveling K , Martinez-Carrera LA , Hölker I , et al. Mutations in BICD2, which encodes a golgin and important motor adaptor, cause congenital autosomal-dominant spinal muscular atrophy. Am J Hum Genet2013;92:946-954.CrossrefPubMed Google Scholar
45 Oates EC , Rossor AM , Hafezparast M , et al. Mutations in BICD2 cause dominant congenital spinal muscular atrophy and hereditary spastic paraplegia. Am J Hum Genet2013;92:965-973.CrossrefPubMed Google Scholar
46 Peeters K , Litvinenko I , Asselbergh B , et al. Molecular defects in the motor adaptor BICD2 cause proximal spinal muscular atrophy with autosomal-dominant inheritance. Am J Hum Genet2013;92:955-964.CrossrefPubMed Google Scholar
47 Synofzik M , Martinez-Carrera LA , Lindig T , Schöls L , Wirth B . Dominant spinal muscular atrophy due to BICD2: a novel mutation refines the phenotype. J Neurol Neurosurg Psychiatry2014;85:590-592.CrossrefPubMed Google Scholar
48 Rossor AM , Oates EC , Salter HK , et al. Phenotypic and molecular insights into spinal muscular atrophy due to mutations in BICD2. Brain2015;138:293-310.CrossrefPubMed Google Scholar
