Analysis and Validation of Biomarkers of immune cell related genes in postmenopausal osteoporosis

doi:10.21203/rs.3.rs-3231585/v1

Download PDF

Research Article

Analysis and Validation of Biomarkers of immune cell related genes in postmenopausal osteoporosis

https://doi.org/10.21203/rs.3.rs-3231585/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Objective

Postmenopausal osteoporosis (PMOP) is a common metabolic inflammatory disease. In conditions of estrogen deficiency, chronic activation of the immune system leads to a hypo-inflammatory phenotype and alterations in its cytokine and immune cell profile, although immune cells play an important role in the pathology of osteoporosis, studies on this have been rare. Therefore, it is important to investigate the role of immune cell-related genes in PMOP.

Methods

PMOP-related datasets were downloaded from the Gene Expression Omnibus database. Immune cells scores between high bone mineral density (BMD) and low BMD samples were assessed based on the single sample gene set enrichment analysis (ssGSEA) method. Subsequently, weighted gene co-expression network analysis (WGCNA) was performed to identify modules highly associated with immune cells and obtain module genes. Differential analysis between high BMD and low BMD was also performed to obtain differentially expressed genes (DEGs). Module genes are intersected with DEGs to obtain candidate genes, and functional enrichment analysis was performed. Machine learning methods were used to filter out the signature genes. The receiver operating characteristic (ROC) curves of the signature genes and the nomogram were plotted to determine whether the signature genes can be used as a molecular marker. Gene set enrichment analysis (GSEA) was also performed to explore the potential mechanism of the signature genes. Finally, RNA expression of signature genes was validated in blood samples from PMOP patients and normal control by real-time quantitative polymerase chain reaction (qRT-PCR).

Results

Our study of PMOP patients identified differences in immune cells (Activated dendritic cell, CD56bright natural killer cell, Central memory CD4 T cell, Effector memory CD4 T cell, Mast cell, Natural killer T cell, T follicular helper cell, Type 1 T helper cell and Type 17 T helper cell) between high and low BMD patients. We obtained a total of 73 candidate genes based on modular genes and differential genes, and obtained five signature genes by least absolute shrinkage and selection operator and random forest model screening. ROC, principal component analysis and t-distributed stochastic neighbor embedding down scaling analysis revealed that the five signature genes had good discriminatory ability between high and low BMD samples. A logistic regression model was constructed based on five signature genes, and both ROC and column line plots indicated that the model accuracy and applicability were good. Five signature genes were found to be associated with proteasome, mitochondria, and lysosome by GSEA. The qRT-PCR results showed that the expression of the signature genes was significantly different between the two groups.

Conclusion

HIST1H2AG, PYGM, NCKAP1, POMP and LYPLA1 might play key roles in PMOP and be served as the biomarkers of PMOP.

Postmenopausal osteoporosis

immune cells

signature genes

diagnosis

Postmenopausal osteoporosis (PMOP) is a systemic disease characterized by low bone mass and destruction of bone microstructure, leading to an increased risk of fractures ¹. According to the World Health Organization (WHO) definition of osteoporosis based on bone mineral density (BMD): osteoporosis is defined as BMD 2.5 SD or more below the mean for premenopausal women ². With the aging global population, the incidence of PMOP has increased dramatically. PMOP is a global public health problem that not only seriously affects the health and quality of life of patients, but also places a huge economic burden on healthcare systems ³. It has been reported that approximately 50% of postmenopausal women worldwide are affected by osteoporosis, and the prevalence of fractures in patients with osteoporosis is as high as 40% ⁴. In addition to causing great suffering to patients, osteoporotic fractures can lead to disability, causing greater pressure on a country's public health system. In addition, associated complications such as pressure ulcers, fall pneumonia, and pulmonary embolism have been noted, resulting in nearly 30% of patients dying within one year of being diagnosed with PMOP. In the United States, it is estimated that by 2020, the population over 50 years of age will increase to 1213 million, of which more than 3 million suffer from some kind of osteoporotic fracture. This means that the estimated cost of treating the disease in the United States alone will approach $25 billion annually ^{2, 4–6}. The incidence of fractures is higher in women than in men, and from the age of 50 onwards, the incidence of fractures is twice as high in women as in men, and low bone mass is an important risk factor for fractures in postmenopausal women ⁷. As life expectancy increases around the world, it is noted that the impact of this age-related disease will continue to expand. Because of its high incidence and severe complications, PMOP has become a worldwide public health problem ⁸.

