Abstract Osteoporosis (OP) is a metabolic bone disease characterized by reduced bone density and fragility, impairing quality of life. Traditional treatments often overlook symptoms like back and joint pain, increasing burden. This study aims to map relationships between natural medicines, targets, and symptom clusters, demonstrating their effectiveness in personalized OP treatment to enhance clinical strategies and self-assessment. We used compounds and targets, applying Summary data-based Mendelian Randomisation (SMR) analysis for biological process and molecular function enrichment. Additionally, we employed Phenome-Wide Association Studies (PheWAS) to select two natural drugs—Rhizoma Drynariae (RD) and Lycii Fructus (LF)—for case analysis. The study found that RD primarily improves symptoms such as indigestion, constipation, fatigue, polyuria, and depression, while LF significantly ameliorates symptoms related to the nervous and muscular systems, such as hoarseness, dizziness, vertigo, and fever symptoms. This analysis successfully differentiated two groups of symptoms and precisely constructed a logical chain among “natural Medicines-molecular tArGets-Illness-symptom Clusters” (MAGIC chain) achieving a refined classification of OP. The results of this study support the effectiveness of implementing personalized medical strategies in the treatment of OP, providing a scientific basis for the clinical application of natural medicines and patient self-management. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-95304-3. Keywords: Natural medicines, Molecular targets, Illness, Symptom clusters, Osteoporosis Subject terms: Computational biology and bioinformatics, Risk factors, Signs and symptoms Introduction Osteoporosis (OP) is a widespread and serious global health issue characterized by significantly reduced bone density and deteriorating bone structure, consequently elevating the risk of fractures^[46]1,[47]2. According to data from the World Health Organization, OP affects millions worldwide, particularly postmenopausal women^[48]3,[49]4. As the global population ages, the prevalence of this condition is expected to increase further, posing substantial challenges to individual health, medical systems, and socio-economic stability^[50]5. Currently, the treatment strategies for OP primarily include lifestyle modifications, nutritional supplementation, and pharmacotherapy^[51]6. Common medications such as bisphosphonates, calcitonin, selective estrogen receptor modulators (SERMs), parathyroid hormone (PTH) analogs, and recent RANKL inhibitors (such as denosumab), can effectively increase bone density and reduce the risk of fractures^[52]7,[53]8. However, long-term use of these drugs may lead to a series of side effects, such as gastrointestinal discomfort, cardiovascular issues, and even an increased risk of certain types of cancer, leading many patients to be cautious about these treatment methods^[54]9–[55]11. Additionally, current pharmacological treatment principles focus only on slowing bone loss, overlooking many clinically significant symptoms such as pain and decreased mobility, which have a profound impact on patients’ quality of life^[56]11. Therefore, an increasing clinical demand indicates that, in addition to improving bone density and slowing bone loss, attention to patients’ own experiences, such as enhancing quality of life and alleviating pain, has also become an urgent need in the clinical treatment of osteoporosis^[57]12. In the face of limitations in contemporary Western medicine for the treatment of OP, Natural medicines, represented by Chinese herbal medicine, offer a unique and multidimensional treatment option, particularly suitable for patients seeking to reduce medication side effects and enhance physiological functions^[58]13,[59]14. For instance, based on the clinical guidelines, medications such as Rhizoma Drynariae (RD) and Lycii Fructus (LF) are not only widely utilized but their therapeutic efficacy has been rigorously validated^[60]15. The therapeutic actions of these herbal medicines extend beyond merely slowing bone loss; current research also indicates that they can significantly alleviate various symptoms associated with osteoporosis, such as pain, fatigue, and muscle spasms^[61]16. However, although natural medicines show significant therapeutic potential and advantages, the effects of these medicines and their active components on specific disease targets, as well as the biological connections between these components and clinical symptoms, have not been fully elucidated. Therefore, there is an urgent need for more scientific research to deeply explore and verify these potential medical correlations. Therefore, this study aims to preliminarily construct a logical chain among “natural Medicines-molecular tArGets-Illness-symptom Clusters” (MAGIC chain) achieving a refined classification of OP from a symptomatic perspective and establishing specific associations with corresponding medications. This approach provides patients with more accurate self-assessment and early warning guidance. We have selected two natural medicines—RD and LF—as our case studies, screening and summarizing the active components and their targets in these herbs. Using Summary data-based Mendelian Randomisation (SMR) analysis, we identified the targets of these drugs in the treatment of OP and further explored the associated biological processes, molecular functions, and signaling pathways to elucidate the underlying mechanisms of their therapeutic effects. By utilizing Phenome-Wide Association Studies (PheWAS) and the SymMap tool, we further constructed mappings between symptoms and corresponding targets, thereby precisely categorizing different treatment groups within natural medicines. This study, under the premise of conventional OP treatment, aims to provide a basis for the selection of natural medicines in clinical treatment by defining different symptom clusters, guiding clinical medication choices, and assisting OP patients in effective self-diagnosis and management (Fig. [62]1). Fig. 1. [63]Fig. 1 [64]Open in a new tab Precision treatment pathway from natural medicines to molecular targets, illness, and symptom clusters. This figure illustrates the comprehensive methodology employed to evaluate the effectiveness of Rhizoma Drynariae and Lycii Fructus in the treatment of osteoporosis through a detailed network of drug ingredients, molecular targets, and symptom clusters. Materials and methods Study design We selected compounds with high oral bioavailability (OB) and drug-likeness (DL) from the Traditional Chinese Medicine Systems Pharmacology Database (TCMSP), setting our thresholds at OB ≥ 30% and DL ≥ 0.18. We then used the PubChem database, Swiss Target Prediction, and Target Net to identify potential targets for these compounds. To pinpoint specific therapeutic targets for OP, we conducted SMR analysis, followed by Colocalization analysis to enhance the robustness of our causal inferences. Furthermore, KEGG pathway and Gene Ontology (GO) enrichment analyses were performed on the identified targets to systematically explore the associated biological processes, molecular functions, and signaling pathways. These enrichment analyses provided comprehensive insights into the molecular mechanisms underlying the therapeutic effects of the selected compounds on OP. For validated targets, we utilize a combination of PheWAS analysis and the SymMap database to screen for target-symptom relationships related to OP and clarify their clinical features. Simultaneously, by analyzing the tissue-specific differential gene expression of RD and LF using the Functional Mapping and Annotation (FUMA) platform, we can determine the effects of these medicines in different tissues, thereby accurately linking them to specific OP symptoms. Data source The expression quantitative trait loci (eQTL) data, encompassing 31,684 individuals, were sourced from the eQTLGen Consortium^[65]17. The tissue-specific expression data of target genes, indicative of a potential causal effect on OP, were obtained from tissue-specific eQTL datasets via the Genotype-Tissue Expression (GTEx) portal ([66]https://gtexportal.org/home/)^[67]18. The GTEx v8 dataset comprises data from 838 donors, including 17,382 samples across 52 tissues and two cell lines. To assess the impact of drug target action on the risk of OP across different tissues, we utilized eQTLs data from tissues where the target genes are highly expressed, publicly available in the GTEx. SNPs were selected as genetic variants, with a false discovery rate-adjusted P-value < 0.05, and further clustering was applied for linkage disequilibrium (LD) at r^2 < 0.20. All GWAS summary data utilized in this study were derived from individuals of European ancestry. Osteoporosis was diagnosed based on the World Health Organization (WHO) criteria, which defines osteoporosis as a bone mineral density (BMD) T-score ≤ -2.5 at the lumbar spine, femoral neck, or total hip, measured by dual-energy X-ray absorptiometry (DXA). The GWAS summary statistics for OP were obtained from the Finngen database ([68]https://www.finngen.fi/en), encompassing genetic factors from 3,203 patients diagnosed with OP and 209,575 control participants^[69]19. Identification of drug target genes Using the TCMSP database ([70]https://tcmsp-e.com/tcmsp.php), we identified compounds with high OB and DL^[71]20. OB is a metric that assesses the rate and extent of a drug being absorbed into the systemic circulation, reflecting the proportion of the drug entering the systemic circulation. DL measures the similarity of a compound to known drugs, which is crucial for optimizing pharmacokinetics and drug properties. We set the screening criteria at OB ≥ 30% and DL ≥ 0.18 to select compounds in the TCMSP database that exhibited higher activity^[72]21,[73]22. To retrieve the isomeric SMILES identifiers for candidate compounds, we searched the PubChem database using each compound’s PubChem Cid. If a PubChem Cid was not available, we downloaded the molecular image from the TCMSP database and used the Virtual Computational Chemistry Laboratory’s structure calculation site ([74]http://www.vcclab.org/web/alogps/) to compute the SMILES identifier. After obtaining the SMILES identifiers for the active metabolites, we imported them into the SwissTargetPrediction ([75]https://www.swisstargetprediction.ch/) and TargetNet ([76]http://targetnet.scbdd.com/) databases to predict the potential targets of these compounds^[77]23,[78]24. Summary data-based Mendelian randomization SMR was used to evaluate the association between methylation, gene expression, and protease protein levels of drug targets and the risk of developing OP. Compared to conventional MR, SMR provides greater statistical power, particularly when leveraging exposures and outcomes from two large, independent cohorts and focusing on top cis-QTLs. We selected a leading cis-QTL within a ± 1000 kb window centered on each gene of interest, applying a genome-wide significance threshold of P-value < 5.0 × 10⁻⁸. This stringent cutoff helps minimize false positives while ensuring that the variants included are highly relevant to gene-level regulation. SNPs demonstrating allele frequency disparities greater than a predetermined threshold (0.2 for this study) across any dataset pairs (including LD reference samples, QTL, and outcome summary data) were excluded. The choice of an LD r^2 threshold of 0.2 was based on the balance between ensuring sufficient independent instrumental variables while minimizing bias from linkage disequilibrium. This threshold prevents the inclusion of highly correlated variants that could confound causal inference while allowing enough SNPs to capture gene-associated variations. The homogeneity in dependent instruments (HEIDI) test was utilized to discern pleiotropy from linkage, with a P-HEIDI value < 0.01 indicating pleiotropy, thereby excluding such instances from the analysis. Both SMR and HEIDI tests were conducted using the SMR software (version 1.3.2). Colocalization analysis In our Colocalization analysis, we report five different posterior probabilities corresponding to five exclusive hypotheses: (1) No causal variant for any trait (H0); (2) Causal variant only for gene expression; (3) Causal variant only for disease risk (H2); (4) Different causal variants for the two traits (H3); (5) The same shared causal variant for the two traits (H4). For drug target-related Colocalization, the Colocalization region window is set to ± 1000 kb based on published articles. The prior probabilities of a causal variant being associated only with Trait 1 (i.e., Target) and only with Trait 2 (i.e., OP) are set at 10^−4. A posterior probability for H4 (PPH4) greater than 0.70 is considered to provide supportive evidence for Colocalization and its significance threshold, corresponding to a false discovery rate of < 0.05, thus strengthening the evidence of causality. Annotation and analysis of target genes via FUMA In our study, we utilized target genes obtained through the SMR method and further analyzed and annotated the results of genome-wide association studies (GWAS) using the FUMA platform ([79]https://fuma.ctglab.nl/). By employing the GENE2FUNC tool and tissue-specific data provided by GTEx V8, we conducted in-depth annotations of these priority genes, aiming to explore their presumed biological mechanisms and expression patterns across various tissues. Enrichment analysis To elucidate the biological functions and relevance of the potential therapeutic targets identified in this study, we performed comprehensive enrichment analyses using both GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) frameworks ([80]www.kegg.jp/kegg/kegg1.html)^[81]25–[82]27. The GO analysis was segmented into its three conventional categories: biological process (BP), molecular function (MF), and cellular component (CC). In parallel, KEGG pathway analysis was employed to pinpoint key pathways enriched with the target proteins. Both analyses were carried out with the R packages ClusterProfiler and Pathview. Enrichment was determined using a hypergeometric test, applying a significance threshold of P < 0.05. Molecular Docking To gain atomic-level insights into the interactions between candidate compounds and their respective target proteins, molecular docking simulations were conducted using AutoDock Vina v.1.2.2 ([83]http://autodock.scripps.edu/)^[84]28. This approach enabled the assessment of ligand–receptor binding affinities and interaction modes, which are critical for prioritizing drug targets for further experimental validation and optimizing potential drug candidates. For each selected compound, docking was performed against the top two target proteins; when only one target was available, docking was conducted solely with that protein. Drug structures were obtained from the PubChem Compound database ([85]https://pubchem.ncbi.nlm.nih.gov/)^[86]29. Protein structural data were sourced from the AlphaFold Protein Structure Database ([87]https://alphafold.ebi.ac.uk/). Prior to docking, water molecules were removed from all protein and ligand files, and polar hydrogen atoms were added. The grid box was centered on each protein’s structural domain to allow unrestricted molecular movement, using a grid point spacing of 0.05 nm and dimensions of 30 Å × 30 Å × 30 Å. The entire docking process was visualized using AutoDock Vina v.1.2.2. Phenome-wide Mendelian randomization To evaluate the potential for horizontal pleiotropy and any unintended effects of the proposed drug targets, we undertook a PheWAS analysis via the AstraZeneca portal ([88]https://azPheWAS.com/)^[89]30,[90]31. This analysis was based on data from the UK Biobank, which included approximately 15,500 binary phenotypes and 1,500 continuous phenotypes derived from around 450,000 participants who had undergone exome sequencing. To mitigate the risk of false positives, we applied multiple corrections and established a significance threshold of 2 × 10^−9. The detailed methodology is documented in the original publication. Finally, by using SymMap ([91]https://www.symmap.org/), the screened targets are correlated with clinical symptoms recognized in modern medicine^[92]32. Data visualization All figures were generated using R software (version 4.3.2). The primary visualization was performed using the ggplot2 package for data plotting. Additional R packages, including ComplexHeatmap for heatmap visualization, ggpubr for statistical annotation, pheatmap for hierarchical clustering heatmaps, and cowplot for figure arrangement, were utilized as needed. Custom scripts were used for further refinement of graphical outputs. Results Acquisition of active ingredients and target prediction of RD and LF We successfully identified key active components in RD and LF using the TCMSP. The primary active components of RD include Aureusidin, Beta-sitosterol, and Digallate, while those in LF include 31-Norcycloartanol, Stigmasterol, and Fucosterol. We further used the PubChem database to retrieve the SMILES identifiers for these components and predicted their potential targets using the SwissTargetPrediction and TargetNet databases. The results of these target predictions were organized along with the active components to preliminarily construct an association network between the active components of RD and LF and their potential drug targets (Tables S1, S2). Molecular targets of RD and LF in treatment of OP In this study, we compiled and organized the components and their targets from RD and LF and conducted systematic SMR analyses to explore their potential effects on OP (Fig. [93]2). The analysis revealed common targets as well as unique targets within these two herbs. Among the common targets, the ESR1 gene was found to be significantly associated with multiple compounds in both herbs. For instance, (R) − naringenin from RD and 24 − ethylcholest − 22 − enol from LF both demonstrated a significant positive modulation of the ESR1 gene (β = 0.429, P = 0.001, OR = 1.536, 95% CI: 1.180-2.000), suggesting that these components might collaboratively participate in the pathophysiology of OP through regulation of the ESR1 gene. Regarding unique targets, RD specific targets included VEGFA and ABCG2, with Aureusidin significantly affecting VEGFA (β = 0.202, P = 0.004, OR = 1.224, 95% CI: 1.066–1.405) and beta-sitosterol impacting ABCG2 (β = 0.153, P = 0.032, OR = 1.165, 95% CI: 1.013–1.340). In LF, 24 − ethylcholesta − 5,22 − dienol significantly affected HMGCR (β = 0.285, P = 0.036, OR = 1.33, 95% CI: 1.02–1.74), and Sitosterol alpha1 had a notable impact on ADRA2A (β = 0.687, P = 0.045, OR = 1.99, 95% CI: 1.02–3.89). Fig. 2. [94]Fig. 2 [95]Fig. 2 [96]Open in a new tab Molecular targets of RD and LF in treatment of OP. (A) Forest plot of RD compound efficacy on target genes associated with osteoporosis. (B) Forest plot of LF compound efficacy on target genes associated with osteoporosis. Colocalization analysis results for shared genetic targets In the colocalization analysis of the gene PARP1, which is a common target between RD and LF, the posterior probability (PPH4) reached 0.87 (Fig. [97]3). This high value suggests strong evidence for the presence of a shared causal variant that influences both the expression of PARP1 and the susceptibility to osteoporosis. This finding indicates a significant overlap in the genetic determinants of PARP1 expression and osteoporosis risk, as highlighted by the shared association in both natural medicines. Fig. 3. [98]Fig. 3 [99]Open in a new tab Colocalization analysis of PARP1 for shared genetic influence on osteoporosis. This figure displays the colocalization analysis results for the PARP1 gene, identified by its Ensembl gene ID ENSG00000143799. Tissue-specific gene expression of RD and LF in treatment of OP In this study, we analyzed the expression of target genes of RD and LF in the treatment of OP across various tissues using the GTEx database. This systematic analysis revealed the specific tissue expression of these herbal targets, elucidating their potential activity in different tissues and their association with specific symptoms. For RD, the expression of target genes is primarily concentrated in tissues related to the digestive system, such as the stomach, intestines, pancreas, and liver. This correlates closely with RD’s specific treatment symptoms such as dyspepsia and constipation. Moreover, the normal functioning of the digestive system is crucial for the absorption of calcium and other minerals, which is especially important for the bone health of osteoporosis patients. For LF, the expression of target genes is significant in the brain, heart, and muscles. This relates to LF’s unique symptoms such as headaches, dizziness, cardiac symptoms (such as angina), and muscle weakness (Fig. [100]4). Fig. 4. [101]Fig. 4 [102]Open in a new tab Tissue-specific gene expression of RD and LF in treatment of OP. (A) Represents the tissue-specific analysis of target sites for active components in RD; (B) Represents the tissue-specific analysis of target sites for active components in LF, with the color gradient from blue (low expression) to red (high expression) indicating the expression levels of targets across different tissues. Comparative pathway and functional enrichment analysis of RD and LF This study conducted KEGG pathway and GO functional category enrichment analyses to explore the molecular-level differences in the actions of RD and LF (Fig. [103]5A and B). The KEGG pathway analysis revealed that RD primarily affects cell cycle regulation, the NF-κB signaling pathway, the VEGF signaling pathway, and the PI3K-Akt signaling pathway. The significant enrichment in these pathways suggests that RD may play a key role in cell growth, metabolic regulation, and signal transduction. In contrast, the pathways affected by LF are markedly different, mainly involving cell apoptosis, cortisol synthesis and secretion, fat digestion and absorption, and base excision repair. Particularly noteworthy is the significant enrichment in the cell apoptosis pathway, highlighting LF’s potential role in regulating cell death and maintaining cellular homeostasis. Fig. 5. [104]Fig. 5 [105]Open in a new tab KEGG pathway and GO enrichment analysis of RD and LF. (A, B) KEGG pathway enrichment analysis of RD (A) and LF (B), showing significantly enriched pathways. (C, D) GO enrichment analysis results of RD (C) and LF (D), categorized into biological processes (BP), cellular components (CC), and molecular functions (MF). The GO enrichment analysis (Fig. [106]5C and D) categorized significant terms into biological process (BP), cellular component (CC), and molecular function (MF). For RD (C), key BP terms included vascular processes in the circulatory system and hormone metabolic processes, while major CC terms involved extracellular region components. Enriched MF terms, such as hormone receptor binding and kinase activity, further indicate the importance of signaling and metabolic regulation. For LF (D), major BP terms included regulation of protein-containing complex assembly and mammary gland alveolus development, implicating both protein metabolism and tissue-specific processes. Meanwhile, the CC category was enriched for membrane-bound organelle components, and MF terms highlighted DNA binding and hormone receptor activity. Molecular Docking results A molecular docking analysis was conducted to investigate the binding interactions between key targets and constituents derived from LF and RD (Fig. [107]6; Table [108]1). Overall, negative binding energies were observed for all docked complexes, indicating thermodynamically favorable binding. Notably, RD-derived aureusidin exhibited the strongest binding affinity to PARP1 (− 9.036 kcal/mol), followed by luteolin (− 8.816 kcal/mol) and kaempferol (− 8.387 kcal/mol). These results suggest that multiple flavonoid compounds, particularly those from RD, have a high potential for PARP1 inhibition. In addition, RD-derived naringenin demonstrated notable binding to both VEGFA (− 5.88 kcal/mol) and YWHAG (− 6.682 kcal/mol), implying possible multifunctional modulation. For LF, atropine bound to BAZ2A with a binding energy of − 5.82 kcal/mol, while ethyl linolenate showed weaker affinities for SCARB1 (− 4.292 kcal/mol) and PARP1 (− 5.444 kcal/mol) relative to other docked ligands. Quercetin from LF also displayed significant affinity for PARP1 (− 8.1 kcal/mol), further highlighting the potential of polyphenolic compounds in modulating this target. Fig. 6. [109]Fig. 6 [110]Open in a new tab Molecular docking of RD and LF active compounds with their respective and common targets. Molecular docking results of key active compounds from LF and RD with their predicted protein targets. (A) Docking of atropine with BAZ2A, a potential target of LF. (B) Docking of ethyl linolenate with SCARB1, another potential LF target. (C) Docking of naringenin with VEGFA, a target of RD. (D) Docking of naringenin with YWHAG, another RD-associated target. (E–I) Docking results of the common target PARP1 with different active compounds from RD and LF, including: (E) Quercetin, (F) Aureusidin, (G) Ethyl linolenate, (H) Kaempferol, and (I) Luteolin. Table 1. Molecular Docking results for common and individual targets of RD and LF. Target Drug Identification Binding energy (kcal/mol) BAZ2A Lycii Fructus atropine -5.82 SCARB1 Lycii Fructus Ethyl linolenate -4.292 VEGFA Rhizoma Drynariae narngenin -5.88 YWHAG Rhizoma Drynariae narngenin -6.682 PARP1 Lycii Fructus Quercetin -8.1 PARP1 Rhizoma Drynariae Aureusidin -9.036 PARP1 Lycii Fructus Ethyl linolenate -5.444 PARP1 Rhizoma Drynariae Kaempferol -8.387 PARP1 Rhizoma Drynariae Luteolin -8.816 [111]Open in a new tab Symptomatic analysis of RD and LF based on specific genetic targets By analyzing the common targets of RD and LF, we identified a range of symptoms potentially influenced by both herbs. Notably, these symptoms include several clinical manifestations closely related to OP, such as back pain, neck pain, joint pain, muscle spasms, and muscle weakness. An analysis of the symptoms associated with specific genes unique to RD revealed that, compared to the common genetic targets, the unique gene expression of RD affects a distinct set of clinical presentations. These symptoms primarily include dyspepsia, constipation, fatigue, polyuria, and depression, which are collectively referred to as “kidney‑yang deficiency syndrome.” In analyzing the symptoms related to specific genes in LF, we identified a broad spectrum of clinical manifestations affecting various body systems, including acute pain, dizziness, headache, fever symptoms, skin itching, chronic fatigue, and muscle weakness, commonly classified under “kidney‑yin deficiency syndrome.” Discussion Study summary Constructing associations among medications, diseases, and symptoms is crucial for understanding a range of severe symptoms caused by osteoporosis, such as pain, muscle weakness, and frailty. However, the current prevention and treatment of OP mainly adopt a single-disease targeted treatment model. For this reason, we used the active components of medications and their targets, employing SMR and PheWAS analysis methods, to build a “natural medicines-molecular-targets-illness-symptom clusters " association network. Guided by this framework, we selected two traditional medicines—RD and LF—as our case studies, both of which have been extensively proven to be clinically effective. After screening and summarizing the active components and their targets in these natural medicines, we further confirmed their effectiveness in treating OP through SMR analysis. Additionally, using the FUMA platform, we determined the tissue specificity of these medications in treating OP. The application of Phewas and SymMap revealed that the differential gene targeting by these two medicines addresses various types of OP symptoms, thus breaking away from the traditional single treatment model and providing a new path for symptom-based personalized treatment. Enhancing osteoporosis diagnosis and treatment through symptom-based classification OP is a serious disease that is often overlooked, with early symptoms such as pain, muscle spasms, and decreased balance frequently not receiving the attention they deserve. This often leads to many patients being diagnosed only after experiencing an osteoporotic fracture, causing irreversible damage. Although previous studies have explored the clinical features and correlations of OP, they often overlook the importance of these common symptoms, which are crucial for the timely detection and treatment of OP^[112]33. This study not only reveals the therapeutic effects of different medications on OP but also further identifies how the specific action targets of these drugs map to two distinct symptom clusters. This symptom-based classification method effectively reveals the subjective sensory conditions of OP patients, addressing an important aspect previously neglected in research. Through detailed study of these symptom clusters, we not only deepen our understanding of the disease’s diverse manifestations but also enhance our ability to detect OP early. Furthermore, the results of this study will optimize the screening, early intervention, and comprehensive management of OP, providing patients with more precise and personalized medical services. Action mechanisms of RD and LF RD and LF are well-known natural medicines, embodying rich concepts of complementary and alternative medicine.They are also key ingredients in commercially available drugs such as Qianggu Capsules and are recommended by clinical practice guidelines as first-class medications for OP^[113]34,[114]35. Therefore, we have chosen these two herbs as representatives in our study, aiming to construct a “natural medicines-molecular targets-illness-symptom clusters” association network. In this study, KEGG and GO enrichment analyses revealed that RD and LF are involved in distinct molecular mechanisms and functional pathways related to osteogenesis and bone homeostasis, indicating their differential potential in preventing and treating osteoporosis. From the perspective of RD, its significant enrichment in the PI3K–Akt signaling pathway highlights its role in promoting cell survival, proliferation, and differentiation, thereby directly enhancing osteoblast activity and facilitating bone tissue regeneration^[115]36,[116]37. In addition, RD showed direct enrichment of the NF-κB signaling pathway, a key mediator of inflammatory responses and bone resorption, suggesting a robust anti-inflammatory effect that could help reduce bone loss^[117]38. The enrichment of pathways related to cell cycle regulation and angiogenesis further implies that RD may improve local blood supply and promote the recruitment of reparative factors, accelerating the bone repair process^[118]39. Complementing these findings, molecular docking analysis provided additional support: key active compounds from RD, such as aureusidin and naringenin, demonstrated favorable binding affinities with critical targets including PARP1, VEGFA, and YWHAG, suggesting that these compounds might contribute to modulating cell signaling, inflammation, and angiogenesis. In contrast to RD, LF was predominantly enriched in pathways related to cell apoptosis, cortisol synthesis and secretion, fat digestion and absorption, and base excision repair, suggesting distinct mechanisms for bone protection. The notable enrichment in the apoptosis pathway aligns with prior findings that controlled cell death is crucial for maintaining bone integrity^[119]40. Meanwhile, LF’s impact on cortisol-related pathways points to a possible modulation of glucocorticoid levels, known to influence bone turnover and osteoporosis risk^[120]41. The involvement in fat digestion and absorption highlights LF’s potential role in optimizing essential micronutrient uptake, which is critical for bone metabolism^[121]42. Furthermore, the enrichment of the base excision repair pathway suggests that LF may bolster DNA repair and mitigate oxidative stress in osteoblasts. Molecular docking analyses support this view, as LF-derived compounds (e.g., quercetin, atropine) showed binding to DNA repair targets such as PARP1, reinforcing LF’s potential to maintain cellular homeostasis and bone health. Comparative analysis of key active compounds in RD and LF RD has shown significant effects on symptoms such as edema, reduced urine output, cold intolerance, lower back pain, and fatigue in patients with OP. These clinical manifestations are closely linked to specific active compounds in RD, including flavonoids. These components have been proven to regulate bone metabolism and enhance bone reconstruction^[122]43,[123]44. Research indicates that RD promotes the proliferation of osteoblasts and inhibits the activity of osteoclasts, thereby maintaining bone balance, effectively enhancing bone density, and improving bone structure^[124]39,[125]45,[126]46. By improving blood circulation and enhancing vascular health, these antioxidants also help alleviate symptoms of edema and cold intolerance caused by poor circulation^[127]47. Additionally, these compounds may help relieve fatigue and weakness by regulating neurotransmitters, thus enhancing the overall vitality and quality of life of patients^[128]48. LF has demonstrated significant improvements in symptoms such as lower back pain, emaciation, insomnia, irritability, and anxiety in patients with OP. These clinical manifestations are closely associated with active components found in LF, particularly polysaccharides, β-carotene, and vitamin C^[129]49. These components have been proven to regulate neurotransmitters and hormone levels, which are crucial for maintaining bone health and neurological stability. Research has shown that the polysaccharides in LF can improve bone density, thereby directly treating the pain caused by OP^[130]50. Furthermore, β-carotene and vitamin C, as potent antioxidants, help reduce oxidative stress, which is often associated with the emaciation symptoms linked to OP^[131]51. Strengths and limitations In this study, we employed a highly selective screening method, ensuring that only compounds with high activity and excellent drug-like properties were selected for in-depth research. This precise screening strategy significantly enhances the success rate of subsequent experiments and their potential in clinical applications. Additionally, our research design included comprehensive target validation and clinical feature correlation analysis. This not only helped us identify specific targets of the compounds but also allowed for a deeper exploration of the relationship between these targets and clinical features. By using the Mendelian randomization approach, utilizing genetic variants as instrumental variables, we could more reliably infer causal relationships. Compared to traditional observational studies, this method more directly demonstrates the true impact of the active ingredients in RD on specific physiological processes, rather than merely superficial associations. Our study has several limitations. Firstly, although MR provides insights into causality, it assumes a linear relationship between exposure and outcome, which may not capture non-linear or U-shaped exposure-response relationships. Secondly, the QTL used in our study may only show minor differences in gene expression levels, which may not fully capture the potential impact of genes. Additionally, this study primarily targeted a European population, limiting the general applicability of the findings to other ethnic groups. Therefore, further research and validation are necessary to generalize the results to other ethnicities. Although this study identified several potentially active compounds through computational methods, these results need to be further validated experimentally to determine their actual effects in biological systems. Future prospects Based on OP, the research strategy of classifying patients according to different symptom types offers a profound new pathway for the treatment with natural medicines. This approach breaks away from the traditional single-treatment model for OP, allowing us to provide more precise treatment plans based on the specific symptoms exhibited by patients. By identifying biomarkers and potential therapeutic targets associated with these symptoms, the complex components of natural medicines such as RD and LF can be utilized more effectively to alleviate symptoms and fundamentally treat the disease. Future research should focus on developing and validating these symptom classification standards, while also exploring how the components of natural medicines interact with these classifications to optimize treatment effects. This will require interdisciplinary collaboration, involving experts from the fields of complementary and alternative medicine, biostatistics, molecular biology, and clinical medicine. Additionally, by conducting large-scale clinical trials and longitudinal studies, not only can the clinical relevance of these classifications be verified, but the comprehensive effects of natural medicines in treating multiple symptoms can also be better understood. The successful implementation of this comprehensive treatment strategy will provide OP patients with more holistic and personalized treatment options. Conclusion In this study, we successfully established a " natural medicines-molecular-targets-illness-symptom clusters " network model for OP using RD and LF, paving new pathways for disease management and drug development. By clearly distinguishing two unique symptom clusters associated with OP, we not only revealed symptoms that are often overlooked in clinical settings but also provided new strategies for personalized treatment of OP. Additionally, these symptom clusters offer a practical self-monitoring tool for OP patients, aiding them in effective symptom management and self-assessment in daily life. Our research opens new vistas for clinical treatment and drug development for OP, demonstrating the importance of integrating drug, pathological, and symptom data. It provides a scientific foundation and practical guidance for precision medicine and patient self-management in the future of OP treatment. Electronic supplementary material Below is the link to the electronic supplementary material. [132]Supplementary Material 1^ (64.1KB, xlsx) [133]Supplementary Material 2^ (60KB, xlsx) [134]Supplementary Material 3^ (19KB, xlsx) [135]Supplementary Material 4^ (46.5KB, xlsx) [136]Supplementary Material 5^ (151.5KB, xlsx) [137]Supplementary Material 6^ (85.8KB, xlsx) Author contributions X.W., Y.L.Z. and L.G.Z. contributed to the conception and design of the study. X.Y.G., J.R.Q. and K.S. contributed to the literature search and data extraction. Y.L., Q.Q.L. and Z.K.J. contributed to data analysis and interpretation. S.Y.H., Y.W.G. and L.W. contributed to the software. X.Y.G., K.S. and Y.L.Z. wrote the first draft of the manuscript. X.W. and L.G.Z. review and edited the manuscript. All authors contributed to critical revision of the manuscript. All authors approved the manuscript. Funding The study was supported by the National Natural Science Foundation of China (grant no. 82205140), the Basic Research Program of Jiangsu Province (Natural Science Foundation; grant no. BK20220468), the Innovation Team and Talents Cultivation Program of National Administration of Traditional Chinese Medicine (grant no. ZYYCXTD-C-202003), the Fundamental Research Funds for the Central public welfare research institutes (grant no. ZZ13-YQ-039) and the self-selected topic project for Wangjing Hospital of China Academy of Chinese Medical Sciences (grant no. WJYY-ZZXT-2023-16). Data availability The datasets generated and/or analyzed during the current study are available from publicly accessible repositories. Specifically: GWAS summary statistics were obtained from the FinnGen study ([138]https://www.finngen.fi/en) and the GWAS Catalog ([139]https://www.ebi.ac.uk/gwas/), accessed on August 15, 2024.eQTL data were sourced from the Genotype-Tissue Expression (GTEx) project v8 ([140]https://www.gtexportal.org/home/), accessed on August 10, 2024, and the eQTLGen Consortium ([141]https://www.eqtlgen.org/), accessed on July 25, 2024.Traditional Chinese Medicine-related data were retrieved from the Traditional Chinese Medicine Systems Pharmacology Database (TCMSP) ([142]https://old.tcmsp-e.com/tcmsp.php), accessed on June 30, 2024, and the SymMap database ([143]https://www.symmap.org/), accessed on July 5, 2024.Phenome-wide association study (PheWAS) ([144]https://atlas.ctglab.nl/) analysis utilized publicly available data, accessed on July 20, 2024. Declarations Competing interests The authors declare no competing interests. Ethics approval and consent to participate This study involves data obtained from a public database, which has been subjected to secondary analysis under the authorization of data use permissions. Given that the study utilizes anonymized data, in accordance with relevant ethical guidelines and data use agreements, no additional ethical approval is required for this research. Furthermore, consent was obtained from participants during the initial data collection phase. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Yili Zhang, Xiangyun Guo, Kai Sun and Liang Wang contributed equally to this work. References