Identification of key genes and immune infiltration in multiple myeloma by bioinformatics analysis

ABSTRACT

Objective:

Multiple Myeloma (MM) is a hematologic malignant disease with unclear molecular mechanisms. This integrated bioinformatic study aimed to identify key genes, pathways and immune cell infiltration pattern in MM.

Methods:

Differentially expressed genes (DEGs) from GSE6477 and GSE16558 dataset were filtrated with R package ‘limma’, whose function were explored by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The key genes were selected from Protein–protein interaction network (PPI) and logistic regression model. The correlation between key genes and survival in MM was evaluated using the survival and survminer package. Additionally, immune filtration analysis was accomplished by CIBERSORT tools.

Results:

118 DEGs (92 up-regulated and 26 down-regulated) from two GSE datasets were identified, which were closely related with B cell receptor signaling pathway and Epstein-Barr virus infection. Furthermore, CD24 and PTPRC of five hub genes identified in PPI network were further screened out by the logistic regression model. Besides, CD24 and PTPRC expression were significantly correlated to the survival time in MM patients. Finally, MM might cause different infiltrating immune cell compositions, including increased infiltrations of B cells memory, Plasma cells, T cells CD4 memory resting, T cells follicular helper, Tregs, NK cells resting, Macrophages(M0/M1), Dendritic cells resting and Mast cells activating, and lower proportions of B cells naïve, T cells CD4 naïve, Macrophages M2 and Neutrophils.

Conclusion:

Targeting CD24 and PTPRC as molecular markers of MM is valuable to MM therapy. Moreover, the immune cell infiltration will provide new insights into MM immunopathology.

KEYWORDS:

• Multiple myeloma
• bioinformatics analysis
• immune filtration
• differential expressed genes
• CD24
• PTPRC
• logistic regression
• survival prognosis

1. Introduction

Multiple Myeloma (MM) is a type of hematologic tumor that affects plasma cells and presents abnormal plasma cell infiltration in the bone marrow and over-production of monoclonal protein in the serum and/or urine [1]. It accounts for approximately 1% of neoplastic diseases and 13% of hematological diseases [2]. The neoplastic cell proliferation is promoted by the cumulative effect of genetic and epigenetic modifications on plasma cells together with the involvement of the bone marrow microenvironment [3]. The typical symptoms of MM include hypercalcemia, renal failure, anemia and bone lesions [4]. Although the development of novel chemotherapeutic agents and cell therapy has improved patient outcomes to some extent, MM is still considered incurable 2 to >10 years median survival time [5]. Consequently, uncovering the pathogenesis and identifying novel biomarkers of MM are crucial for early diagnosis, prevention, and targeted therapy.

In order to elucidate the etiopathogenesis of MM, many studies have been conducted via gene expression profiling based on high-throughput sequencing technology [6,7]. Related research data are available in public databases such as the Gene Expression Omnibus (GEO) database [8], but inconsistent results and contradictory conclusions have often been observed due to the different study methods and sample selection. Therefore, a better approach combining expression profiling techniques with the integrated bioinformatics method might eliminate the inaccuracy and poor reproducibility in MM-related research. Furthermore, it is possible to identify novel diagnosing biomarkers and therapeutic targets for MM.

In this study, the differentially expressed genes (DEGs) between patients with MM and healthy controls were identified by using integrated bioinformatics analysis on data obtained from the GEO database. The messenger RNA (mRNA) gene expression profiles of MM including GSE6477 and GSE16558 have been selected and the common DEGs have been screened using the R package with the standard data processing. Then, gene ontology (GO) and KEGG pathway enrichment analysis, along with the protein–protein interactions (PPI) network establishment integrating modular analysis was performed to identify hub genes in MM. In addition, CIBERSORT was used to conduct the immune infiltrates. Our findings may highlight some initial signs for the diagnosis and treatment of MM.

2. Materials and methods

2.1. Microarray data information and DEG identification

GSE6477 and GSE16558 were selected and downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo). The GSE6477 dataset, which included purified plasma cell samples from 125 diagnosed patients with MM and 15 healthy individuals, was performed on GPL96 platforms (HG-U133A Affymetrix Human Genome U133A Array). The GSE16558 data with 65 MM patients and 5 healthy controls of purified plasma cell samples was obtained from GPL6244 platforms (HuGene-1_0-st Affymetrix Human Gene 1.0 ST Array). The R software was utilized to screen DEGs across experimental conditions between MM patients and healthy adults. In the present study, the p-value < .05 and [log2FoldChange (log2FC)] > 1 were the criteria set for screening. The volcanic maps were used to visually show the genes identified and the top 50 remarkable varied genes with statistical significance identified in a heatmap drawn by R package.

