Bioinformatics Analysis of the Key Genes and Pathways in Multiple Myeloma

Introduction

Multiple myeloma is a B-cell malignancy characterized by the malignant proliferation of clonal plasma cells in bone marrow,¹ It is the second most common hematological malignancy after lymphoma.² The age-standardised rate (ASR) of multiple myeloma incidence was 1·78 (95% UI 1·69-1·87) per 100 000 people globally and mortality was 1·14 (95% UI 1·07-1·21) per 100 000 people globally in 2020.³ Although the treatment strategies for multiple myeloma,^4,5 such as proteasome inhibitors, immunomodulatory drugs, autologous stem cell transplantation and cellular immunotherapy, continue to improve, due to the characteristics of drug resistance and frequent recurrence of multiple myeloma, multiple myeloma is still incurable and the 5-year survival rate is approximately 50%^6,7. At present, there is no clinical treatment plan to cure multiple myeloma, mainly through early chemotherapy and other methods to control the disease, but the treatment effect is not ideal. Studies have shown that 80–90% of patients with multiple myeloma have been troubled by active bone diseases, such as bone pain, spontaneous fracture, hypercalcemia, etc.⁸ because the delicate balance between bone formation and bone destruction has been broken.⁹ This seriously affects the quality of life of myeloma patients and increases their mortality. Although some studies^10,11 have found that there is a bidirectional cross talk between bone cells and immune cells. Immune cells participate in the physiological and pathological activities of bone,¹² but the effects on bone of immune cells (such as macrophages, granulocytes and innate immune pathways) remain unclear.¹³ Therefore, it is urgent to implement effective prevention, treatment and prognosis monitoring for patients with multiple myeloma. With the development of the biological big data era, it is possible to analyze biological information containing a large amount of data through the combination of biology, computer science and information technology. High throughput microarray is widely used in molecular diagnosis, tumor classification, prognosis and predictive analysis of tumors.¹⁴ It is expected to provide a new way to study the pathogenesis and molecular diagnosis of multiple myeloma (MM). In this study, bioinformatics technology was used to analyze the gene expression microarray data of patients with multiple myeloma (MM) and healthy control group, mine the differentially expressed genes of patients with multiple myeloma (MM) and healthy control group, screen the key genes through differential expression gene analysis, functional enrichment and survival analysis, and explore the changes of gene pathways and functions related to multiple myeloma (MM), The aim is to provide new research ideas and theoretical basis for the early diagnosis and gene therapy of multiple myeloma (MM).

Methods and Materials

Microarray Data