PMOP is a common skeletal disorder associated with reduced estrogen in postmenopausal women, and estrogen deficiency is a key pathogenesis. The pathogenesis of PMOP may be closely related to endocrine disorders, genetic factors, oxidative stress, and nutritional imbalance, with decreased estrogen being the key factor ⁸. After menopause, estrogen levels decrease in women, and the inhibitory effect of estrogen on osteoclasts (OCs) is weakened. As a result, the number of OCs increases, apoptosis decreases, and their lifespan is extended, thereby disrupting the balance between OCs and osteoblasts (OBs) and enhancing bone resorption ⁹. PMOP is characterized by a high rate of bone turnover. Bone resorption is increased, as is bone formation, but at a higher rate in the former than in the latter. In healthy bone, bone homeostasis is regulated by a complex coordination between the interacting components of BRC, receptor activator of NF-κB ligand (RANKL) and osteoprotegerin (OPG). OPG is a decoy receptor for RANKL. RANKL secreted by OBs interacts with rank, triggering the differentiation of precursor osteoclasts into larger multinucleated active OCs. However, to control bone resorption, OBs also secrete OPG, which competitively inhibits the RANKL-receptor activator of NF-κB (RANK) interaction ^{10, 11}. In addition to regulating endocrine, paracrine and mechanical factors and modulating bone metabolism, the immune system plays a regulatory role in bone remodeling and maintenance of bone mass. Innate immune cells regulate osteoporosis in large part by influencing the RANK/RANKL/OPG axis by producing a variety of pro-inflammatory mediators and regulating cells that are critical in causing osteoporosis ¹². There is increasing evidence that both innate and adaptive immune cells promote the pathogenesis of osteoporosis by producing proinflammatory mediators ^{13, 14}. Estrogen affects the proliferation, differentiation, activation, and homing of many immune cells. Under conditions of estrogen deficiency, chronic activation of the immune system can lead to a low inflammatory phenotype characterized by changes in cytokines and immune cell profiles, both in rodents and humans. It has been reported that the RANK/RANKL/OPG pathway is the most important signal mediating the effects of estrogen in bone. Many immune cells pass through OPG/RANKL, inflammatory cytokines such as interleukin-6 (IL-6) and tissue necrosis factors α (TNF-α). And other mediators secreted by immune cells directly or indirectly affect bone cells. Some experimental studies have shown that immune cell-mediated factors have a significant impact on bone cells, and blocking the interaction between immune cells and bone cells can eliminate or alleviate ovariectomy (OVX)-induced bone loss and OVX-induced delayed bone repair. In addition to the direct negative effects of estrogen deficiency on bones, changes in immune cells are also involved in the pathogenesis of postmenopausal osteoporosis ¹⁵.

Despite the important role of immune cells in the pathology of osteoporosis, relatively few studies are available. Therefore, this study aimed to screen immune cell related genes in postmenopausal osteoporosis from the perspective of immune cell infiltration by bioinformatics analysis and successfully constructed a diagnostic model to provide new biomarkers for the diagnosis and treatment of osteoporosis.

2.1. Data source

GSE56814, GSE56815 and GSE13850 datasets were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE56814, based on the GPL5175 platform, contains gene expression profiles of blood monocytes from 26 subjects with high bone mineral density (BMD) and 16 subjects with low BMD, all of postmenopausal females. GSE56815, based on the GPL96 platform, contains gene expression profiles of blood monocytes from 40 postmenopausal females, 20 with high BMD and 20 with low BMD. GSE13850, based on the GPL96 platform, contains gene expression profiles of whole blood samples from 40 females with postmenopausal females, of whom 20 had high BMD and 20 low BMD. To reduce batch effects from non-biological technical biases involved in the expression profiles of GSE56814 and GSE56815 as a combined dataset, the “ComBat” algorithm25 and the SVA package were applied. Principal component analysis (PCA)on batch pre-calibration and post-calibration data.

2.2. Single sample gene set enrichment analysis (ssGSEA)

To investigate the variability of immune infiltration between high BMD and low BMD, this study was implemented to estimate the level of infiltration of immune cell types using ssGSEA in the R package gsva (version 1.42.0). ssGSEA allows the definition of an enrichment score that represents the absolute enrichment of gene sets in each sample within a given dataset. The level of infiltration of 28 immune cells is estimated using the combined dataset as input. The immune cell ssGSEA score heat map was plotted for each sample using the pheatmap R package (version 1.0.12). A wilcoxon rank sum test was used to compare the differences between low BMD and high BMD for different immune cells in the combined dataset, and box plots were drawn using the ggpubr package (version 0.4.0) to visualize the comparative results ¹⁶.

2.3. Weighted gene co-expression network analysis (WGCNA)

WGCNA is a bioinformatics tool for building gene co-expression networks and mining network modules closely associated with disease ¹⁷. In this study, the R package WGCNA (version 1.7-3) was used to filter the combined data for outliers and then construct a co-expression network. The differential immune cell infiltration abundance of high BMD and low BMD were used as trait data for WGCNA to search for genes associated with BMD ¹⁷.

First, the optimal soft-thresholding power (β) value was calculated based on the expression data of the mRNAs in all samples. Hierarchical clustering and dynamic tree cutting functions were used to detect modules (minimum size = 100). The name of each co-expression module was named as the color. Next, each first principal component of each gene module was identified as module migengenes (MEs). We then looked for correlations between traits and MEs and identified the modules associated with BMD as key modules (r > 0.6 and p < 0.05) ¹⁸.

Gene significance (GS) was used to identify relationships between genes and traits, and module membership (MM) represented the correlation between MEs and gene expression profiles. Genes with high GS and MM in key modules are highly correlated with traits. In this study we selected genes with |GS|>0.2 and |MM|>0.8 as pivotal genes in the key module ¹⁸.

2.4. Identification of differentially expressed genes

The R language limma package (version 3.52.2) was used to compare the differentially expressed genes (DEGs) between the combined data low BMD and high BMD groups. The screening condition was p.adjust < 0.05.

2.5. Identification of candidate genes and functional enrichment analysis