2.2. GO and pathway enrichment analysis

Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were used to annotate gene functions and signaling pathways involved. The GO enrichment represents the protein functions in the three hierarchies including Molecular Functions (MF), Biological process (BP), and Cellular Component (CC). The items with p-values < .05 were considered as significantly enriched.

2.3. PPI network establishment and modular analysis

STRING database was applied for the establishment of protein–protein interaction (PPI) network, providing an insight into the functional relationships of interacting genes or proteins. The data for all proteinprotein interactions is summarized, analyzed and calculated. The PPI network was visualized by Cytoscape software 3.7.2. Moreover, Molecular Complex Detection (MCODE), an application in Cytoscape for finding clusters (highly interconnected regions) in a PPI network, was utilized to identify the most significant module (based on the MCODE score and node number).

2.4. Gene set enrichment analysis

GSEA is a computational way to interpret the expression data of one gene set and evaluate the significant difference in various biological states. It was conducted for enrichment analysis. The GSEA was performed with the following settings: (1) default weighted enrichment method, (2) p < .05, (3) false discovery rate(FDR)q < 0.25. The enrichment results were visualized by using the plugin in Cytoscape.

2.5. Logistic regression model

The logistic regression model is a relationship between predictor variables and a categorical response variable. In the present study, all samples were classified into two sets- MM samples and the normal control samples. Subsequently, the logistic regression analysis was adopted to prioritize genes followed by multivariate logistic regression analysis to identify the target genes contained in the model by using the Generalized Linear Models (GLM) function of R language. The expression of mRNA was the continuous variable, and the sample set was the binary categorical response variable.

2.6. Survival analysis

The correlation between the selected gene expression and survival in MM patients was analyzed using the ‘survival’ and ‘survminer’ packages. Based on the GSE4581 data set (n = 414), the patients were classified into two groups according to the median of each selected gene expression for overall survival.

2.7. Immune infiltration analysis

