Hub genes and associated drugs for multiple myeloma with 1q21+: identified by bioinformatic analysis

ABSTRACT

While 1q21+ was common genetic alteration and found to have adverse effect on prognosis, the underlying genes remain unclear. Identification of related genes may provide additional help for rational intervention. The microarray dataset GSE2658 associated with MM was downloaded from the Gene Expression Omnibus (GEO) database. The differentially expressed genes (DEGs) were obtained, and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to annotate their functions. The hub genes were derived from the combined results of up-regulated DEGs and weighted gene coexpression network analysis (WGCNA). The receiver operating characteristic (ROC) curves of hub genes were plotted to evaluate correlation with 1q21+. Survival analysis and drug-gene interaction of hub genes were performed separately to find the prognostic value and potential targeted drugs. A total of 55 DEGs were identified. GO and KEGG pathway analyses suggested that the DEGs were related to several pathways of cell proliferation. NVL, IL6R, DUSP23 were proven to be highly correlated with 1q21+ and have adverse effects on prognosis. IL6R, DUSP23 were matched to known interaction-drug. This study revealed potential roles of hub genes in the pathogenesis and progression of MM with 1q21+, further investigations are needed to elucidate the mechanisms.

KEYWORDS:

• Multiple myeloma
• 1q21+
• hub genes
• bioinformatic analysis

Introduction

Multiple myeloma (MM) is the second most prevalent hematological malignancy and characterized by marked biological and genetic heterogeneity [1]. Additional copies of chromosome 1q (1q21+) are defined as a gain of 1q21 (3 copies) and amplification of 1q21 (≥four copies) which can be detected by fluorescence in situ hybridization (FISH). This is 1q21+ was an important genetic alteration in MM, approximately 40% of newly diagnosed MM (NDMM) were identified with 1q21+ [2]. 1q21+ were also associated with clinical features, with the disease progression from smoldering myeloma to newly diagnosed MM to relapsed/refractory disease, there was an increased tendency of the frequency of 1q21+.

1q21+ often coexisted with high-risk chromosomal abnormalities at high frequency, and was found to be an independent prognostic factor for inferior outcomes in the era of highly-active Proteasome Inhibitor (PI) and Immunomodulatory Drug (IMiD) [3]. For better evaluation of prognosis of MM, a simplified and powerful prognostic staging system called revised International Staging System (R-ISS) has been constructed, it consists high-risk chromosome abnormalities [17p-, t(14;16), t(14,20)], and clinical parameters [4]. Patients with R-ISS third-stage disease and 1q21+ have an inferior outcome comprising an ultra-high-risk group [5]. Considering the significance of 1q21+, they were incorporated into R-ISS to form a revised staging system (R2-ISS), allowing a better stratification of patients with intermediate-risk NDMM [6].

Though the important role of 1q21+ was identified, the mechanisms underlying the development of 1q21+ are controversial, whether it occurs as a byproduct of high-risk biology and genomic instability or increased expression of genes on 1q, such as CKS1B, MCL1, ADAR1, and PDZK1 [3], which is also important for effective therapy. Therefore, the purpose of this study was to explore the potential genetic drivers of 1q21+ and to identify potential therapeutic drugs targeted against these abnormalities.

Material and method

Microarray data and DEGs identification

The gene expression profile of GSE2658 datasets was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). A total of 248 pre-treatment samples of bone marrow from MM patients, which were collected from the datasets after deleting the missing data of 1q21 status, including 114 samples with 1q21+ (amplification or gain) and 134 standard samples. 1q21 statuses were detected by FISH in the dataset.