In this study, the R package (VennDiagram, version 1.7.1) was used to obtain candidate genes by taking intersections of the DEGs and key module genes obtained above. To understand the function of the candidate genes, Gene Ontology (GO) process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using the clusterProfiler R package (version 4.4.4) and the enrichment results were visualized using the enrichplot package (version 1.16.1) to create histograms and bubble plots ¹⁹.

2.6. Selection of signature genes

We used gene expression matrices from high BMD and low BMD samples for candidate genes and constructed least absolute shrinkage and selection operator (Lasso) models and random forest (RF) models to screen for disease-critical signature genes. We use the glmnet R package (version 4.1-4) with the parameters family for binomial and type.measure for class to implement LASSO logistic regression ¹⁹. The RF method in caret R package (version 6.0–92) was also used to screen candidate genes for signature genes ²⁰. Using the pROC (version 1.17.0.1) package, the ability of the Lasso and RF models to diagnose patients in the training and validation (GSE13850) sets was assessed by calculating the area under the receiver operating characteristic curve (AUC) ²⁰. The final signature genes were obtained by intersecting the signature genes of lasso and RF for subsequent analysis.

2.7. Construction and evaluation of diagnostic models

The distribution of signature genes between high and low BMD samples was analyzed using PCA and t-distributed stochastic neighbor embedding (t-SNE) dimensionality reduction algorithms on the merged dataset ²¹. Receiver Operating Characteristic Curves (ROC) analysis was performed using the pROC (version 1.18.0) R package, and ROC validation was done on the signature genes of the combined dataset to determine the accuracy of the genes in distinguishing between high and low BMD.

After combining the clinical data with the risk values obtained in the previous step, we A multi-factor logistic regression model (LRM) of risk values and clinical factors was constructed using the R package rms (version 6.3-0) ²¹. Nomogram was drawn from the LBM, and calibration curves were plotted to assess the predictive power of the nomogram model ²². In addition, the accuracy of the logistic regression model was verified by plotting the ROC using the GSE13850 dataset.

2.8. Gene set enrichment analysis (GSEA)

To understand the function of the signature genes, GSEA was performed using the clusterProfiler R package (version 4.4.4) to look for functional and related pathways between high and low risk groups. |NES|>1, p.adjust < 0.05 and q.value < 0.2 were considered significant enrichment ²³.

2.9. RNA extraction and real-time quantitative polymerase chain reaction (qRT-PCR)

This study was in accordance with the Declaration of Helsinki (as revised in 2013), and it was approved by the Medical Ethics Committee of Shenzhen Bao'an Hospital of Traditional Chinese Medicine (KY-2023-058-01). All participants signed the relevant consent forms. Postmenopausal women older than 60 years admitted to the hospital's Department of Rehabilitation between December 2022 and January 2023 were included in this study. Patients were excluded from the study if they were found to have any of the following: severe cardiopulmonary disease; untreated coagulopathy; systemic or spinal infection, history of malignancy, untreated bone metabolic disease other than osteoporosis, digestive dysfunction, severe liver and / or kidney dysfunction, history of spine or hip surgery, 25 hydroxyvitamin d < 23 nmol/L, small amounts of serum calcium, and abnormal parathyroid levels. A total of 10 PMOP patients and 10 healthy subjects were recruited. Total RNA was extracted from the participants' peripheral blood samples (10 PMOP vs 10 normal control) using the TRIzol kit. Reverse transcription was conducted using SureScript First strand cDNA synthesis kit (Genecopoeia, USA). The signature genes level was assessed using Universal Blue SYBR Green qPCR Master Mix (Yeasen, Shanghai, China) following the instructions. Relative expression levels were normalized to GAPDH. Primers are available in TABLE S1. The relative expression of signature genes was calculated by 2 − ΔΔCt method with the normalization to GAPDH.

2.10. Statistical analysis

All analyses were conducted using the R programming language, and the data from different groups were compared by the Wilcoxon test. If not specified above, a p value less than 0.05 was considered statistically significant.

3.1. Data Preprocessing

We performed background correction and quantile normalization on dataset GSE56814 and dataset GSE56815, and then combined the 2 datasets. Gene expression matrices for 82 samples (high BMD = 46, low BMD = 36) were obtained after removing the batch effect, and box plots and PCA results showed good results in removing the batch effect (Fig. 1A-B).

3.2. Immune cell infiltration landscape analysis

By examining the differences in immune infiltration between high BMD and low BMD, a total of nine (Activated dendritic cell, CD56bright natural killer cell, Central memory CD4 T cell, Effector memory CD4 T cell, Mast cell, Natural killer T cell, T follicular helper cell, Type 1 T helper cell and Type 17 T helper cell) out of 28 immune cells were found to be significantly different between high BMD and low BMD (p < 0.05) (Fig. 2A-B).

3.3. WGCNA

Firstly, the cluster analysis showed that the samples clustered well overall and No need to eliminate samples (Fig. 3A). The soft thresholding power β was set at 11 in the subsequent analysis, because the scale independence reached 0.85 (Fig. 3B). A total of 7 modules were obtained (Fig. 3C), of which the red module (404 genes) and turquoise module (3158 genes) were the most relevant modules for immune cells (Fig. 3D). The criteria of |MM|>0.8 and |GS|>0.2 were used to obtain 68 key module genes in the red module and 244 key module genes in the turquoise module (Fig. 3E-F).

