Identification of a three-gene-based prognostic model in multiple myeloma using bioinformatics analysis

Data preparation

The raw data of GSE6477 and GSE136337 were downloaded from the gene expression omnibus (GEO) website (http://www.ncbi.nlm.nlh.gov/geo/). There were 101 case and 15 normal control in GSE6477, was based on GPL96Platforms ([HG-U133A] Affymetrix Human Genome U133A Array). There were 401 cases with overall survival (OS) and other clinical information was chosen to identify prognostic genes and construct the prognostic model in GSE136337. All the patients were randomly divided into 2 groups (60% as the training set and 40% as the test set). The Robust multi-array average (RMA) method was used to normalize data.

Analysis of differentially expressed genes (DEGs)

Limma package was used to identify differentially expressed genes (DEGs) (Ritchie et al., 2015), and the threshod value was |log2foldchange| >1.0 and adj.P.Val <0.05.

GO term and KEGG pathway enrichment analysis

To better explore the biological significance of DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed (Yu et al., 2012). q value <0.05 was used as the cutoff value.

A protein-protein interaction (PPI) network

A PPI network of these DEGs was constructed by using the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/) (Szklarczyk et al., 2015).

Co-expression network analysis

The co-expression network was constructed by WGCNA to explore the correlation of genes (Top50%mad) and search for important interacted gene modules. PickSoftThreshold function of the WGCNA was used to calculate the soft thresholding power β value, and the β value was set as 6. Then, we constructed the co-expression modules using WGCNA algorithms with R software package. Next, Gene significance (GS) and Module membership measures were used to identify genes with high clinical significance as well as high Module membership. Module membership was the correlation between gene expression profile and module eigengene. While, GS was the absolute value of the correlation between gene expression and module traits.

Identification of prognostic key genes and construction of the prognostic model

Key genes were identified both from hub genes and the DEGs (Xia et al., 2018). All the patients were randomly divided into two groups (60% as the training set and 40% as the test set) and univariate and multivariate Cox regression analyses were used to investigate the correlation between OS and the expression level of key genes. P-value <0.05 was considered of significant difference in the univariate analysis. Key genes with P < 0.05 based on the univariate analysis were further included in the multivariate Cox regression analysis. Finally, a prognostic gene signature of MM patients was established based on a linear combination of the regression coefficient derived from the multivariate Cox regression model (β) multiplied with its expression level. The patients were separated into low- and high-risk groups based on the optimal cut-off value (Camp, Dolled-Filhart & Rimm, 2004) Kaplan–Meier (K–M) survival curves, time-dependent receiver operating characteristic (ROC) curve analyses and area under the curve (AUC) were used to evaluate the prediction accuracy of our model (Heagerty & Zheng, 2005). We also validated the prognostic model in the test set and whole set by K–M curves, ROC curve analyses, and AUC.

Identification of independent prognostic parameters of MM

To investigate independent prognostic parameters and to identify whether the prognostic model was the independent prognostic parameter, univariate and multivariate Cox regression analyses were performed using the Cox regression model method with forwarding stepwise procedure. Other clinical parameters including age, gender, number of transplantation (NOT), tumor staging (ISS and R-ISS), cytogenetic abnormalities, lactate dehydrogenase (LDH) level, albumin level, β2-microglobulin level, and percentage of abnormal plasma cell in bone marrow. P-value <0.05 was considered of significant difference. Also, we compared the AUC of our model with other conventional models such as age, gender, ISS, and R-ISS models.

Identification of DEGs in MM and the enrichment of these DEGs

The DEGs of GSE6477 were screened by the limma package. A total of 178 DEGs were identified, including 59 upregulated genes and 119 downregulated genes in MM case samples compared to normal control samples (Fig. 1A). The DEGs were shown on a heatmap (Fig. 1B). GO analysis showed that these genes were significantly enriched in neutrophil activation, humoral immune response, and structural constituent of ribosome (Fig. 1C). In KEGG analysis, the DEGs were mainly enriched in Asthma, Staphylococcus aureus infection, and Ribosome (Fig. 1D). Next, A PPI network, as demonstrated in Fig. 2, was constructed using the online STRING website to understand the interaction among the related DEGs in MM.

Co-expression analysis and hub genes identification