In order to explore the important genes related to the occurrence and development of multiple myeloma, the original gene expression profile of GSE146649 data set was obtained from the comprehensive gene expression database (GEO) (https://www.ncbi.nlm.nih.gov/geo/). The mRNA expression profile of GSE146649 data set was detected using GPL570 platform (Affymetrix human genome U133 Plus 2.0 array [hg-u133_plus_2] Affymetrix human genome U133 Plus 2.0 array). In the study, 40 samples from the GSE146649 data set were selected for analysis, including samples initially diagnosed as multiple myeloma, samples of recurrent multiple myeloma, samples of fully remitted multiple myeloma and samples of normal people.

Identification of Differentially Expressed Genes

The original data set of GSE146649 is analyzed by using GEO2R software (https://www.ncbi.nlm.nih.gov/geo/geo2r/), a network tool using R’s “limma” software package to screen DEGs between multiple myeloma and normal tissues. Gene annotation and corresponding data files of DEGs were extracted by R software. The threshold value was |log2foldchange| >1.0 and adj.P.Val <0.05.

GO Analysis and KEGG Analysis

In order to further systematically analyze the functional annotation of differentially expressed genes (DEGs), go enrichment analysis and KEGG pathway analysis of differentially expressed genes (DEGs) were performed with David, an online biological information annotation database. The q value <0.05 was used as the cutoff value.

Construction of Protein-Protein Interaction Network

In order to explore the interaction between DEGs, the protein interaction (PPI) network of multiple myeloma differentially expressed genes (DEGs) was constructed by using string¹⁵ and Cytoscape online analysis tools, so as to find the hub gene.

Screening of Key Modules

The data set obtained from string analysis is imported into Cytoscape, and Cytohubba, a plug-in of Cytoscape,¹⁶ is used to select the central gene of biological network. In this study, the first 15 genes were selected as the central genes.

Coexpression and Module Function Analysis

Weighted gene co-expression network analysis (WGCNA) is a biological network that describes the correlation between differentially expressed genes (DEG).¹⁷ Co-expression analysis is a multidirectional network in which each node represents a gene. After evaluating the expression profile, the co-expression network was constructed by using the DEGs of GSE146649 dataset. The co-expression network of DEGs was constructed by R-package WGCNA. In addition, the Pearson correlation coefficient of genes was calculated to obtain the similarity matrix. In order to make the network conform to the scale-free network distribution, the appropriate weight value is selected and calculated. Select the appropriate soft threshold to measure the connectivity between genes. The adjacency matrix is transformed into a topological overlap matrix (TOM); At the same time, WGCNA package is used to cluster the matrix. For the generated cluster tree, the dynamic tree cutting method is used to cut the gene cluster tree. Genes with similar expression patterns were assigned to a branch, and each branch represented a co-expression module. In this study, the soft threshold is 6 and the minimum module size is set to 30. After analyzing the correlation between modules and clinical features, we selected the hub module which is most relevant to multiple myeloma. Gene significance (GS) represents the association between genes and clinical features. Module members (mm) reveal the relationship between module characteristic genes and gene expression.

Venn Diagram of Hub Genes

The crossing of central genes obtained by Cytohubba and key network modules constructed by WCGNA is considered to be a real hub gene, which is more related to the occurrence and development of multiple myeloma (MM). In order to achieve this intersection, the Lian chuan biological cloud platform (https://www.omicstudio.cn/tool) draw a Venn diagram.

Survival Analysis of Hub Genes

Then, the RNA sequencing data and matched clinical characteristics of multiple myeloma obtained from the Cancer Genome Atlas (TCGA) were used to analyze the impact of different expressions of hub genes on the survival of patients with multiple myeloma (MM).

Verification and Analysis of Hub Genes

Patient and Tissue Samples

(P)–Patient: From October 2020 to October 2021, 19 newly diagnosed multiple myeloma (MM) patients were selected from the Department of Hematology, Zhongshan Hospital Affiliated to Xiamen University, including 7 males and 12 females, with a median age of 61 (42 ~ 71) years. All the selected patients were diagnosed by clinical, myelogram, immunofixation electrophoresis, M protein quantification, whole-body flat bone radiography, bone scan, CT, MRI and other examinations. The diagnosis and staging were based on the Chinese guidelines for the diagnosis and treatment of multiple myeloma.¹⁸

(I)– Intervention: All patients in our hospital were treated with chemotherapy based on lenalidomide, bortezomib and dexamethasone.

(C)– Comparison: 10 healthy donors were selected as control group one. At the same time, 19 newly diagnosed multiple myeloma (MM) patients served as control group two after three months of standardized treatment.

(O)– Outcome: Compare the changes of RNA expression level (RT qPCR).

(T)–Time: The outcome would be measured after standardized treatment for three months.

In this study, 3mL heparin anticoagulation peripheral blood samples were collected before and after treatment. Cells were extracted by centrifugation within 24h and frozen at −80 °C. In this study, the use of patient organization and the acquisition of patient information were approved by the ethics committee of Zhongshan Hospital Affiliated to Xiamen University.(We ensured that all subjects obtained informed consent before the experiment was enrolled. The study was conducted in accordance with the principles of the declaration of Helsinki.)

RT qPCR Detection of RNA Expression Level

RNA extraction and quantitative polymerase chain reaction. RNA isolation of fresh multiple myeloma tissue samples was performed with Trizol reagent (Roche). Complementary DNA (cDNA) was synthesized by reverse transcription using hifiscript Kit (Invitrogen). Quantitative polymerase chain reaction was performed on TNFSF11 (RANKL) and GAPDH using SYBR Green qPCR master kit (Roche). GAPDH was selected as the internal reference gene. The primer sequence is as follows: STMN2 forward 5- tctcttgctgttgctgatgg-3; STMN2 reverse 5-ctcgtgagaccttcgctctt-3; GAPDH forward 5-CAAGGTCATCCATGACAACTTTG-3; GAPDH reverse 5- GTCCACCACCCTGTTGCTGTAG-3.

Results

Identification and Enrichment of DEGs in Multiple Myeloma

The GEO2R tool is used to analyze the DEGs in GSE146649. We screened the 424 DEGs between multiple myeloma and normal tissues, and found 350 up-regulated genes and 74 down regulated genes. DEGs are displayed in the volcano map and heat map (Figures 1 and 2). In GO analysis (Figure 3), biological process (BP) is mainly enriched in transcription regulation, cell adhesion and cell differentiation of RNA polymerase II promoter. Cellular components (CC) are mainly extracellular exosomes, cytoplasm, cytoplasmic membrane, etc. Molecular function (MF) is significantly enriched in protein binding rate. Further analysis of KEGG pathway enrichment of differentially expressed genes showed that the differentially expressed genes were mainly enriched in PI3K Akt signaling pathway, actin cytoskeleton regulation, AGE-RAGE signaling pathway in diabetes complications, HIF-1 signaling pathway, TGF-β Signaling pathway, Rap1 signaling pathway, and tumor necrosis factor signaling pathway. (Figure 4).

Figure 1 Volcano map of DEGs between MM case group and normal control group. The red plots in the volcano represent up-regulation and the green points represent down-regulation.

Figure 2 Heatmap of the DEGs according to the value of |log FC|. The color in heat maps from blue to red shows the progression from low expression to high expression cell log. Abbreviation: FC, log fold change.

Figure 3 GO term enrichment analysis of module genes.

Figure 4 Bar chart of KEGG enrichment. The horizontal axis represented ratio of DEGs under the pathway to total number of DEGs. The vertical axis was the description information of enriched pathways. The color of points corresponded to the value of P. adjust, and the length of bar corresponded to the number of DEGs under the pathways.

Selection of Key Modules

Through string analysis, 424 DEGs are input into the PPI network, including 302 nodes and 1165 edges (Figure 5). Cytohubba is a plug-in that ranks nodes according to their network capabilities.¹⁶ It provides 11 analysis methods, including maximum neighborhood component density, maximum Group Centrality, maximum neighborhood component, degree, edge penetration component and six centralities (bottleneck, heart rate deviation, proximity, radius, medium and stress). This study used all 11 analysis methods in the PPI network to identify the top 15 genes (Figure 6). According to 11 analysis methods of Cytohubba, 48 genes were screened. This finding may indicate that 48 genes play an important role in multiple myeloma.

Figure 5 424 DEGs were filtered into the DEGs PPI network complex that contained 302 nodes and 1165 sides.

Figure 6 (A) Maximal Clique Centrality methods in Cytohubba. (B) Radiality methods in Cytohubba. (C) Betweenness methods in Cytohubba. (D) Bottleneck methods in Cytohubba. (E) Closeness methods in Cytohubba. (F) Degree methods in Cytohubba. (G) Density of maximum neighborhood component methods in Cytohubba. (H) EcCentricity methods in Cytohubba. (I) Edge percolated component methods in Cytohubba. (J) Maximum neighborhood component methods in Cytohubba. (K) Stress methods in cytoHubba.

Construction and Analysis of the Key Gene Co-Expression Network

A total of 40 clinical samples of GSE146649 were collected for analysis. WGCNA performs R language analysis and divides highly related genes into one module choice β= 6 (scale-free r2=0.90) as the appropriate soft threshold to ensure scale-free analysis. Six modules were identified. It was found that the blue module was mainly related to multiple myeloma (Figures 7 and 8). In addition, the relationship between modules and clinical features was also studied for gene significance (GS) and module members (MM) in the six modules. The heat map of the phenotypic correlation between modules and clinical features is shown in the figure (Figure 9). We can see that the blue module has the most significant correlation with multiple myeloma (diagnosis), which is 0.94, and the p value is statistically significant. Later, the scatter diagram of the correlation between gene significance (GS) and module members (MM) (Figures 10 and 11) was carried out. We can see that the blue module has the most significant correlation with multiple myeloma (diagnosis), while the blue module has the least significant correlation with health. Finally, according to the critical values: module member value >0.80 and GS value >0.2, the gene in the blue module is selected as the central gene. A total of 295 central genes from the blue module were selected for further analysis.

Figure 7 Hierarchical cluster analysis was conducted to detect co-expression clusters with corresponding color assignments. Each color represents a module in the constructed gene co-expression network by WGCNA.

Figure 8 Heatmap plot of the adjacencies of modules. Red represents high adjacency (positive correlation) and blue represents low adjacency (negative correlation).

Figure 9 Module–trait associations. Each row corresponds to a module, and each column corresponds to a trait. Each cell contains the corresponding correlation and P value. The table is color-coded by correlation according to the color legend.

Figure 10 Scatter diagram of blue module and initially diagnosis.

Figure 11 Scatter diagram of blue module and healthy.

Identification of Key Genes

295 genes constructed by WGCNA were selected as central genes. In the PPI network, 48 genes identified by Cytohubba were defined as hub genes. Finally, the intersection of Cytohubba and WGCNA is four genes (TNFSF11, FGF2, SGMS2, IGFBP7), which are considered the “real” hub genes (Figure 12).

Figure 12 Venn diagram of hub genes.

Survival Analysis of Hub Gene

Next, the survival condition of four hub genes were analyzed. Survival is affected by many factors, including treatment strategies, patients’ physical status, tumor characteristics and so on. The markers related to these factors are often the key to the prognosis evaluation of patients with multiple myeloma. Targeted therapy against these markers will become a breakthrough in the treatment of multiple myeloma in the future. In this study, the TCGA database was used to analyze the relationship between key genes and prognosis. It was found (Figure 13) that the long-term overall survival rate of multiple myeloma with TNFSF11, FGF2, SGMS2 and IGFBP7 mutations was significantly lower, and the difference was statistically significant; This indicates that the differential expression of target genes has adverse effects on the survival and overall survival of patients with multiple myeloma. Therefore, the screening of four target genes (TNFSF11, FGF2, SGMS2, IGFBP7) is the key molecule for the occurrence and development of multiple myeloma, and can be used as an important target for the diagnosis and treatment of multiple myeloma.

Figure 13 Survival analysis of TNFSF11, FGF2, SGMS2, IGFBP7. The red represented the MM group in which the key gene was highly expressed, and green represented the normal group in which the target gene was lowly expressed. P <0.05 was considered statistically significant.

Verification of Key Genes

qPCR Detection

Multiple myeloma whole blood samples and normal whole blood samples were evaluated by RT-qPCR to verify whether TNFSF11 was differentially expressed in the samples.

The expression of TNFSF11 and internal reference GAPDH in 29 samples was detected by qPCR. The gene amplification curve and dissolution curve are shown in Figures 14 and 15. The expression level of TNFSF11 in whole blood samples of patients with multiple myeloma was significantly higher than that of healthy donors (Figure 16). After treatment, the expression level of TNFSF11 in whole blood samples of patients with multiple myeloma decreased (Figure 17).

Figure 14 TNFSF11 (RANKL) amplification curve.

Figure 15 TNFSF11 (RANKL) dissolution curve.

Figure 16 The green plots in the Point-fold Line Chart indicates the expression of TNFSF11 (RANKL) in multiple myeloma, and The blue plots in the Point-fold Line Chart indicates the expression of TNFSF11 (RANKL) in healthy (normal).

Figure 17 The green plots in the bar chart represent the expression of TNFSF11 (RANKL) in patients with multiple myeloma before treatment, and the blue plots in the bar chart represent the expression of TNFSF11 (RANKL) in patients with multiple myeloma after treatment.

Discussion

Multiple myeloma is the second most common hematological malignancy after lymphoma. There is no clinical treatment plan for multiple myeloma. The disease is mainly controlled by early chemotherapy, but the treatment effect is not ideal. Therefore, it is very important to screen early diagnostic markers and regulatory pathways related to multiple myeloma. In this study, we used the gene expression profile of GSE146649 from 31 patients with multiple myeloma and 9 normal subjects to find some biomarkers and clarify the molecular mechanism of multiple myeloma.

The study found that 424 DEGs were associated with multiple myeloma, including 350 up-regulated genes and 74 down regulated genes. PPI and WGCNA analyses were used to select central genes associated with multiple myeloma. In addition, the related functions and pathways of multiple myeloma were analyzed. In GO analysis, DEG is mainly enriched in the transcription regulation, cell adhesion, cell differentiation of RNA polymerase II promoter, significant enrichment in protein binding rate, cell adhesion is the key to the spread and progress of multiple myeloma. Research shows that MM plasma cell trafficking and homing to the BM are regulated by soluble factors, mainly chemokines such as CXCL12, as well as by direct cell–cell interactions through adhesion molecules expressed by plasma and bone marrow-resident cells.¹⁹ The final phases in MM involve the egress of myeloma cells from the BM to the bloodstream to colonize different organs, once they become independent from growth signals provided by the BM milieu, a condition named extramedullary disease.^20,21 Extramedullary multiple myeloma remains challenging from a therapeutic and biological perspective.²²

Further analysis of KEGG pathway enrichment of differentially expressed genes showed that the differentially expressed genes were mainly enriched in PI3K Akt signaling pathway, actin cytoskeleton regulation, AGE-RAGE signaling pathway in diabetes complications, HIF-1 signaling pathway and TGF-β Signaling pathway, Rap1 signaling pathway, and tumor necrosis factor signaling pathway. We know that the ability of bone formation depends on the activity of osteoblasts and the osteogenic differentiation ability of bone marrow stromal cells (BMSCs). Therefore, in order to fully understand the pathogenesis of multiple myeloma and explore its characteristics, it is feasible to study from two aspects. Pi3k/akt signaling pathway is the key signaling pathway for osteogenic differentiation of bone marrow mesenchymal stem cells.²³ It can promote chondrocyte proliferation and inhibit autophagy. Relevant studies have confirmed that pi3k/akt signaling pathway is involved in pathological bone diseases such as osteoporosis, osteoarthritis and osteosarcoma,^24,25 and regulates the proliferation, differentiation and apoptosis of osteoclasts and osteoblasts.^26,27 It is reported²⁸ that AGEs modified collagen can promote the release of intracellular reactive oxygen species, interfere with the adhesion between osteoblasts and matrix, and inhibit the differentiation and proliferation of osteoblasts. Rage and its ligands play an important physiological role in bone. The increase of rage and its ligands is obviously related to various bone related diseases.²⁹ HIF-1 signaling pathway plays an important role in bone formation and bone resorption, and is closely related to the promotion of osteoblast bone formation, this signaling pathway can increase the number of osteoclasts and osteoclastic activity.³⁰ TGF- β signaling pathway can make progenitor cells enter adipogenic or osteoblastic lineage by recombining cytoskeleton, so as to improve the osteogenic ability of mesenchymal stem cells.³¹

WGCNA analysis showed that the six modules had related expression patterns. 296 hub genes related to metastasis were selected from the blue module. Based on the plug-in Cytohubba in Cytoscape, four key genes were screened. TNFSF11, FGF2, SGMS2 and IGFBP7. FGF2 is a multifunctional growth factor belonging to the FGF family. FGF2 is an efficient mitogen, which can promote the proliferation of osteoblast precursor cells.³² In addition, FGF2 can inhibit the aging of BMSCs and promote cell proliferation through pi3k/akt-mdm2 pathway, to maintain the stem cell characteristics of BMSCs.³³ Sgms2 is a sphingomyelin synthetase SMS2 that is significantly expressed in cortical bone and encodes on the plasma membrane. Some research results show³⁴ that SMS2 has bone regulatory properties. The mutation of SGMS2 is the basis of a series of bone diseases, including simple osteoporosis and complex bone dysplasia, indicating that plasma membrane combined with sphingomyelin metabolism plays a key role in the balance of the bone environment. IGFBP7 belongs to the insulin-like growth factor binding proteins (IGFBPs) family, which consists of six high affinity and several low affinity members,³⁵ and plays an important role in bone metabolism. It has been reported³⁶ that the role of IGFBP7 involves insulin resistance, which is considered to be closely related to the development of osteoporosis. Ye et al. Found³⁷ that IGFBP7 inhibits osteoclast formation in vitro and in vivo, suggesting that IGFBP7 is a negative regulator of osteoclast formation and plays a protective role in osteoporosis. Zhang et al. Found³⁸ that IGFBP7 increased OPG expression in vivo and enhanced osteogenic ability. We know osteoprotegerin (OPG) and nuclear factor κ B receptor activating factor (RANK) are members of the tumor necrosis factor-α receptor superfamily,³⁹ TNFSF11 (RANKL) is a transmembrane receptor, belonging to the tumor necrosis factor (TNF) superfamily. TNFSF11 (RANKL) is mainly distributed on the surface of osteoblasts and bone marrow stromal cells, and plays an important role in bone formation.⁴⁰

The combination of TNFSF11 (RANKL) and RANK can promote the differentiation of osteoclasts, while OPG can inhibit the differentiation and activity of osteoclasts,⁴¹ delaying the activation of osteoclast precursor cells and maintain the balance of bone metabolism. Therefore, RANK, RANKL (TNFSF11) and OPG are the key regulatory mechanisms to control the differentiation and activity of osteoblasts and osteoclasts.⁴² The ratio of RANKL (TNFSF11)/OPG is an important lever index to evaluate bone health and regulate the balance of bone resorption.⁴³ This study found that the expression of TNFSF11 (RANKL) in multiple myeloma was up-regulated compared with the normal group. TNFSF11 is highly expressed in whole blood samples of patients with multiple myeloma. The expression level of TNFSF11 (RANKL) in peripheral blood of patients with multiple myeloma is closely related to osteolytic osteopathy, which can reflect the short-term bone metabolism.⁴⁴

In this study, TNFSF11, FGF2, SGMS2 and IGFBP7 were selected as important key genes in the pathogenesis of mm by bioinformatics method, and then the levels of TNFSF11 (RANKL) in peripheral blood of patients with multiple myeloma diagnosed in our hospital before and after treatment were analyzed and compared by PCR technology. We found that the expression level of TNFSF11 (RANKL) in peripheral blood of patients with multiple myeloma decreased significantly after treatment, and compared with healthy patients, TNFSF11 (RANKL) was highly expressed in patients with multiple myeloma, suggesting that TNFSF11 (RANKL) may play a pathogenic role in multiple myeloma. In addition, in this study, we used comprehensive analysis from WGCNA and DEGs to screen potential biomarkers, which may also make our results more credible. The biological functions and potential molecular mechanisms of the key genes TNFSF11, FGF2, SGMS2 and IGFBP7 explored in this study still need to be explored in future research.

Conclusion

In conclusion, through comprehensive bioinformatics analysis, TNFSF11, FGF2, SGMS2 and IGFBP7 were screened as key genes related to multiple myeloma, and their roles and pathways in multiple myeloma were discussed. At the same time, RT-qPCR was used to preliminarily confirm that TNFSF11 was highly expressed in whole blood samples of patients with multiple myeloma. Combined with the clinical conditions of patients and current relevant research progress, we believe that the above four key genes are still worth further exploration and analysis. The results of this study increase our understanding of biomarkers and potential therapeutic targets of multiple myeloma. However, we need to further study the role and mechanism of TNFSF11, FGF2, SGMS2 and IGFBP7 in the future development of multiple myeloma.