3.4. Identification of DEGs

A total of 1394 differentially expressed genes were obtained from the combined data between the Low BMD and High BMD groups, of which 873 genes were up-regulated in the Low BMD group and 521 genes were down-regulated in the Low BMD. Volcano and heat maps showed the expression of differential genes (Fig. 4A-B).

3.5. Identification of candidate genes and functional enrichment analysis

A total of 73 candidate genes were obtained by taking the intersection of the differential genes and key module genes obtained above (Fig. 5A). A total of 85 (p.adjust < 0.05, count > 2) GO processes were enriched for 73 candidate genes (p.adjust < 0.05, count > 2), of which 26 were biological process (BP) terms, 51 were cellular component (CC) terms and 7 were molecular function (MF) terms (Fig. 5B). Candidate genes were enriched with GO terms such as establishment of protein localization to membrane, Arp2/3complex-mediated actin nucleation and proteasome assembly. Candidate genes in KEGG were mainly enriched to the proteasome pathway (p.adjust < 0.05 & count > 2) (Fig. 5C).

3.6. Identification of signature genes

We used gene expression matrix data from 82 samples (high BMD = 46, low BMD = 36), for 73 candidate genes, to construct the Lasso and RF models. Lasso model was filtered to obtain 29 signature genes (Fig. 6A), and the ROC of the training and validation sets were plotted based on the 29 signature genes, and it was found that the AUC values of both the combined data and the validation set were equal to 1, indicating the accuracy and applicability of the model (Fig. 6B-C). The RF model obtained six signature genes and plotted the ROC of the training and validation sets based on the six signature genes, and found that the AUC values of both the combined and validation sets were equal to 1, indicating that the accuracy and applicability of the model were also good (Fig. 6D-G). We selected the intersection genes of the Lasso model and the RF model as the identified signature genes. A total of five signature genes (HIST1H2AG, PYGM, NCKAP1, POMP and LYPLA1) were obtained in this study (Fig. 6H).

3.7. Construction and evaluation of diagnostic models

PCA and t-SNE analyses revealed that five signature genes could significantly distinguish patients in the high BMD and low BMD (Fig. 7A-B). The AUC values for all five signature genes in the combined data were greater than 0.7, suggesting that these signature genes have some accuracy in distinguishing high BMD from low BMD samples and that they can be used as diagnostic genes for our disease (Fig. 7C).

Based on the risk values obtained in the previous step and combined with the clinical data, we constructed a multi-factor logistic regression model with risk values and clinical factors (Fig. 7D). The c-index of the nomogram is 0.847, which indicates that the prediction effect of the nomogram is good, and the AUC value of the validation set is greater than 0.8, which also indicates that the model has high accuracy (Fig. 7E-G).

3.8. GSEA

We performed GO and KEGG single gene GSEA enrichment analysis on five signature genes and found that five signature genes were associated with immune macromolecules such as proteasome, mitochondria and lysosome. POMP and PYGM genes were enriched in the proteasome complex pathway in both the high and low BMD groups. LYLPLA1 gene and POMP gene were enriched in the mitochondrial protein-containing complex in both the high and low BMD groups. PYGM gene was enriched in the mitochondrial protein-containing complex in the high BMD group. The figures show the score curves for the top 10 GO enrichment results and the top 10 KEGG enrichment results for the five signature genes (Fig. 8A-D, FIGURE S1-2).

3.9. Expression levels of signature genes

The qRT-PCR results showed that the expression of the signature genes was significantly different between the two groups. POMP was expressed up-regulated in the PMOP group and other genes (HIST1H2AG, PYGM, NCKAP1 and LYPLA1) were expressed down-regulated in the PMOP group compared to the normal control (Fig. 9A-E).

PMOP is a chronic systemic metabolic disease associated with decreased estrogen levels. The symptoms of PMOP are characterized by imbalances in bone metabolism caused by various factors such as genetics, epigenetics, and the environment. The decline in estrogen levels after menopause exacerbates bone resorption ¹. Estrogen deficiency is a multifaceted disease affecting the whole body that not only affects postmenopausal osteoporosis but also interacts with various immune cells, leading to a chronic low-grade proinflammatory phenotype under estrogen deficiency ^{24, 25}. Clinical data on the immunophenotype of postmenopausal patients suggested that postmenopausal women had elevated levels of inflammatory cytokines, such as IL-1 β, IL-6 and TNF-α ^26–29. Estrogen deficiency is a driver of increased cytokine expression. In addition to higher levels of inflammatory mediators, changes in immune cell numbers occurred in postmenopausal women. One study found that in PMOP patients with bone fractures, T cells were more likely to express TNF-α ³⁰. Another study showed that circulating T cells and monocytes were increased in postmenopausal women ²⁸. During PMOP, B cells activated by estrogen deficiency and proinflammatory conditions contribute to increased bone resorption by secreting enhanced levels of granulocyte colony stimulating factor (G-CSF) and RANKL^31–34. G-CSF can promote neutrophil differentiation, potentially contributing to the increased number of neutrophils in estrogen deficient subjects ³². In addition, patients with PMOP had a higher neutrophil / lymphocyte ratio in peripheral blood ³⁵ and an increased number of mast cells in the bone marrow ³⁶. In summary, PMOP patients presented with a chronic low-grade inflammatory phenotype, altered cytokine expression and immune cell profiles. In recent years, research on the cellular and molecular mechanisms of postmenopausal osteoporosis has found that a large number of genes and proteins are related to the occurrence and development of the disease. However, the impact of these molecular mechanisms on the occurrence and development of osteoporosis is unclear ³⁷. Therefore, the research field needs to further explore the genes and/or proteins associated with postmenopausal osteoporosis and examine the molecular mechanisms of their interaction.