WGCNA was used to construct co-expression network and to find the modules of highly correlated genes. Firstly, we calculated the soft thresholding power β, and the power of β was set at 6 to ensure a scale-free network in our research (Figs. 3A–3D). Next, we identified 15 modules for further analysis (Fig. 4A). The heatmap and meta-modules were used to visualize the gene network (Figs. 4B and 4C). Among these 15 modules, module black and module purple were the most relevant modules with cancer traits (Fig. 4D). Subsequently, intramodular analysis of the genes in the 15 modules was followed. And genes in the black and purple modules were found to have high correlations with MM (Figs. 4E and 4F). Finally, genes in these two modules were selected as hub genes according to the cut-off criteria: module membership values >0.80 and GS values >0.2. Totally, there were 92 hub genes from black and purple modules that were chosen for next analysis.

Identification of survival related key genes and construction of three-gene prognostic signature

Key genes were identified both from hub genes and the DEGs. Totally, 47 key genes were identified (Fig. 5). To evaluate the prognostic values of these key genes in MM patients, we next performed a univariate Cox regression, and identified 8 genes significantly related to OS in MM patients. Then, a stepwise multivariate Cox regression analysis was performed, and three genes (LYVE1, RNASE1, and RNASE2) were eventually selected to construct a prognostic model (Fig. 6). Risk score = (−0.59714 * expression level of LYVE1) + (−1.35971* expression level of RNASE1) + (0.28240 * expression level of RNASE2). Next, we calculated the optimal cut-off value for the three-gene expression risk score with X-Tile software, and divided patients in training set (n = 241) into high risk group and low risk group. The K–M curve showed that high risk group had poorer prognosis than low risk group (P < 0.0001) (Fig. 7). The AUCs of time-dependent ROC curves were calculated to assess the prognostic capacity of the three-gene signature. The AUCs of risk scores in training set at 8-,10-, and 12-year survival times were 0.717, 0. 720 and 0.800 (Fig. 8A).

To verify the predictive value of the three-gene signature, we used the internal test set (n = 160) and whole set (n = 401) to assess the results from the training set. Consistent of the findings of the training set, the K–M survival curves of test set and whole set revealed that worse OS existed in the high risk groups (Fig. 7). The AUCs of risk scores in test set at 8-,10-, and 12-year survival times were 0.709, 0.768 and 0.766 (Fig. 8B), and the whole set were 0.714, 0.741 and 0.783, respectively (Fig. 8C). We also compared the prognostic capacity of our model with other conventional models, including age, gender, ISS, and R-ISS stage in whole set. The AUCs for 8-, 10-, and 12-year OS predictions for age were 0.500, 0.518, and 0.478; and for gender were 0.516, 0.506, and 0.450 (Figs. 8D–8F). The AUCs of ISS stage at 8-, 10-, and 12-year survival times were 0.503, 0.547, and 0.633; and of R-ISS stage were 0.613, 0.616, and 0.634 (Figs. 8D–8F). Our model had the highest AUCs, indicating a superior prognostic value. Together, these results demonstrated that our model was capable of predicting OS of MM patients.

The prognostic model is independent for MM patients

The univariate and multivariate Cox regression analyses were utilized to investigate independent prognostic parameters and identify whether the three-gene prognostic model was the independent prognostic parameter in 401 MM patients with complete clinical information from GSE136337. Univariate Cox regression analysis showed that NOT, tumor stage (ISS and R-ISS), t(11,14), albumin level, β2-microglobulin level, and the prognostic model had prognostic values (Fig. 9A). After the multivariate Cox regression analysis, NOT, R-ISS, t(11,14), LDH level and the prognostic model were proved to be independent prognostic factors for OS of MM patients (Fig. 9B).

Correlation between the prognostic model with clinical characters

Correlation between risk score and various clinical parameters were shown in Fig. 10. Our results indicated that risk score was significantly correlated with ISS, and R-ISS stage; patients in phase III of ISS or R-ISS had the highest risk score when compared with patients in phase I and II (Fig. 10). However, no differences were found between risk score and age or gender or NTO (Table 1). Besides, the relationships between the three genes in the prognostic model with clinical characters were also analyzed. We found that as LYVE1 expression, which was higher in males than in females, decreased, ISS stage in MM patients increased significantly (Fig. 10A, Table 1). Additionally, as LYVE1 and RNASE1 expression decreased, R-ISS stage in MM patients increased (Figs. 10B and 10C). These findings showed that LYVE1 and RNASE1 exerted protective effects against disease progression.