Objective: To design a weighted co-expression network and build gene expression signature-based nomogram (GESBN) models for predicting the likelihood of bone metastasis in breast cancer (BC) patients. Methods: Dataset GSE124647 was used as a training set, and GSE14020 was taken as a validation set. In the training cohort, limma package in R was adopted to obtain differentially expressed genes (DEGs) between BC non-bone metastasis and bone metastasis patients, which were used for functional enrichment analysis. After weighted co-expression network analysis (WGCNA), univariate Cox regression and Kaplan-Meier plotter analyses were performed to screen potential prognosis-related genes. Then, GESBN models were constructed and evaluated. Further, the expression levels of genes in the models were explored in the training set, which was validated in GSE14020. Finally, the prognostic value of hub genes in BC was explored. Results: A total of 1858 DEGs were obtained. WGCNA result showed that the blue module was most significantly related to bone metastasis and prognosis. After survival analyses, GAJ1, SLC24A3, ITGBL1, and SLC44A1 were subjected to construct a GESBN model for overall survival. While GJA1, IGFBP6, MDFI, ITGFBI, ANXA2, and SLC24A3 were subjected to build a GESBN model for progression-free survival. Kaplan-Meier plotter and receiver operating characteristic analyses presented the reliable prediction ability of the models. Besides, GJA1, IGFBP6, ITGBL1, SLC44A1, and TGFBI expressions were significantly different between the two groups in GSE124647 and GSE14020. The hub genes had a significant impact on patient prognosis. Conclusion: Both the four-gene signature and six-gene signature could accurately predict patient prognosis, which may provide novel treatment insights for BC bone metastasis.