First of all, this study of PMOP patients identified differences in immune cells (Activated dendritic cell, CD56 bright natural killer cell, Central memory CD4 T cell, Effector memory CD4 T cell, Mast cell, Natural killer T cell, T follicular helper cell, Type 1 T helper cell and Type 17 T helper (Th17) cell) between high and low BMD patients. One study showed that activated dendritic cells can develop functional OC in a RANKL/RANK dependent manner after interaction with T cells (CD4 T cell cells or T helper cell) ³⁸. Several clinical studies have shown that levels of circulating T cells are normal or increased in postmenopausal osteoporosis patients with altered cytokine expression profiles ¹⁵. In addition, it was found that Th17 cells were increasingly present in the bone marrow of OVX mice ³⁹. Th17 cells in bone marrow were proved to increase the recruitment of inflammatory monocytes as osteoclast precursor somatic cell to bone marrow ⁴⁰. There was indirect osteoclastic effect of IL-17 secreted by Th17 cells, and IL-17 can upregulate osteoclast precursor cells, thereby increasing their sensitivity to RANKL stimulation ^{41, 42}. In the inflammatory state, osteoclasts can affect CD4 T cells, thereby affecting the production of TNFα ^{43, 44}.

Based on modular genes and differential genes, we determined that the top five signature genes were HIST1H2AG, PYGM, NCKAP1, POMP and LYPLA1. We collected serum from PMOP patients and postmenopausal subjects with normal BMD to detect differential genes based on bioinformatics analysis using qRT-PCR. The experimental results showed that the genes were statistically significantly different from postmenopausal healthy women, consistent with our bioinformatics findings. In order to further elucidate the significance of this study, it is important to elucidate the role of each of these five genes. The main role of PYGM is to provide sufficient energy for muscle contraction. It is expressed in tissues other than muscle, such as brain, lymphoid tissues, and blood. PYGM is important not only in glycogen metabolism but also in various processes such as insulin and glucagon signaling pathways, insulin resistance, necroptosis, immune response, and phototransduction. PYGM has been associated with several pathological conditions such as McArdle disease, schizophrenia and cancer. Among them, retinopathy, one of the symptoms of McArdle disease, is associated with Ca2 + levels in the retinal pigment epithelium and cone photoreceptors, but the exact mechanism remains unclear ⁴⁵. Additional experimental analysis is needed to understand the potential role of pygm in regulating aspects of cellular Ca2 + levels ⁴⁶. In addition, the gene PYGM was found to play a central role in metabolic changes in response to heat stress ⁴⁷. PYGM is also expressed in T lymphocytes ⁴⁸. PYGM pathway plays an important role in regulating the immune function of T cells ⁴⁹. Activation of PYGM through the pygmil-2 pathway is key to IL-2-stimulated T cell migratory and proliferative responses ⁵⁰. A web-based meta-analysis has identified important differentially expressed genes in women with osteoporosis, including PYGM ⁵¹. Lycium barbarum (LC) has a powerful role in inhibiting bone loss. Kukoamine A (KuA) is a biologically active compound extract of LC that has an anti-osteoporosis effect. PYGM is the target gene with the strongest correlation with the active components of LC and is regulated by KuA ⁵². Protein corresponding to POMP gene directly related to the proteasome precursor complex, which acts on the assembly and final maturation steps of eukaryotic proteasomes ⁵³. A study performed in vitro and in vivo experiments showing that pomp knockdown has a robust antiproliferative effect in cancer cells ⁵⁴. In terms of immunity and inflammation, POMP can be upregulated by interferon to enhance major histocompatibility complex (MHC) - I dependent antigen presentation, thereby enhancing specific T cell responses ⁵⁵. POMP is a chaperone for proteasome assembly, and AD mutations in POMP cause a form of PRAAS with prominent immunodeficiency referred to as POMP-related autoinflammation and immune dysregulation (PRAID) ⁵⁶. The protein expressed by NCKAP1 is found in sporadic Alzheimer’s disease (AD) as part of the WAVE complex along with ABI1-2, BRK1, CYFIP1-2, and WASF1-2 proteins ⁵⁷. NCKAP1 regulates various intracellular processes such as apoptosis, migration, and invasion, and plays an essential role in disease pathogenesis ⁵⁷. Both PYGM and NCKAP1 genes are strongly implicated in carcinogenesis. NCKAP1 may function as an oncogene in a variety of cancer types. NCKAP1 expression is highly tissue-specific, and its expression has been found in multiple cancer types ^{22, 58, 59}. Studies in HIST1H2AG and LYPLA1 were less frequently reported. The expression level of lnc-HIST1H2AG-6 gene were significantly correlated with some features of T2DM ¹⁶. High expression of LYPLA1 was associated with a poor prognosis in lung adenocarcinoma (LUAD) ⁵⁹. In addition, we did GSEA analysis and found that these five signature genes were associated with proteasome, mitochondria and lysosome. Although some of the above genes are related to immunity, there have been no reports of direct correlation between these 5 genes in PMOP. Therefore, we found for the first time that these 5 immune cell related genes are differentially expressed in PMOP. In order to identify the discriminative ability, accuracy, and applicability of these five signature genes for PMOP, we conducted many additional analyses, such as ROC, principal component analysis, t-distributed stochastic neighbor embedding down scaling analysis and logistic regression analysis. These analyses indicated that the five model genes have good disease diagnosis capabilities.