Immune cell infiltration characteristics of MM were assessed by uploading gene expression matrices to the Cell Type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) website (https://cibersort.stanford.edu/), and the immune cell infiltration matrix was obtained. The ‘corrplot’ package was used to visualize the interaction infiltrating immune cells of 22 different immune types, and a p-value less than .05 was used as a criterion to filter qualified samples. The difference in immune infiltration extents represented by immune infiltration scores among subgroups was compared by the Wilcoxon test.

3. Results

3.1. Identification of DEGs

The gene expression profiles were obtained based on plasma cell samples from patients of multiple myeloma and healthy controls from the GSE6477 and GSE16558 datasets. We applied ‘limma’ R package to identify DEGs associated with MM and set p < .05 and [logFC] ≥ 1 as the cut-off value (Figure 1(A)). The top 50 genes with the most significant variabilities for each dataset between MM patients and healthy controls were clustered in the heatmap separately (Figure 1(B)). The following analysis on GSE6477 identified 979 DEGs in MM patients compared with the healthy donors, including 609 upregulated ones and 370 downregulated ones. Moreover, 903 DEGs were identified including 730 upregulated and 173 downregulated genes in GSE16558 datasets. Finally, we obtained 118 common DEGs from two datasets including 92 upregulated and 26 downregulated genes (Figure 1(C)).

Figure 1. Identification of DEGs in GSE6477 and GSE16558. (A) Volcano plots of the distribution of DEGs in two datasets. (B) Heatmap of the top 50 DEGs in two datasets. (C) Commonly upregulated or downregulated DEGs in two datasets.

3.2. GO and KEGG enrichment analysis

To further explore the potential biological function of the 118 DEGs summarized above, we performed GO analysis and KEGG pathway enrichment analysis. The GO terms for DEGs are shown in Figure 2(A,B). DEGs in biological process (BP) were significantly associated with lymphocyte proliferation, mononuclear proliferation, leukocyte proliferation, B cell activation and B cell differentiation. For molecular function (MF) GO items, DEGs were mostly observed in the proteoglycan binding, heparan sulfate proteoglycan binding, growth factor binding and integrin binding. The cellular components (CC) analysis showed that DEGs were mainly enriched in the external side of the plasma membrane, secretory granule lumen, cytoplasmic vesicle lumen, vesicle lumen and focal adhesion. The top KEGG pathway showed that the related DEGs were significantly enriched in the B cell receptor signaling pathway, Epstein–Barr virus infection and coronavirus disease pathway (Figure 2(C,D)).

Figure 2. Functional enrichment analysis. (A) The GO analysis of DEGs, including biological process (BP), molecular function (MF) and cellular component (CC). (B) The chord diagram of enriched GO terms. (C) The bubble map of the KEGG enrichment analysis of DEGs. (D) The chord diagram of enriched KEGG terms.

3.3. GSEA analysis

In order to identify the possible mechanism in MM, GSEA was carried out to screen out significant differential biological events between MM and control cases. The analysis indicated that the most significant-enriched gene sets included the NOD-like receptor signaling pathway, the chemokine signaling pathway, the B cell receptor signaling pathway, the oxidative phosphorylation, and the Ribosome and P53 signaling pathway (Figure 3(A–F)).

Figure 3. GSEA showing the most enriched gene sets: (A) NOD-like receptor signaling pathway. (B) Chemokine signaling pathway. (C) B cell receptor signaling pathway. (D) oxidative phosphorylation. (E) Ribosome (F) P53 signaling pathway.

3.4. PPI network analysis and hub gene selection

To predict the functional connections and interactions of proteins encoded by DEGs in MM, we constructed the protein–protein interaction (PPI) network of DEGs based on the STRING database and Cytoscape software. After removing the separated nodes, we described a complex PPI network of DEGs shown in Figure 4(A). Subsequently, we identified three significant modules by utilizing cluster analysis of the PPI network in Cytotype MCODE (Figure 4(B–D)). Additionally, key genes of the PPI network were screened by using the cytoHubba. Hub genes selected contained PTPRC, CD79A, CD24, CR2 and SYK, and their interactive network is shown in Figure 4(E).

Figure 4. PPI network and modular analysis: (A) PPI network of DEGs. (B–D) Three modules generated from PPI by the MCODE. (E) Top five hub genes from the PPI network.

3.5. Logistic regression model for prediction of MM

In order to confirm the diagnostic role of the above key genes in MM, logistic regression models were constructed in GSE6477 datasets based on these five mRNAs separately. In the model, mRNAs were employed as the independent variables and the model followed the normal distribution. Meanwhile, the variables including mRNA of PTPRC and CD24 in the model had a significant linear relationship with the response variables, indicating that these two independent variables were fit for use in the model (Figure 5(A,B)). Furthermore, the reliability and validity of this model were evaluated by ROC analysis. The AUC value of the model in the GSE6477 dataset was 0.9635, showing the good explanatory power of the model (Figure 5(C)). Similarly, the AUC value of the GSE16558 dataset was 0.9133 (Figure 5(D)). Overall, the logistic regression model involving the two mRNAs provided the possibility of differentiating MM sample type from healthy controls, and PTPRC and CD24 were potential targets for further research on MM.

Figure 5. The construction of logistic regression model: (A) Component plus residual plot of PTPRC in the model. (B) Component plus residual plot of CD24 in the model. (C) The ROC curve of predicted outcomes of GSE6477 by logistic regression model. (D) The ROC curve of predicted outcomes of GSE16558 by logistic regression model.

3.6. Survival analysis

We explored the potential prognostic value of CD24 and PTPRC for patients with MM. Based on the GSE4581 data set (n = 414), the correlation between gene expression of CD24 and PTPRC and survival time was analyzed by survival and survminer package. The survival analysis showed that there was a significant correlation between the expression of CD24 (p < .01) and PTPRC (p < .01) and survival time in MM patients (Figure 6(A,B)).

Figure 6. Correlation between the expression of CD24 and PTPRC and MM survival. (A) Overall survival analysis of CD24 for MM patients. (B) Overall survival analysis of PTPRC for MM patients.

3.7. Immune cell infiltration results

To investigate the discrepancy in immune infiltration between 22 immune cell subsets in the MM vs. NC group, the CIBERSORT algorithm was employed to analyze the ratios of the infiltrating immune cells in GSE6477 datasets. Then, Immune filtration analysis revealed that MM samples had increased infiltration of B cells memory (p < .0001), Plasma cells (p < .0001), T cells CD4 memory resting (p < .01), T cells follicular helper (p < .1), T cells regulatory (Tregs) (p < .0001), NK cells resting (p < .01), Macrophages M0 (p < .001), Macrophages M1 (p < .0001), Dendritic cells resting (p < .1) and Mast cells activated (p < .1), and lower proportions of B cells naïve (p < .0001), T cells CD4 naïve (p < .01), Macrophages M2 (p < .1) and Neutrophils (p < .0001), which exhibited different infiltrating immune cell composition and may be related to the occurrence of MM (Figure 7(A,B)). As a consequence, it was suggested that the change in the distributing ratios of immune cells was associated with the occurrence and development of MM, which will provide new insights into the immunopathology of MM.

Figure 7. Immune cells filtration in MM samples and heathy controls: (A) The fractions of twenty-two immune cell types in MM samples and healthy controls. (B) Violin diagram of comparisons between 22 types of immune cells in MM samples and healthy controls. The threshold for p-value is .05 for significance. *: p < .05, **: p < .01, ***: p < .001, ****: p < .0001.

4. Discussion

MM is a common hematological malignancy owing to complex outcomes of the tumor development and progression affected by genetic factors [9]. In spite of the progress achieved in the knowledge and new treatments of MM [10,11], the pathogenesis of MM remains unrevealed and this disease is still incurable for most patients. Thus, it is urgent and indispensable to investigate the potential mechanisms of MM at the molecule level. Recently, the rapid development of microarray technology has been extensively adopted to examine the general genetic alterations of MM and identify the possible MM disease predictors [12].

In the present research, the expression of genes in purified plasma cells of MM patients and donors from two microarrays based on GSE6477 and GSE16558 datasets were analyzed. In total, 118 common DEGs were identified between patients and healthy control groups, consisting of 92 upregulated genes and 26 downregulated genes with p < .05 and [logFC] ≥ 1. Thereafter, bioinformatic analysis was used to further annotate the DEGs identified, including GO terms enrichment, KEGG signaling pathway enrichment and PPI network.

GO enrichment analysis in BPs indicated that the functional enrichment of DEGs is mainly involved in lymphocyte proliferation, mononuclear cell proliferation, leukocyte proliferation, B cell activation and differentiation. This finding validated that immune cell dysfunction played a crucial part in the pathogenesis of MM. KEGG analysis showed enrichment of B cell receptor signaling pathway, Epstein–Barr virus infection, Coronavirus disease - COVID-19, Hematopoietic cell lineage and NF-kappa B signaling pathway. Consistent with our results, B cell receptor signaling pathway such as B-cell maturation antigen (BCMA) signaling has been identified to be the most significant pathway in other studies of MM [13]. In addition, the Epstein–Barr virus infection has been reported to be strongly associated with the clinical characteristics and prognosis of MM patients, and MM cells infected by Epstein–Barr virus have been found to produce drug resistance-related proteins [14,15]. Furthermore, we built a PPI network of DEGs and identified 5 candidate hub genes, PTPRC, CD79A, CD24, CR2 and SYK. Subsequently, the logistic regression models based on these hub genes suggested that CD24 and PTPRC had an outstanding performance in the diagnosis of MM among the five hub genes. These results implied the important roles of CD24 and PTPRC in the MM progression.

CD24 is a heavily glycosylated protein localized in membrane lipid raft domains [16], which is mainly expressed by immune cells and often overexpressed in several types of tumors and appears oncogenic [17–19]. In cancer, CD24 plays an important role in regulating the cell proliferation, clonal and migration [20]. Recently, CD24 was identified as a signal of antiphagocytisis which enables the prevention of cancer cells being attacked by binding to Siglec-10-expressing on macrophages [21,22]. It was reported that blocking of CD24–Siglec-10 signaling by monoclonal antibodies strongly intensified the elimination of CD24-positive tumor cells [21], suggesting that CD24 is a promising target for cancer immunotherapy. What’s more, Mensah and colleagues have investigated the connection between the expression of CD24 and the differentiation of B cells, the results have shown that CD24 expression in B cells might be associated with the energy metabolism during all stages of differentiation by the metabolic pathway involving phosphorylation of AMPK [23]. In our study, CD24 expression was significantly increased in MM patients compared with healthy controls from two microarray datasets. In consistent with our results, the findings of Alaterre and colleagues have shown that CD24 was significantly overexpressed in malignant plasma [24]. Combined with the previous studies, CD24 might be a valuable therapeutic target in patients with MM.

Protein tyrosine phosphatase receptor type C (PTPRC), also known as CD45, is a leucocyte common antigen and is expressed in almost all hematopoietic cells excluding mature erythrocytes [25]. The deficiency or abnormal expression of PTPRC is related to various diseases including hematologic tumor [26–28]. For example, Keeney et al found that adult Acute lymphoblastic leukemia (ALL) patients with elevated PTPRC had a reduced rate of survival [29]. A previous study by Saint-Paul has revealed that CD45 is involved in the development of acute myeloid leukemia (AML) and targeting of CD45 positioning among lipid rafts provides a novel strategy in the treatment of AML [30]. Our study results also showed that PTPRC expression was upregulated in MM patients compared with healthy donors from two microarray datasets. The survival analysis showed that high expression of PTPRC was related to worse overall survival in MM patients. Additionally, the research of Kumar et al has demonstrated that the expression levels of PTPRC on the myeloma cells have important clinical and biological effects [31,32]. Another study showing the RNA seq data on purified MM cells from 767 diagnostic patient samples in the online database MM Research Foundation supports that the patients with the reduced PTPRC expression have longer overall and progression-free survival [33].

Taken together, CD24 and PTPRC may participate in the pathogenesis of MM and be considered as possible diagnosis and prognostic markers for MM. Moreover, the study focusing on these two candidate genes might provide a new therapeutic strategy for MM.

Along with in-depth research on MM, the bone marrow microenvironment is believed to play a crucial role in the disease progression, which provides the basis for immunotherapy of MM. Immune cells are important components of this microenvironment by interacting with bone marrow stroma and the myeloma cells. In order to understand the role of immune cell infiltration in MM, CIBERSORT was utilized to conduct an analysis on the MM immune infiltration in statistics. The immune cell composition was significantly different between the samples of MM patients and healthy donors. The overall proportion of NK cell, M0 macrophages and M1 macrophages were strikingly increased in MM patients, compared to that in healthy controls. Alfarra et al have described the key function of enhancing NK cell activity in MM and have anticipated the possibility of blocking NK cell inhibitory receptors as a potential therapeutic strategy for MM treatment [34]. Moreover, it has been found that the higher level of myeloma-associated macrophages correlates with poor prognosis in MM [35]. Interestingly, Zheng et al have indicated that myeloma-associated macrophages prevent myeloma cells from apoptosis [36], which could provide evidence to support our previous prediction that myeloma cells may escape immune surveillance from macrophage through CD24–Siglec-10 signal (‘do not eat me’ signal). Out immune filtration analysis results showed that changes in the proportion of T cells had a critical effect on MM pathology, in which four types of T cells revealed to obviously change in MM patient groups. As reported by the previous study, myeloma cells might express tumor antigens and further initiate the T cell responses. However, T cells promptly collected from patients with MM unsuccessfully recognize autologous tumor cells to secrete IFN-a, indicating the failure of playing the antimyeloma role in vivo [37]. Additionally, a previous study also demonstrated that CD25 + T cells from MM patients were unable to inhibit the proliferation of PBMC under the condition of anti-CD3 stimulation and deduced the dysfunction of MM Tregs [38]. All the study results including the previous literature and analytical data, when taken together, have manifested that immune cell infiltration is of great significance and value in MM and should be focused on in further research in future.

In summary, critical genes for MM were identified and the role of immune cell infiltration in MM was explored in our study. However, there were many limitations in our study: (1) The key genes were predicted by PPI and logistic regression models, lacking experimental validation, and the exact mechanisms of key genes in MM required to be further investigated. (2) The accuracy of CIBERSORT analysis on the basis of insufficient genetic data is probably impacted by the heterotypic interactions of cells, phenotypic plasticity or disease-induced disorders. However, this study provides new insights into the molecular mechanism and immunopathology of MM.

4.1. Conclusions

In this study, based on bioinformatical analyses, we identified key genes and significant pathways from DEGs in MM compared with healthy controls. GO and KEGG enrichment analyses and GSEA were performed to help understand relevant gene functions and mechanisms. CD24 and PTPRC were identified as key genes in MM according to the PPI and logistic regression model. Immune filtration analysis was simultaneously applied utilizing CIBERSORT and differences in immune cell composition were found including B cells, T cells, macrophages, NK cells resting, Dendritic cells resting and Neutrophils. More related studies are awaiting to improve the understanding of their exact roles in MM pathology and to explore new treatment targets.

Acknowledgements

We thank HPS Gene Technology Co., Ltd. for their excellent support in the article. All authors contributed to the study’s conception and design. F.X., L.L. and LM.J. contributed to the data collection, analysis and visualization. F.X. and J.Z. performed the original drafting, writing-reviewing and editing. All authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Ningxia Provincial Natural Science Foundation [grant number 2021AAC05017].