Differentially expressed genes (DEGs) between 1q21+ and standard MM were analyzed using the online analytical tool SangerBox [7] (http://sangerbox.com/), which is based on R package ‘limma(version 3.40.6)’. The cutoff criteria were adjusted to P < 0.05 and |log2FoldChange| > 1 (|log2FC| > 1). All genes were visualized by volcano plot, and the DEGs were shown by heatmap.

GO and KEGG pathway analysis

Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway analysis were performed for DEGs (P < 0.05) using SangerBox, which is based on R package ‘clusterProfiler(version 3.14.4)’. GO analysis includes three terms: biological process (BP), molecular functions (MF), and cellular components (CC). P-values < 0.05 and counts ≥2 were considered significantly enriched.

WGCNA of candidate genes

WGCNA (Weighted correlation network analysis) was also performed on Sangerbox. First, the samples were clustered, and outliers were removed to ensure the accuracy of the analysis. The soft threshold parameter that satisfied the scale-free co-expression network relationship was set to 5. Second, to reduce interference from spurious correlations, the adjacency matrix was replaced by a topological overlap matrix (TOM). The corresponding dissimilarity was calculated using the dynamics cut tree algorithms to identify hierarchical clustering genes, and the minimum module size was set to 30 to obtain large modules. Third, different module eigenvectors (MEs) were achieved to merge modules with high similarity at a height cut of 0.75. The relevance between modules and clinical traits is shown with a heatmap to determine the modules most closely associated with 1q21+. The |correlation coefficient|>0.5 and P-value < 0.05 were cutoffs for module screening. Finally, the module membership (mM) and gene significance (GS) of individual genes were measured, and the genes in the key modules with |mM| > 0.6 and |GS| > 0.2 were selected for further analysis.

Hub gene identification and evaluation

The intersection of candidate genes from WGCNA and DEGs was constructed to identify the hub genes, and a Venn diagram was used to display the overlapping genes. The expression data of genes were extracted to evaluate the ability of hub genes to recognize 1q21+ MM samples, and receiver operator characteristic (ROC) curves were drawn using the R package ‘pROC (version 1.17.0.1)’.

Survival analysis of hub gene

The survival curve was calculated using Kaplan–Meier analysis and the log-rank test. The difference was considered statistically significant when P < 0.05. High and low expression group of hub genes were divided according to the optimal threshold determined by the ROC curve. The influences of hub genes on survival were conducted by R packages ‘survival (version 3.2.11)’.

The construction of drug-gene interaction

The QuartataWeb [8] (http://quartata.csb.pitt.edu/) supports searching for information on drug-gene interaction. The identified hub genes were uploaded into the database to screen the existing top 20 drugs or compounds. The flow chart of this study is shown in Figure 1.

Figure 1. Flowchart of the study. GSE, GEO Series; DEGs, differentially expressed genes; WGCNA, weighted gene co-expression network analysis; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ROC curve, receiver operator characteristic curve.

Results

DEGs screening

According to the status of 1q21, there were 55 DEGs (supplementary Table 1) between 1q21+ and standard MM samples, including 28 up-regulated genes and 27 down-regulated genes. All DEGs were described by the volcano plot (Figure 2(a)), and the heatmap shows the expression of the DEGs (Figure 2(b)).

Figure 2. Screening of DEGs. a Volcano plot with cutoff criteria set to adjusted P-value < 0.05 and |log2FC| > 1. Red triangles indicate Up-regulated genes, and blue triangles indicate Down-regulated genes. b Heatmap of the 55 DEGs according to the values of |log2FC|. DEGs, differentially expressed genes; FC, fold change.

Function and pathway analysis

GO and KEGG pathway analyses were performed to investigate the functions of DEGs. The top 10 KEGG and GO pathways were selected for display by bubble diagram according to the rank of enrichment count. The KEGG analysis of up-regulated DEGs showed that the pathways associated with transcriptional misregulation in cancer, pathways in cancer and p53 signaling pathway were significantly enriched (Figure 3(a)). The GO analysis of up-regulated DEGs was mainly enriched in the regulation of molecular function, endomembrane system, and system development (Figure 3(b)). While down-regulated DEGs were significantly associated with Kaposi sarcoma-associated herpesvirus infection, human cytomegalovirus infection, and PI3K-Akt signaling pathway in KEGG analysis (Figure 3(c)). For GO analysis, the results indicated that down-regulated DEGs were mostly located in negative regulation of cellular process, cell surface receptor signaling pathway, and immune system process (Figure 3(d)).

Figure 3. KEGG pathway and GO analysis of DEGs. The bubble diagrams show the top 10 functional and results with significant differences. a KEGG of up-regulated DEGs. b GO of up-regulated DEGs. c KEGG of down-regulated DEGs. d GO of down-regulated DEGs. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes.

Candidate genes in WGCNA

The overall clustering and clinical characteristic heatmap of the training set are shown in Figure 4(a). As shown in Figure 4(b), when the soft threshold was set to 3 (R² = 0.85), the interaction between genes conformed to the scale-free distribution to the maximum extent, and the mean value of the adjacency function gradually approached 0. A total of 18 modules were identified by average linkage hierarchical clustering based on the expression values of 10,275 genes, shown in Figure 4(c). The heatmap in Figure 4(d) shows the correlation between different modules. The heat map of pertinence between modules and clinical traits is shown in Figure 4(e). We found that the black module (|cor| = 0.7, P-value = 2.3e-37) was significantly positively associated with 1q21+. The black module contained 475 genes. GS and mM were calculated for each gene to draw scatterplots (Figure 4(f)). According to the |mM| > 0.6 and |GS| > 0.2 conditions, a total of 126 genes were obtained.

Figure 4. The process of WGCNA. a Sample clustering with the trait heatmap. b the left panel shows the scale-free fit index (y-axis) as a function of the soft thresholding power (x-axis); the right panel displays the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis); c Clustering dendrograms for the 10275 genes with dissimilarity based on the topological overlap together with the assigned module colors; eighteen co-expression modules were constructed with various colors. d Module-trait relationships; each row corresponds to a module eigengene, each column corresponds to a trait, and each cell consists of the corresponding correlation and P-value, which are color-coded by correlated according to the color legend. e The heatmap shows the eigengene adjacency. f Scatterplot of 475 genes in the black module; the correlation and P-value are in the upper left corner. WGCNA, weighted gene co-expression network analysis.

Hub gene identification and evaluation

Finally, the intersections of DEGs and WGCNA are three genes (NVL, IL6R, DUSP23) (supplementary Figure 1), which are considered as the candidate genes. Subsequently, the ROC curves of these genes were plotted. As shown in Figure 5(a), the AUCs (area under the curve) of candidate genes were all greater than 0.7, with good performance. Therefore, all are regarded as the hub genes.

Figure 5. a ROC curves of candidate genes. b Survival curves of hub genes.

Survival analysis of hub gene

It was found that the overall survival rate of multiple myeloma with high expression of DUSP23, IL6R, and NVL were significantly lower (Figure 5(b)). Therefore, the three target genes were not only highly correlated with 1q21+, but also consistently had adverse effects on prognosis.

Drug-gene interaction of hub gene

The three hub genes were uploaded into the database for drug-gene interaction analysis. Only DUSP23 and IL6R were identified as interaction drugs (supplementary Figure 2). DUSP23 was matched to one known interaction-drug, IL6R has two known drugs approved by FDA (Table 1), and the others are predicted drugs.

Table 1. The known interaction-drugs targeted to hub genes.

Gene ID	Drugs ID	Drug name	Drug type	Drug group
DUSP23	DB04214	4-Nitrophenyl Phosphate	Small Molecule-Drug	Experimental
IL6R	DB06273	Tocilizumab	Biotech-Drug	Approved
	DB11767	Sarilumab	Biotech-Drug	Approved; Investigational

Discussion

Detection of 1q21 by FISH has been included in the International Myeloma Working Group (IMWG) consensus for risk stratification in multiple myeloma [9]. Through genomic analysis, patients with amplification of CKS1B (at 1q21) in the context of ISS-III or bi-allelic TP53 inactivation were referred to as Double-Hit and had a dire prognosis [10]. Though 1q21+ has been considered as a prognostic risk factor for MM, there wasn’t clear evidence about consistent dysregulations of genes or pathways among MM patients with 1q21+. Hence, we screened the DEGs between 1q21+ and standard MM samples and verified the ability of hub genes as biomarkers for 1q21+. The biological functions of DEGs and their significantly related signaling pathways were further explored, and potentially interactional drugs for hub genes were identified. In this study, three target genes (NVL, IL6R, DUSP23 located in 1q region) associated with 1q21+ in MM patients were obtained by overlapping with DEGs-up regulation.

Interleukin-6 (IL6) is a potent pleiotropic cytokine that regulates plasma cell growth via the IL-6 receptor (IL6R). IL6/IL6R autocrine activity is implicated in the development and progression of MM and other cancers [11]. Myeloma cells release soluble IL6R about immortalization in vitro [12]. In the further analysis of KEGG enrichment of IL6R, we found that it was mainly enriched in pathways of pathways in cancer, PI3K-Akt signaling pathway, viral protein interaction with cytokine and cytokine receptor, human cytomegalovirus infection, Jak-STAT signaling pathway, hematopoietic cell lineage, HIF-1 signaling pathway, and cytokine-cytokine receptor interaction. A study confirmed that IL6R and ADAR1 collaborated to induce a hyper-activation of the oncogenic STAT3 pathway in myeloma, in association with 1q21 amplification [13]. In addition, soluble IL6R levels were increased, significantly correlated to 1q21 amplification, and associated with lower overall survival [14]. MM patients with strong IL6R expression by immunohistochemistry have poor responses to thalidomide-dexamethasone treatment [11]. The standard of care in MM has evolved without much evidence regarding the relative benefit of novel combinations on patients with 1q21+, including PI/IMiDs, Daratumumab and auto stem cell transplantation (ASCT) [3,15]. By searching the database, Sarilumab and Tocilizumab were potential therapeutic drugs of IL6R that have been used in rheumatoid arthritis [16]. The expression of IL6R alpha is restricted to neoplastic plasma cells, not normal plasma cell [17]. These distinct expression patterns suggested that the use of antibody against IL6R may be a promising approach for treatment of 1q21+ MM.

DUSP23 is a low molecular weight DUSP (Dual-specificity protein phosphatase) that participates in interferon and interleukin signaling pathways by dephosphorylating signal transducer and activator of transcription (STAT) proteins [18]. Upregulation of VH1-like (member Z) dual specific protein phosphatase encoded by DUSP23 gene is associated with progression of various human solid tumors, in particular breast cancer [19]. NVL encodes nuclear protein designated NVLp (nuclear VCP-like protein), which is involved in an ATP-dependent nuclear process [20]. It has been proved to be prognostic outlier genes for metastatic progression and overall survival in prostate cancer [21]. Little research has been done about DUSP23 and NVL in MM. We found that DUSP23 and NVL were related to 1q21+ and affected the prognosis of MM. So the mechanism and signal pathway in MM with 1q21+ needs further exploration.

Despite these meaningful findings, our study still has several limitations need to be acknowledged. Firstly, the small sample size may restrict the statistical power and generalizability of these findings. Secondly, this is an exploratory study, and external validation, especially further experimental and clinical evidences are needed to validate the mechanism and clinical value of the hub genes obtained in this study. Moreover, MM is a heterogeneous disease characterized by diverse symptoms, subtypes, and clinical stages, it is important to explore mechanisms of specific MM subtypes for better understanding and treatment for MM.

Contributions

ZX and JN wrote the manuscript and performed the data analysis; YC revised the manuscript and supervised the project. ZX and JY contributed equally to this work and share first authorship. All authors read and approved the final manuscript.

Ethics statement

Because the data were accessed in the public database and identified information is marked, so no ethics approval is needed.

Acknowledgements

The authors thank Xingyue Jin for writing support.

Disclosure statement

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

Data availability statement

The datasets GSE2658 presented in this study can be found in [GEO database] [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE2658].