In summary, by combining various biological information mining analyses, our study found abnormal expression of genes in postmenopausal osteoporosis patients. This study identified immune cell related signature genes for the first time in PMOP, and validated these genes. Additionally, there is a drawback to this study, as the clinical application of characteristic genes requires more sample data support. In future research, our researchers will continue to explore and pay attention to the roles of these genes. This discovery brings researchers closer to determining the molecular mechanism of osteoporosis in postmenopausal women.

5.1 Acknowledgments

Our sincere gratitude must be expressed to Shenzhen Science and Technology Innovation Committee and all contributors who offered any assistance in the article.

5.2 Authors' contributions

Conceptualization, visualization and material preparation were performed by LHC, JJQ, and YZ. Data collection and formal analysis were performed by LHC, JJQ, and YZ. The first draft of the manuscript was written by LHC, JJQ, and YZ. XSL and YZ reviewed andedited the original draft and supervised over the study. All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

5.3 Funding

This work was funded and supported by Shenzhen Science and Technology Innovation Committee (grant No. JCYJ20190808100818959).

5.4 Ethics Approval and Consent to Participate

The study involving human participants was approved by the Medical Ethics Committee of Shenzhen Bao'an Hospital of Traditional Chinese Medicine (KY-2023-058-01). The patients provided their written informed consent to participate in this study.

5.5 Availability of data and materials

The datasets analyzed during the current study are available in the GEO (https://www.ncbi.nlm.nih.gov/geo/) [GSE56814, GSE56815 and GSE13850] repository.

5.6 Consent for publication

Not applicable.

5.7 Competing interests

No competing interests exist.

Avioli LV. Senile and postmenopausal osteoporosis. Adv Intern Med 1976; 21:391-415.
Lewiecki EM. Romosozumab, clinical trials, and real-world care of patients with osteoporosis. Ann Transl Med 2020; 8 (15):974.
Häussler B, Gothe H, Göl D, Glaeske G, Pientka L, Felsenberg D. Epidemiology, treatment and costs of osteoporosis in Germany--the BoneEVA Study. Osteoporos Int 2007; 18 (1):77-84.
Rachner TD, Khosla S, Hofbauer LC. Osteoporosis: now and the future. Lancet 2011; 377 (9773):1276-87.
Black DM, Rosen CJ. Clinical Practice. Postmenopausal Osteoporosis. N Engl J Med 2016; 374 (3):254-62.
Burge R, Dawson-Hughes B, Solomon DH, Wong JB, King A, Tosteson A. Incidence and economic burden of osteoporosis-related fractures in the United States, 2005-2025. J Bone Miner Res 2007; 22 (3):465-75.
Goulding A, Jones IE, Taylor RW, Manning PJ, Williams SM. More broken bones: a 4-year double cohort study of young girls with and without distal forearm fractures. J Bone Miner Res 2000; 15 (10):2011-8.
Li J, Chen X, Lu L, Yu X. The relationship between bone marrow adipose tissue and bone metabolism in postmenopausal osteoporosis. Cytokine Growth Factor Rev 2020; 52:88-98.
Liu S, Chen T, Wang R, et al. Exploring the effect of the "quaternary regulation" theory of "peripheral nerve-angiogenesis-osteoclast-osteogenesis" on osteoporosis based on neuropeptides. Front Endocrinol (Lausanne) 2022; 13:908043.
McNamara LM. Osteocytes and Estrogen Deficiency. Curr Osteoporos Rep 2021; 19 (6):592-603.
Walsh MC, Takegahara N, Kim H, Choi Y. Updating osteoimmunology: regulation of bone cells by innate and adaptive immunity. Nat Rev Rheumatol 2018; 14 (3):146-56.
Saxena Y, Routh S, Mukhopadhaya A. Immunoporosis: Role of Innate Immune Cells in Osteoporosis. Front Immunol 2021; 12:687037.
Mundy GR. Osteoporosis and inflammation. Nutr Rev 2007; 65 (12 Pt 2):S147-51.
Zupan J, Jeras M, Marc J. Osteoimmunology and the influence of pro-inflammatory cytokines on osteoclasts. Biochem Med (Zagreb) 2013; 23 (1):43-63.
Fischer V, Haffner-Luntzer M. Interaction between bone and immune cells: Implications for postmenopausal osteoporosis. Semin Cell Dev Biol 2022; 123:14-21.
Jiang H, Lou P, Chen X, Wu C, Shao S. Deregulation of lncRNA HIST1H2AG-6 and AIM1-3 in peripheral blood mononuclear cells is associated with newly diagnosed type 2 diabetes. BMC Med Genomics 2021; 14 (1):149.
Zhai X, Xue Q, Liu Q, Guo Y, Chen Z. Colon cancer recurrence‑associated genes revealed by WGCNA co‑expression network analysis. Mol Med Rep 2017; 16 (5):6499-505.
Li A, He J, Zhang Z, et al. Integrated Bioinformatics Analysis Reveals Marker Genes and Potential Therapeutic Targets for Pulmonary Arterial Hypertension. Genes (Basel) 2021; 12 (9):1339.
Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021; 2 (3):100141.
Zhao S, Bao Z, Zhao X, Xu M, Li MD, Yang Z. Identification of Diagnostic Markers for Major Depressive Disorder Using Machine Learning Methods. Front Neurosci 2021; 15:645998.
Yang M, Zheng H, Xu K, et al. A novel signature to guide osteosarcoma prognosis and immune microenvironment: Cuproptosis-related lncRNA. Front Immunol 2022; 13:919231.
Fang Q, Chen H. Development of a Novel Autophagy-Related Prognostic Signature and Nomogram for Hepatocellular Carcinoma. Front Oncol 2020; 10:591356.
Zhao X, Zhang L, Wang J, et al. Identification of key biomarkers and immune infiltration in systemic lupus erythematosus by integrated bioinformatics analysis. J Transl Med 2021; 19 (1):35.
Pacifici R. Estrogen, cytokines, and pathogenesis of postmenopausal osteoporosis. J Bone Miner Res 1996; 11 (8):1043-51.
Ralston SH. Analysis of gene expression in human bone biopsies by polymerase chain reaction: evidence for enhanced cytokine expression in postmenopausal osteoporosis. J Bone Miner Res 1994; 9 (6):883-90.
Cohen-Solal ME, Graulet AM, Denne MA, Gueris J, Baylink D, de Vernejoul MC. Peripheral monocyte culture supernatants of menopausal women can induce bone resorption: involvement of cytokines. J Clin Endocrinol Metab 1993; 77 (6):1648-53.
Zheng SX, Vrindts Y, Lopez M, et al. Increase in cytokine production (IL-1 beta, IL-6, TNF-alpha but not IFN-gamma, GM-CSF or LIF) by stimulated whole blood cells in postmenopausal osteoporosis. Maturitas 1997; 26 (1):63-71.
Abildgaard J, Tingstedt J, Zhao Y, et al. Increased systemic inflammation and altered distribution of T-cell subsets in postmenopausal women. PLoS One 2020; 15 (6):e0235174.
Zupan J, Komadina R, Marc J. The relationship between osteoclastogenic and anti-osteoclastogenic pro-inflammatory cytokines differs in human osteoporotic and osteoarthritic bone tissues. J Biomed Sci 2012; 19 (1):28.
Pietschmann P, Grisar J, Thien R, et al. Immune phenotype and intracellular cytokine production of peripheral blood mononuclear cells from postmenopausal patients with osteoporotic fractures. Exp Gerontol 2001; 36 (10):1749-59.
Toni R, Di Conza G, Barbaro F, et al. Microtopography of Immune Cells in Osteoporosis and Bone Lesions by Endocrine Disruptors. Front Immunol 2020; 11:1737.
Zhang Z, Yuan W, Deng J, et al. Granulocyte colony stimulating factor (G-CSF) regulates neutrophils infiltration and periodontal tissue destruction in an experimental periodontitis. Mol Immunol 2020; 117:110-21.
Weitzmann MN. The Role of Inflammatory Cytokines, the RANKL/OPG Axis, and the Immunoskeletal Interface in Physiological Bone Turnover and Osteoporosis. Scientifica (Cairo) 2013; 2013:125705.
Li Y, Terauchi M, Vikulina T, Roser-Page S, Weitzmann MN. B Cell Production of Both OPG and RANKL is Significantly Increased in Aged Mice. Open Bone J 2014; 6:8-17.
Fang H, Zhang H, Wang Z, Zhou Z, Li Y, Lu L. Systemic immune-inflammation index acts as a novel diagnostic biomarker for postmenopausal osteoporosis and could predict the risk of osteoporotic fracture. J Clin Lab Anal 2020; 34 (1):e23016.
Fallon MD, Whyte MP, Craig RB, Jr., Teitelbaum SL. Mast-cell proliferation in postmenopausal osteoporosis. Calcif Tissue Int 1983; 35 (1):29-31.
Compston JE, McClung MR, Leslie WD. Osteoporosis. Lancet 2019; 393 (10169):364-76.
Alnaeeli M, Penninger JM, Teng YT. Immune interactions with CD4+ T cells promote the development of functional osteoclasts from murine CD11c+ dendritic cells. J Immunol 2006; 177 (5):3314-26.
Yu M, Pal S, Paterson CW, et al. Ovariectomy induces bone loss via microbial-dependent trafficking of intestinal TNF+ T cells and Th17 cells. J Clin Invest 2021; 131 (4):e143137.
Ciucci T, Ibáñez L, Boucoiran A, et al. Bone marrow Th17 TNFα cells induce osteoclast differentiation, and link bone destruction to IBD. Gut 2015; 64 (7):1072-81.
Kotake S, Udagawa N, Takahashi N, et al. IL-17 in synovial fluids from patients with rheumatoid arthritis is a potent stimulator of osteoclastogenesis. J Clin Invest 1999; 103 (9):1345-52.
Sato K, Suematsu A, Okamoto K, et al. Th17 functions as an osteoclastogenic helper T cell subset that links T cell activation and bone destruction. J Exp Med 2006; 203 (12):2673-82.
Ibáñez L, Abou-Ezzi G, Ciucci T, et al. Inflammatory Osteoclasts Prime TNFα-Producing CD4(+) T Cells and Express CX(3) CR1. J Bone Miner Res 2016; 31 (10):1899-908.
Madel MB, Ibáñez L, Ciucci T, et al. Dissecting the phenotypic and functional heterogeneity of mouse inflammatory osteoclasts by the expression of Cx3cr1. Elife 2020; 9:e54493.
Migocka-Patrzałek M, Elias M. Muscle Glycogen Phosphorylase and Its Functional Partners in Health and Disease. Cells 2021; 10 (4):883.
Vaclavik V, Naderi F, Schaller A, Escher P. Longitudinal case study and phenotypic multimodal characterization of McArdle disease-linked retinopathy: insight into pathomechanisms. Ophthalmic Genet 2020; 41 (1):73-8.
Sammad A, Luo H, Hu L, et al. Joint Transcriptome and Metabolome Analysis Prevails the Biological Mechanisms Underlying the Pro-Survival Fight in In Vitro Heat-Stressed Granulosa Cells. Biology (Basel) 2022; 11 (6):839.
de Luna N, Brull A, Lucia A, et al. PYGM expression analysis in white blood cells: a complementary tool for diagnosing McArdle disease? Neuromuscul Disord 2014; 24 (12):1079-86.
Llavero F, Luque Montoro M, Arrazola Sastre A, et al. Epidermal growth factor receptor controls glycogen phosphorylase in T cells through small GTPases of the RAS family. J Biol Chem 2019; 294 (12):4345-58.
Llavero F, Urzelai B, Osinalde N, et al. Guanine nucleotide exchange factor αPIX leads to activation of the Rac 1 GTPase/glycogen phosphorylase pathway in interleukin (IL)-2-stimulated T cells. J Biol Chem 2015; 290 (14):9171-82.
He H, Cao S, Niu T, et al. Network-Based Meta-Analyses of Associations of Multiple Gene Expression Profiles with Bone Mineral Density Variations in Women. PLoS One 2016; 11 (1):e0147475.
Luo L, Guan Z, Jin X, Guan Z, Jiang Y. Identification of kukoamine a as an anti-osteoporosis drug target using network pharmacology and experiment verification. Mol Med 2023; 29 (1):36.
Burri L, Höckendorff J, Boehm U, Klamp T, Dohmen RJ, Lévy F. Identification and characterization of a mammalian protein interacting with 20S proteasome precursors. Proc Natl Acad Sci U S A 2000; 97 (19):10348-53.
Zhang X, Schulz R, Edmunds S, et al. MicroRNA-101 Suppresses Tumor Cell Proliferation by Acting as an Endogenous Proteasome Inhibitor via Targeting the Proteasome Assembly Factor POMP. Mol Cell 2015; 59 (2):243-57.
Kloetzel PM. Antigen processing by the proteasome. Nat Rev Mol Cell Biol 2001; 2 (3):179-87.
Boyadzhiev M, Marinov L, Boyadzhiev V, Iotova V, Aksentijevich I, Hambleton S. Disease course and treatment effects of a JAK inhibitor in a patient with CANDLE syndrome. Pediatr Rheumatol Online J 2019; 17 (1):19.
Innocenti M, Zucconi A, Disanza A, et al. Abi1 is essential for the formation and activation of a WAVE2 signalling complex. Nat Cell Biol 2004; 6 (4):319-27.
Xiong Y, He L, Shay C, et al. Nck-associated protein 1 associates with HSP90 to drive metastasis in human non-small-cell lung cancer. J Exp Clin Cancer Res 2019; 38 (1):122.
Rai A, Greening DW, Xu R, Suwakulsiri W, Simpson RJ. Exosomes Derived from the Human Primary Colorectal Cancer Cell Line SW480 Orchestrate Fibroblast-Led Cancer Invasion. Proteomics 2020; 20 (14):e2000016.

No competing interests reported.

SupportingInformation.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Analysis and Validation of Biomarkers of immune cell related genes in postmenopausal osteoporosis

Status:

Version 1

Abstract

Objective

Methods

Results

Conclusion

Figures

1. Introduction

2. Materials and methods

2.1. Data source

2.2. Single sample gene set enrichment analysis (ssGSEA)

2.3. Weighted gene co-expression network analysis (WGCNA)

2.4. Identification of differentially expressed genes

2.5. Identification of candidate genes and functional enrichment analysis

2.6. Selection of signature genes

2.7. Construction and evaluation of diagnostic models

2.8. Gene set enrichment analysis (GSEA)

2.9. RNA extraction and real-time quantitative polymerase chain reaction (qRT-PCR)

2.10. Statistical analysis

3. Results

3.1. Data Preprocessing

3.2. Immune cell infiltration landscape analysis

3.3. WGCNA

3.4. Identification of DEGs

3.5. Identification of candidate genes and functional enrichment analysis

3.6. Identification of signature genes

3.7. Construction and evaluation of diagnostic models

3.8. GSEA

3.9. Expression levels of signature genes

4. Discussion

5. Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1