Multiple Myeloma: Bioinformatic Analysis for Identification of Key Genes and Pathways.

Multiple myeloma (MM) is a hematological malignancy in which monoclonal plasma cells multiply in the bone marrow and monoclonal immunoglobulins are overproduced in older people. Several molecular and cytogenetic advances allow scientists to identify several genetic and chromosomal abnormalities that cause the disease. The comprehension of the pathophysiology of MM requires an understanding of the characteristics of malignant clones and the changes in the bone marrow microenvironment. This study aims to identify the central genes and to determine the key signaling pathways in MM by in silico approaches. A list of 114 differentially expressed genes (DEGs) is important in the prognosis of MM. The DEGs are collected from scientific publications and databases (https://www.ncbi.nlm.nih.gov/). These data are analyzed by Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) software (https://string-db.org/) through the construction of protein-protein interaction (PPI) networks and enrichment analysis of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, by CytoHubba, AutoAnnotate, Bingo Apps plugins in Cytoscape software (https://cytoscape.org/) and by DAVID database (https://david.ncifcrf.gov/). The analysis of the results shows that there are 7 core genes, including TP53; MYC; CDND1; IL6; UBA52; EZH2, and MDM2. These top genes appear to play a role in the promotion and progression of MM. According to functional enrichment analysis, these genes are mainly involved in the following signaling pathways: Epstein-Barr virus infection, microRNA pathway, PI3K-Akt signaling pathway, and p53 signaling pathway. Several crucial genes, including TP53, MYC, CDND1, IL6, UBA52, EZH2, and MDM2, are significantly correlated with MM, which may exert their role in the onset and evolution of MM.

Introduction

Multiple myeloma (MM) is one of the incurable hematological diseases, characterized by the proliferation of abnormal plasma cells with distinct cytogenetic characteristics in the bone marrow, representing about 10% to 15% of all hematopoietic tumors. The International Agency for Research on Cancer (IARC) estimates a worldwide incidence of 160 000 cases and a worldwide mortality of 106 000 patients in 2018. However, men are much more affected by this disease than women. Moreover, the median age of patients at diagnosis is about 66 to 70 years old, with 37% of patients being younger than 65 years old.

Multiple myeloma is a complex genomic landscape and it is characterized by many types of chromosomal aberrations. Hyperdiploidy and immunoglobulin heavy chain (IGH) translocations are included as early occurrences. They are present in the precursor stages of monoclonal gammopathy of undetermined significance (MGUS) and latent multiple myeloma and are completely clonal in the majority of cases. Hyperdiploidy is considered as the first type of copy number alteration commonly seen in MM. It is usually associated with a standard risk prognosis. It is also characterized by trisomies of 3 or more of the odd-numbered chromosomes, namely, 3, 5, 7, 9, 11, 15, 19, and 21, whereas the most common IGH translocations are t(4; 14), t(11; 14), t(14; 16), t(14; 20), and t(6; 14) and its prognostic impact depends largely on the partner chromosome. Other alterations in copy number, gains, and losses of the whole chromosome or part of the chromosome occur in the later stages of the disease trajectory and are not considered triggering events. The exception is the gain of the first copy of 1q, which appears to be an early event, whereas subsequent additional gains of 1q are later events. Common gains and losses include del 1p, gain 1q, del 13/13q, del 17p, in addition to del 16q and del 12p. Moreover, point mutations occur later in MM and commonly affect the MAPK pathway, NFKB pathway, DNA damage and repair, plasma cell differentiation, MYC activation, regulation of gene expression, and the cell cycle pathway. These mutations can arise in various subclones, and their influence on disease progression varies depending on environmental factors. KRAS, NRAS, and BRAF were among the most commonly mutated genes involved in the MAPK pathway, where hot-spot activating mutations have been discovered and were present in roughly 40% to 50% of patients with newly diagnosed multiple myeloma. Recurrent mutations have also been observed in FAM46C, CYLD, DIS3, IRF4, TRAF3, TP53, RB1, LTB, SP140, and more, and mutations in therapeutic targets, namely, lenalidomide targets (CRBN, IKZF1, and IKZF3), proteasome subunits, and steroid receptor (NR3C1), can induce treatment resistance. The combination of all these chromosomal aberrations and gene mutations leads to a differential gene expression profile in the plasma cells of patients with MM.

Bioinformatics is one of the newest fields of biological research, which should be broadly considered as the use of mathematics, statistics, and informatics to process and analyze biological data.[8,9] In addition, it is an important element in laboratories that generate, analyze, store, and interpret data from molecular genetic tests. The present study aims to predict interactions between differentially expressed genes (DEGs) to identify central genes and determine key signaling pathways in MM to understand its genetic heterogeneity.

Materials and Methods

To accomplish this integrative analysis, a pipeline was used (Figure 1).

Figure 1.

Pipeline chart of all study analysis steps. DEG indicates differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction.

A set of 114 DEGs in myeloma plasma cells involved in the progression of the disease were collected from different scientific publications[11 -47] and databases (https://www.ncbi.nlm.nih.gov/). The following search terms were used: multiple myeloma, gene expression, heterogeneity, mutational profiles, and genetic predisposition.

We used the GenBank database, available at the following link (https://www.ncbi.nlm.nih.gov/gene) and GeneCards (https://www.genecards.org/) to match each gene to its ID and functional annotation (Table 1).

Table 1.

Differentially expressed genes (DEGs) in multiple myeloma selected from various databases and scientific publications.

Gene ID	Gene symbols	Gene names	Locus	Description
2	A2M	Alpha-2-macroglobulin	12p13.31	https://www.ncbi.nlm.nih.gov/gene/2
84517	ACTRT3	Actin-related protein T3	3q26.2	https://www.genecards.org/cgi-bin/carddisp.pl?gene=ACTRT3
103	ADAR	Adenosine deaminase RNA specific	1q21.3	https://www.ncbi.nlm.nih.gov/gene/103
57379	AICDA	Activation induced cytidine deaminase	12p13.31	https://www.ncbi.nlm.nih.gov/gene/57379
242	ALOX12B	Arachidonate 12-lipoxygenase, 12R type	17p13.1	https://www.ncbi.nlm.nih.gov/gene/242
81611	ANP32E	Acidic nuclear phosphoprotein 32 family member E	1q21.2	https://www.genecards.org/cgi-bin/carddisp.pl?gene=ANP32E
328	APEX1	Apurinic/apyrimidinic endodeoxyribonuclease 1	14q11.2	https://www.ncbi.nlm.nih.gov/gene/328
27350	APOBEC3C	Apolipoprotein B mRNA editing enzyme catalytic subunit 3C	22q13.1	https://www.ncbi.nlm.nih.gov/gene/27350
140564	APOBEC3D	Apolipoprotein B mRNA editing enzyme catalytic subunit 3D	22q13.1	https://www.ncbi.nlm.nih.gov/gene/140564
200316	APOBEC3F	Apolipoprotein B mRNA editing enzyme catalytic subunit 3F	22q13.1	https://www.ncbi.nlm.nih.gov/gene/200316
60489	APOBEC3G	Apolipoprotein B mRNA editing enzyme catalytic subunit 3G	22q13.1	https://www.ncbi.nlm.nih.gov/gene/60489
164668	APOBEC3H	Apolipoprotein B mRNA editing enzyme catalytic subunit 3H	22q13.1	https://www.ncbi.nlm.nih.gov/gene/164668
596	BCL2	B-cell lymphoma 2	18q21.33	https://www.ncbi.nlm.nih.gov/gene/596
607	BCL9	B-cell lymphoma 9	1q21.2	https://www.ncbi.nlm.nih.gov/gene/607
29760	BLNK	B-cell linker	10q24.1	https://www.ncbi.nlm.nih.gov/gene/29760
23476	BRD4	Bromodomain containing 4	19p13.12	https://www.ncbi.nlm.nih.gov/gene/23476
716	C1S	Complément C1S	12p13.31	https://www.ncbi.nlm.nih.gov/gene/716
23492	CBX7	Chromobox 7	22q13.1	https://www.ncbi.nlm.nih.gov/gene/23492
100507056	CCAT1	Colon cancer–associated transcript 1	8q24.21	https://www.ncbi.nlm.nih.gov/gene/100507056
90835	CCDC189 orC16orf93	Coiled-coil domain containing 189	16p11.2	https://www.genecards.org/cgi-bin/carddisp.pl?gene=CCDC189
595	CCND1	Cyclin D1	11q13.3	https://www.ncbi.nlm.nih.gov/gene/595
894	CCND2	Cyclin D2	12p13.32	https://www.ncbi.nlm.nih.gov/gene/894
896	CCND3	Cyclin D3	6p21.1	https://www.ncbi.nlm.nih.gov/gene/896
928	CD9	CD9 molecule	12p13.31	https://www.ncbi.nlm.nih.gov/gene/928
948	CD36	CD36 molecule	7q21.11	https://www.ncbi.nlm.nih.gov/gene/948
973	CD79A	CD79a molecule	19q13.2	https://www.ncbi.nlm.nih.gov/gene/973
975	CD81	CD81 molecule	11p15.5	https://www.ncbi.nlm.nih.gov/gene/975
55536	CDCA7L	Cell division cycle associated 7 like	7p15.3	https://www.genecards.org/cgi-bin/carddisp.pl?gene=CDCA7L
153241	CEP120	Centrosomal protein 120	5q23.2	https://www.ncbi.nlm.nih.gov/gene/153241
54480	CHPF2	Chondroitin polymerizing factor 2	7q36.1	https://www.genecards.org/cgi-bin/carddisp.pl?gene=CHPF2
91851	CHRDL1	Chordin-like 1	Xq23	https://www.ncbi.nlm.nih.gov/gene/91851
1163	CKS1B	CDC28 protein kinase regulatory subunit 1B	1q21.3	https://www.ncbi.nlm.nih.gov/gene/1163
1164	CKS2	CDC28 protein kinase regulatory subunit 2	9q22.2	https://www.ncbi.nlm.nih.gov/gene/1164
4094	c-MAF or MAF	MAF bZIP transcription factor	16q23.2	https://www.ncbi.nlm.nih.gov/gene/4094
1380	CR2	Complement C3d receptor 2	1q32.2	https://www.ncbi.nlm.nih.gov/gene/1380
55790	CSGALNACT1	Chondroitin sulfate N-acetylgalactosaminyltransferase 1	8p21.3	https://www.ncbi.nlm.nih.gov/gene/55790
1521	CTSW	Cathepsin W	11q13.1	https://www.ncbi.nlm.nih.gov/gene/1521
1545	CYP1B1	Cytochrome P450 family 1 subfamily B member 1	2p22.2	https://www.ncbi.nlm.nih.gov/gene/1545
1634	DCN	Decorin	12q21.33	https://www.ncbi.nlm.nih.gov/gene/1634
79961	DENND2D	DENN domain containing 2D	1p13.3-p13.2	https://www.genecards.org/cgi-bin/carddisp.pl?gene=DENND2D
22894	DIS3	DIS3 homolog, exosome endoribonuclease	13q21.33	https://www.genecards.org/cgi-bin/carddisp.pl?gene=DIS3
1788	DNMT3A	DNA methyltransferase 3 alpha	2p23.3	https://www.ncbi.nlm.nih.gov/gene/1788
27335	EIF3K	Eukaryotic translation initiation factor 3 subunit K	19q13.2	https://www.ncbi.nlm.nih.gov/gene/27335
22936	ELL2	Elongation factor for RNA polymerase II 2	5q15	https://www.genecards.org/cgi-bin/carddisp.pl?gene=ELL2
2071	ERCC3	ERCC excision repair 3, TFIIH core complex helicase subunit	2q14.3	https://www.ncbi.nlm.nih.gov/gene/2071
2146	EZH2	Enhancer of zeste 2 polycomb repressive complex 2 subunit	7q36.1	https://www.ncbi.nlm.nih.gov/gene/2146
2167	FABP4	Fatty acid binding protein 4	8q21.13	https://www.ncbi.nlm.nih.gov/gene/2167
2200	FBN1	Fibrillin 1	15q21.1	https://www.ncbi.nlm.nih.gov/gene/2200
55294	FBXW7	F-box and WD repeat domain containing 7	4q31.3	https://www.ncbi.nlm.nih.gov/gene/55294
2261	FGFR3	Fibroblast growth factor receptor 3	4p16.3	https://www.ncbi.nlm.nih.gov/gene/2261
3006	H1-2	H1.2 linker histone, cluster member	6p22.2	https://www.ncbi.nlm.nih.gov/gene?term=(hist1h1c[gene])%20AND%20(Homo%20sapiens[orgn])%20AND%20alive[prop]%20NOT%20newentry[gene]&sort=weight
3105	HLA-A	Major histocompatibility complex, class I, A	6p22.1	https://www.ncbi.nlm.nih.gov/gene/3105
3113	HLA-DPA1	Major histocompatibility complex, class II, DP alpha 1	6p21.32	https://www.ncbi.nlm.nih.gov/gene/3113
3213	HOXB3	Homeobox B3	17q21.32	https://www.ncbi.nlm.nih.gov/gene/3213
3479	IGF1	Insulin-like growth factor 1	12q23.2	https://www.ncbi.nlm.nih.gov/gene/3479
3488	IGFBP5	Insulin-like growth factor–binding protein 5	2q35	https://www.genecards.org/cgi-bin/carddisp.pl?gene=IGFBP5
3514	IGKC	Immunoglobulin kappa constant	2p11.2	https://www.genecards.org/cgi-bin/carddisp.pl?gene=IGKC
3569	IL6	Interleukin 6	7p15.3	https://www.ncbi.nlm.nih.gov/gene/3569
3570	IL6R	Interleukin 6 receptor	1q21.3	https://www.ncbi.nlm.nih.gov/gene/3570
3608	ILF2	Interleukin enhancer binding factor 2	1q21.3	https://www.ncbi.nlm.nih.gov/gene/3608
3662	IRF4	Interferon regulatory factor 4	6p25.3	https://www.ncbi.nlm.nih.gov/gene/3662
55818	KDM3A	Lysine demethylase 3A	2p11.2	https://www.ncbi.nlm.nih.gov/gene/55818
10365	KLF2	Kruppel-like factor 2	19p13.11	https://www.ncbi.nlm.nih.gov/gene/10365
3936	LCP1	Lymphocyte cytosolic protein 1	13q14.13	https://www.ncbi.nlm.nih.gov/gene/3936
4023	LPL	Lipoprotein lipase	8p21.3	https://www.ncbi.nlm.nih.gov/gene/4023
151827	LRRC34	Leucine-rich repeat containing 34	3q26.2	https://www.genecards.org/cgi-bin/carddisp.pl?gene=LRRC34
344657	LRRIQ4	Leucine-rich repeats and IQ motif containing 4	3q26.2	https://www.genecards.org/cgi-bin/carddisp.pl?gene=LRRIQ4
389692	MAFA	MAF bZIP transcription factor A	8q24.3	https://www.ncbi.nlm.nih.gov/gene/389692
9935	MAFB	MAF bZIP transcription factor B	20q12	https://www.ncbi.nlm.nih.gov/gene/9935
9500	MAGED1	MAGE family member D1	Xp11.22	https://www.ncbi.nlm.nih.gov/gene/9500
4170	MCL1	Apoptosis regulator, BCL2 family member	1q21.2	https://www.ncbi.nlm.nih.gov/gene/4170
4193	MDM2	Murine double minute 2	12q15	https://www.ncbi.nlm.nih.gov/gene/4193
4582	MUC1	Mucin 1, cell surface associated	1q22	https://www.ncbi.nlm.nih.gov/gene/4582
4609	MYC	MYC proto-oncogene	8q24.21	https://www.ncbi.nlm.nih.gov/gene/4609
55892	MYNN	Myoneurin	3q26.2	https://www.ncbi.nlm.nih.gov/gene/55892
7468	NSD2	Nuclear receptor–binding SET domain protein 2	4p16.3	https://www.ncbi.nlm.nih.gov/gene/7468
5174	PDZK1	PDZ domain containing 1	1q21.1	https://www.ncbi.nlm.nih.gov/gene/5174
10957	PNRC1	Proline-rich nuclear receptor coactivator 1	6q15	https://www.genecards.org/cgi-bin/carddisp.pl?gene=PNRC1
57580	PREX1	Phosphatidylinositol-3,4,5-trisphosphate-dependent Rac exchange factor 1	20q13.13	https://www.ncbi.nlm.nih.gov/gene/57580
78994	PRR14	Proline-rich 14	16p11.2	https://www.ncbi.nlm.nih.gov/gene/78994
339105	PRSS53	Serine protease 53	16p11.2	https://www.genecards.org/cgi-bin/carddisp.pl?gene=PRSS53
5710	PSMD4	Proteasome 26S subunit ubiquitin receptor, non-ATPase 4	1q21.3	https://www.ncbi.nlm.nih.gov/gene/5710
23475	QPRT	Quinolinate phosphoribosyltransferase	16p11.2	https://www.ncbi.nlm.nih.gov/gene/23475
5888	RAD51	RAD51 recombinase	15q15.1	https://www.ncbi.nlm.nih.gov/gene/5888
25780	RASGRP3	RAS guanyl releasing protein 3	2p22.3	https://www.ncbi.nlm.nih.gov/gene/?term=25780
1102	RCBTB2	RCC1 and BTB domain–containing protein 2	13q14.2	https://www.ncbi.nlm.nih.gov/gene/1102
55159	RFWD3	Ring finger and WD repeat domain 3	16q23.1	https://www.genecards.org/cgi-bin/carddisp.pl?gene=RFWD3
9810	RNF40	Ring finger protein 40	16p11.2	https://www.ncbi.nlm.nih.gov/gene/9810
6152	RPL24	Ribosomal protein L24	3q12.3	https://www.ncbi.nlm.nih.gov/gene/6152
6161	RPL32	Ribosomal protein L32	3p25.2	https://www.ncbi.nlm.nih.gov/gene/6161
11224	RPL35	Ribosomal protein L35	9q33.3	https://www.ncbi.nlm.nih.gov/gene/11224
6181	RPLP2	Ribosomal protein lateral stalk subunit P2	11p15.5	https://www.ncbi.nlm.nih.gov/gene/6181
6203	RPS9	Ribosomal protein S9	19q13.42	https://www.ncbi.nlm.nih.gov/gene/6203
6223	RPS19	Ribosomal protein S19	19q13.2	https://www.ncbi.nlm.nih.gov/gene/6223
710	SERPING1	Serpin family G member 1	11q12.1	https://www.ncbi.nlm.nih.gov/gene/710
51548	SIRT6	Sirtuin 6	19p13.3	https://www.ncbi.nlm.nih.gov/gene/51548
6635	SNRPE	Small nuclear ribonucleoprotein polypeptide E	1q32.1	https://www.ncbi.nlm.nih.gov/gene/6635
11262	SP140	SP140 nuclear body protein	2q37.1	https://www.ncbi.nlm.nih.gov/gene/11262
6850	SYK	Spleen-associated tyrosine kinase	9q22.2	https://www.ncbi.nlm.nih.gov/gene/6850
54855	TENT5C or FAM46C	Terminal nucleotidyltransferase 5C	1p12	https://www.genecards.org/cgi-bin/carddisp.pl?gene=TENT5C
7018	TF	Transferrin	3q22.1	https://www.ncbi.nlm.nih.gov/gene/7018
10043	TOM1	Target of myb1 membrane trafficking protein	22qf12.3	https://www.ncbi.nlm.nih.gov/gene/10043
7157	TP53	Tumor protein p53	17p13.1	https://www.ncbi.nlm.nih.gov/gene/7157
57212	TP73-AS1	TP73 ARN antisense 1	1p36.32	https://www.genecards.org/cgi-bin/carddisp.pl?gene=TP73-AS1
7295	TXN	Thioredoxin	9q31.3	https://www.ncbi.nlm.nih.gov/gene/7295
7311	UBA52	Ubiquitin A-52 residue ribosomal protein fusion product 1	19p13.11	https://www.ncbi.nlm.nih.gov/gene/7311
9898	UBAP2L	Ubiquitin-associated protein 2 like	1q21.3	(https://www.genecards.org/cgi-bin/carddisp.pl?gene=UBAP2L).
55585	UBE2Q1	Ubiquitin conjugating enzyme E2 Q1	1q21.3	https://www.ncbi.nlm.nih.gov/gene/55585
29089	UBE2T	Ubiquitin conjugating enzyme E2T	1q32.1	https://www.ncbi.nlm.nih.gov/gene/29089
7398	USP1	Ubiquitin-specific peptidase 1	1p31.3	https://www.ncbi.nlm.nih.gov/gene/7398
7874	USP7	Ubiquitin-specific peptidase 7	16p13.2	https://www.ncbi.nlm.nih.gov/gene/7874
7412	VCAM1	Vascular cell adhesion molecule 1	1p21.2	https://www.ncbi.nlm.nih.gov/gene/7412
1462	VCAN	Versican	5q14.2-q14.3	https://www.ncbi.nlm.nih.gov/gene/1462
10413	YAP1	Yes1-associated transcriptional regulator	11q22.1	https://www.ncbi.nlm.nih.gov/gene/10413

Abbreviation: DEG, differentially expressed genes.

Analysis Software

STRING software analysis

Analysis by STRING version 11.5 (https://string-db.org/) allowed us to reveal key genes and co-expressed genes in MM, to identify enriched biological terms, GO terms, and finally to determine key signaling pathways of core genes. FDR value ⩽ 0.05 was considered as the cut-off criterion, with a confidence score of 0.400 set as the cut-off criterion.

Determination of key genes

The input file, composed of 114 DEGs collected, has been submitted for analysis by STRING. The parameters retained for the analysis were the type of organism (Homo sapiens), type of network (evidence network), and confidence score (<0.400). The analysis was performed without enrichment. The core genes were selected based on their number of interactions with other genes (⩾20 interactions) and their location in the PPI network (at the center).

GO analysis

The GO enrichment analysis includes biological process, cellular component, and molecular function. GO terms have been analyzed from files corresponding to the GO enrichment, then downloaded from the software in tabular separated values (TSV) format, which can be opened in Excel as simple tabular text output. The top GO analysis results were selected with their significant values of FDR.

Identification of co-expressed genes

To determine the co-expressed genes, a file corresponding to the string interaction was downloaded from the software in TSV format, which can be opened in Excel as a simple tabular text output. We allow classifying the co-expressed genes by their scores, which indicates the level of association of the expression data and can also be determined by the presence of a black line between the nodes in the PPI network.

Determination of key signaling pathways

They were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) process in String software and selected with their significant FDR value ⩽ 0.05.

Cytoscape software analysis

To verify the results obtained by String software, Cytoscape version 2.3.8 (https://cytoscape.org/) was used, considering the P value < .05 and a confidence score of 0.400 as cut-off criteria.

Key genes identification

The TSV file of the gene’s interaction was imported into Cytoscape software. The CytoHubba Cytoscape plugin was used to classify the network nodes according to their characteristics. The degree is one of the topological analysis methods provided by CytoHubba, allowing us to observe target genes with higher degrees, which could constitute key genes.

GO enrichment

The Cytoscape BiNGO plugin was used to provide the GO enrichment. The most significant GO terms were selected with their P value < .05 and their number of annotated genes.

Network clusters

Cytoscape software groups PPI networks generated by STRING into clusters, taking the P value < .05 and a confidence score of 0.400 as cut-off criteria. ClusterMaker Apps from AutoAnnotate application was used to perform Markov Clustering (MCL) of the protein network. Markov Clustering makes it possible to visualize the enriched terms in the form of circular plots on the nodes of the network.

DAVID functional annotation bioinformatics microarray analysis

The DAVID database version 6.8 (https://david.ncifcrf.gov/) was used to perform KEGG pathway enrichment, to determine the most significant signaling pathways and compare them with those obtained by STRING. The data were pasted as a list of official gene symbols into the DAVID database.

Results

Network analysis

The network analysis of the set of 114 DEG revealed that 110 were annotated by STRING (Figure 2). A central network was obtained encompassing 99 genes, and a small network at the periphery, grouping together 2 genes, APOBEC3H and APOBEC3C. In addition, 9 genes have been outside the core network: RASGRP3; CEP120; MAGED1; CHRDL1; PNRC1; ALOX12B; PREX1; RCBTB2; DENND2D, and 7 core genes constitute the network engine: TP53; MYC CDND1; IL6; UBA52; EZH2; MDM2. In this functional network, all 101 genes represented various interactions, which may be explained by functional links between them.

Figure 2.

Network of protein-protein interaction. The network view (evidence view) summarizes the set of predicted associations for a group of 110 genes. The nodes of the network are the gene product, and the edges represent the predicted functional associations. The edges are represented by lines of different colors that indicate the type of interaction to predict the associations. Clicking on a node will give detailed information about the protein and clicking on an edge will display a detailed breakdown of the evidence.

Co-expression results

From the first network, we determined the co-expression related to the core target genes of MM (Table 2). Key genes have been identified by their scores, which indicate the level of association of the gene expression data. The nodes represented genes and the black line shown between the nodes in the PPI network indicated co-expression (Figure 2). A co-expression score existing between 2 genes represented it. The score indicated the level of association of expression data during a process. If 2 genes showed similar expression under different conditions, it was likely that they were jointly involved in the same process (one requiring the other).

Table 2.

Co-expressed between core genes.

Node 1	Node 1	Score of co-expression
CCND1	MYC	0.069
CCND1	MUC1	0.071
CCND1	CKS1B	0.068
CCND1	FGFR3	0.065
CCND1	CKS2	0.065
CCND1	YAP1	0.317
CCND1	IGFBP5	0.098
EZH2	MYC	0.091
EZH2	TP53	0.074
EZH2	WHSC1	0.120
EZH2	USP1	0.152
EZH2	RAD51	0.212
IL6	CYP1B1	0.086
IL6	FABP4	0.058
IL6	VCAM1	0.063
MYC	TP53	0.070
MYC	RAD51	0.069
MYC	CDCA7L	0.096
MYC	MUC1	0.065
MYC	APEX1	0.070
MYC	MCL1	0.065
MYC	TXN	0.063
TP53	RFWD3	0.073
TP53	RAD51	0.073
TP53	CKS1B	0.076
TP53	RPS19	0.087
TP53	UBE2T	0.062
UBA52	RPL32	0.531
UBA52	RPL35	0.323
UBA52	RPLP2	0.259
UBA52	RPS9	0.112
UBA52	CKS1B	0.065
UBA52	RPL24	0.211
UBA52	RPS19	0.548

The genes marked in bold signify co-expression between 2 key genes.

Among the identified core genes, some of them represented co-expression with each other, including CCND1 with MYC, EZH2 with MYC, EZH2 with TP53, and MYC with TP53. In addition, they also represented co-expression with other genes in the functional network, which were presented in Table 1.

GO enrichment

GO terms were classified in the STRING software according to the value of the false discovery rate and the number of core genes they contained. Moreover, they were classified according to the number of annotated genes in each GO term. Thus, the most important ones had a high number of key genes involved in the development of MM. The GO analysis allowed the obtaining of a total of 448 GO items, including 392 biological process (BP) entries, 23 molecular function (MF) entries, and 33 cellular component (CC) entries (Figure 3).

Figure 3.

Gene Ontology enrichment. The diagram represents the most significant GO terms according to the number of genes involved in the network, which are indicated in parenthesis in the diagram.

GO indicates Gene Ontology.

The most significant biological process determined with FDR values of 7.55e-09; 1.26e-08; 1.26e-08 were as follows: Cellular macromolecule metabolic process (GO:0044260) (64/114 genes), regulation of gene (70/114 genes), respectively.

The most significant molecular function observed with FDR values of 0.00011; 0.00019; 0.00019; were as follows: RNA binding (GO:0003723) (28/114 genes), nucleic acid binding (GO:0003676) (46/114 genes), and protein binding (GO:0005515) (66/114 genes), respectively.

Potentially important target genes have been expressed in the membrane-bounded organelle (GO:0043227), organelle (GO:0043226), and intracellular organelle (GO:0043229), which have a number of annotated genes of 97/114; 100/114; 96/114 and FDR values of 2.80e-06; 1.45e-05; and 1.45e-05, respectively.

Identification of KEGG pathways

KEGG pathway enrichment analyses were performed to reveal potential signaling pathways in the 114 DEGs (Table 3). They were significantly enriched in Epstein-Barr virus (EBV) infection (hsa05169), MicroRNAs in cancer (hsa05206), PI3K-Akt signaling pathway (hsa04151), and p53 signaling pathway (hsa04115), which were considered as the major pathways involved during the development of MM.

Table 3.

KEGG pathways enrichment.

KEGG pathway	Description	Count in network	Number of central genes	False discovery rate
hsa05169	Epstein-Barr virus infection	14	5	4.49e-09
hsa05206	MicroRNAs in cancer	11	5	7.12e-07
hsa04151	PI3K-Akt signaling pathway	13	5	1.85e-05
hsa04115	p53 signaling pathway	7	3	3.33e-05

Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

The table represents the most important KEGG pathways, selected with their significant FDR value ⩽ 0.05 and their number of genes annotated in the functional network.

Identification of key genes

Cytoscape analysis showed 110 annotated genes. CytoHubba apps from the Cytoscape program allowed users to determine the top genes with the “degree” method, where the target genes have the highest degrees, and which are often key target genes (Figure 4). Each node showed a different color depth depending on its own degree (the darker the color, the higher the degree). The network provided 7 important key target genes involved in the regulation of many other genes in the PPI network, which may have potential therapeutic targets in MM.

Figure 4.

Network of genes interaction degree. Proteins with a higher degree of importance are more likely to be essential.

GO enrichment by BINGO plugin

The analysis by the BINGO apps in Cytoscape allowed identifying a set of GO terms with a high number of annotated genes and a significant P value (Supplementary Table S1). These results supported the results obtained by STRING software, and the most GO common were:—Biological process: GO ID 48518 (positive regulation of biological process); GO ID 10467 (gene expression); GO ID 44260 (cellular macromolecule metabolic process)—Molecular function: GO ID 5515 (protein binding)—cellular component: GO ID 43227 (membrane-bounded organelle); GO ID 43229 (intracellular organelle); GO ID 43226 (organelle), and they were marked in bold in Table S1. All GO terms identified in this analysis could be implicated in the development of MM.

The 110 annotated genes were involved in 13 clusters, of which 12 were linked, forming a large network, and 1 remained outside this large network named Apolipoprotein deaminase single, with 2 genes (APOBEC3H and APOBEC3C) (Figure 5). The subgroups of this large network were named according to their protein annotations, and they were classified according to the number of genes involved: Susceptibility apoptosis cell is the largest subgroup and the most important, containing 51 genes, of which 6 were top genes (TP53, c-MYC, CCND1, EzH2, IL6, and MDM2);

Figure 5.

Clustered protein association network. The cluster network provides 13 groups (clusters) of all 114 myeloma gene sets. Each node presents a gene or gene product, the color of the node indicates the 3-dimensional structure of the protein. Each gene cluster indicates a biological function.

Protein us4 ribosomal with 10 genes of which one (UBA52) is the top gene;

Reattaching heterochromatin mitotic with 8 genes;

Chondroitin proteoglycan microfibrils with 5 genes;

Proteinase c1s c1r with 4 genes;

na chloride scarb1 with 4 genes;

Somatomedin insulin higher with 3 genes;

Single deaminase independent with 3 genes;

Removes h2a h2afz with 3 genes;

Exosome 3’ untranslated with 3 genes;

dm mhc class with 2 genes;

Cytolytic release cytochrome with 2 genes.

DAVID database results

To validate the results obtained by STRING and Cytoscape, an analysis by the DAVID software was carried out. The results obtained showed 10 GO terms representing 10 signaling pathways classified according to their P value and the number of genes involved in them (Table 4). The most significant signaling pathways were EBV infection, PI3K-Akt signaling pathway, MicroRNAs in cancer, and the p53 signaling pathway. These results consolidated those obtained by STRING software.

Table 4.

KEGG pathways enrichment by DAVID.

Category	Term	Genes	Count	%	P value	Benjamini
KEGG_PATHWAY	Epstein-Barr virus infection		10	9,1	2,3E-6	3,2E-4
KEGG_PATHWAY	PI3K-Akt signaling pathway		13	11,8	1,0E-4	7,2E-3
KEGG_PATHWAY	MicroRNAs in cancer		11	10,0	4,0E-4	1,5E-2
KEGG_PATHWAY	p53 signaling pathway		6	5,5	4,6E-4	1,5E-2
KEGG_PATHWAY	B cell receptor signaling pathway		6	5,5	5,2E-4	1,5E-2
KEGG_PATHWAY	Bladder cancer		5	4,5	6,3E-4	1,5E-2
KEGG_PATHWAY	Pathways in cancer		12	10,9	1,3E-3	2,4E-2
KEGG_PATHWAY	Small-cell lung cancer		6	5,5	1,4E-3	2,4E-2
KEGG_PATHWAY	FoxO signaling pathway		7	6,4	1,8E-3	2,8E-2
KEGG_PATHWAY	Ribosome		7	6,4	2,0E-3	2,8E-2

Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

Pathways found by DAVID database, selected with their significant P value and their number of genes annotated.

Discussion

In this study, we performed a bioinformatics analysis of a biological data set to identify the key genes and signaling pathways involved during the development of MM. The functional PPI network resulting from the STRING database and Cytoscape software identified that among 114 DEGs introduced, 7 are considered as key genes that represented a high connectivity in the network. Including TP53; MYC, CDND1, IL6, UBA52, EZH2, and MDM2. The CCND1 with MYC, EZH2 with MYC, EZH2 with TP53, and MYC with TP53 are co-expressed, which can be explained by their complementarity and involvement in the same or different processes during the transformation of normal plasma cells into myeloma plasma cells (MM evolution).

Indeed, the study of the co-expression of cell cycle genes has shown the presence of a dynamic balance between the co-expression of sets of genes activating and inhibiting the cell cycle (CDK/cyclin-dependent kinases and CDK1) in MM. The clustering analysis showed that the susceptibility apoptosis cell cluster (the largest subgroup and the most important), containing 6 top genes (TP53, c-MYC, CCND1, EzH2, IL6, and MDM2), could be implicated in the development of MM.[5,16,20,49,50] Nevertheless, the UBA52 gene was found in the protein us4 ribosomal cluster. The ribosomal us4 protein has multiple functions, including mRNA decoding, initiation of small ribosomal subunit aggregation by binding directly to the 16S rRNA and translation repression. It was also involved in the anti-termination activities of transcription.

The ontological analysis revealed the involvement of these 7 hub genes of the functional network in BP related to cellular macromolecular metabolic process, regulation of gene expression, macromolecular metabolic process, and positive regulation of the biological process. These hub genes were the major regulators of these 4 important processes that can be altered during the development of MM.[16,49,52 -54] The CC results suggested that the 7 hub genes were mainly scattered throughout the cell, including the nucleus, mitochondria, plastids, vacuoles, vesicles, ribosomes, and the cytoskeleton. Then, the MF reflected the involvement of hub genes in RNA binding, nucleic acid binding, and protein binding. The alteration in the expression of one or more of these key genes may lead to deregulation of these processes, which can be essential for the cell life cycle and their deregulation may play a role in the development of MM.[16,20,49,50]

was significantly down-regulated in MM, identified as the first hub gene in the PPI network, and it represented the highest degree of connectivity (43 interactions), including the most important interactions with the other 6 major genes identified. It is well known that TP53 is a critical tumor suppressor gene that encodes a tumor suppressor protein with domains for transcriptional activation, DNA binding, and oligomerization. When dysregulated, TP53 plays a key and multifaceted role in cancer development and cancer therapy, including MM. It is implicated in a variety of biological functions, with canonical roles including cell cycle arrest, DNA repair, senescence, and apoptosis, as well as noncanonical roles, such as regulation of cell metabolism, and autophagy. This gene is activated by various stress stimuli and is maintained at a low level by a variety of regulators in normal cells. Its dysregulation was presented in 3 forms in newly diagnosed patients with MM: monoallelic deletion of a part of the chromosome 17p (~8%), monoallelic mutations (~6%), and biallelic inactivation (~4%). Del 17p was still considered a high-risk characteristic in MM and was a part of the current illness staging standards.

The gene was significantly up-regulated in MM. In the PPI network, it was represented in the center with a high degree of connectivity, which involves 36 interactions, including TP53, MYC CDND1, IL6, UBA52, EZH2, and MDM2. It is one of the regulatory and proto-oncogenes that encode transcription factors. It is involved in the regulation of many biological functions, including functions that affect cell growth and proliferation-replication, transcription and translation, cell metabolism, and apoptosis. It has also been recognized as 1 of 4 genes, including Oct4, KLF4, and Sox2, which could together reprogram fibroblasts into a pluripotent stem cell state. This gene is considered a key regulator of MM. Its dysfunction was one of the main features of disease progression, being a trigger for MGUS to MM transition. These alterations include translocations, which have been observed in 50.1% of patients with newly diagnosed MM. These translocations involved immunoglobulin (IG) loci (IGH, IGL, IGK) and certain non-IG partners, in particular FAM46C, FOXO3, and BMP6, where the hypermutations associated with APOBEC deregulation are located. These mutational hotspots are often found at MYC breakpoints, indicating their roles in the generation of the MYC translocation in MM. In addition, MYC gains and duplications at the 8q24.21 locus were present in almost 15% of newly diagnosed patients with MM, which was related to shorter survival in univariate analysis.

, a gene encoding Cyclin D1, was also a gene that was upregulated in MM. It has also been identified as a core gene based on its importance in the network. It interacted with 30 genes in the PPI network. Cyclin D1 was a key cell cycle regulator and was involved in the pathogenesis of several cancers. Its overexpression and translocation are frequent events in MM, suggesting that it might be the cause of the initiation and development of this malignancy. During mitosis, the transition of the cell cycle from G1 to S phase is controlled by the cyclin-dependent kinases (CDKs), CDK4 and CDK6, which form protein complexes with cyclin D1. Cyclin D1 catalyzes the phosphorylation of the tumor suppressor protein retinoblastoma (RB) while releasing the transcription factor E2F and triggering the downstream gene transcription required for cell cycle progression. As a result, cyclin D1 inhibition was conducted to cell cycle arrest, while overexpression of the protein accelerates the G1 phase transition. Due to chromosomal duplication, translocations, and alteration of normal intercellular trafficking and proteolysis, abnormal expression of cyclin D1 has been reported in a variety of human tumors. In particular, 15% to 20% of MM samples carry a t(11;14) chromosomal translocation, resulting in abnormal transcriptional activation of CCND1. This most common translocation in MM, t(11; 14) (q13; q32), involves an abnormal fusion of the IGH locus with cyclin D1. Overexpression of cyclin D1 was also observed in 25% to 50% of MM samples, indicating that the deregulation of cyclin D1, which was a critical regulator of the G1/S transition and therefore of the cell cycle, may be a key event in MM development.

In the PPI network, Interleukin-6) represented the connectivity of 25 interactions. This gene was significantly upregulated in MM. IL6 is a pleiotropic cytokine with extensive functions in inflammation and immunity. It has been extensively studied for its role in normal antibody-producing plasma cells. It promotes the normal B-cell differentiation into cells producing antibodies without triggering the proliferation of B cells. It has been shown that IL6 is a key factor in the proliferation of MM, promoting the growth, proliferation, and survival of myeloma cells. Bone marrow stromal cells presented the major source of IL6 in MM cells, and its overexpression was linked to poor prognosis and larger tumor cell mass in MM.

(Enhancer of Zeste 2 Polycomb Repressive Complex 2) represented a high degree of connectivity, with 22 interactions in the PPI network. It is one of the DEGs upregulated in MM. The EZH2 gene was a histone methyltransferase that acts primarily on H3K27 and catalyzes the conversion to a trimethylated marker (H3K27me3). This gene played a critical role in normal B cell development since H3K27me3 expression and levels influence differentiation decisions. Increased EZH2 expression in germinal center B cells resulted in cell cycle checkpoints disappearing and allowing B cells to expand. Subsequently, EZH2 decreased, as a result allowing the cells to differentiate into plasma cells. EZH2 was known to be deregulated in MM. Pawlyn et al discovered for the first time a link between high expression of EZH2 and myeloma survival. This association was robust across different datasets, persisted regardless of the treatment method used, and was independent of other factors known to affect survival in patients with myeloma. This reinforces the importance of EZH2 expression in the pathogenesis of myeloma and suggested that its inhibition may be a potential therapeutic approach for myeloma therapy and should be studied in clinical trials.

(ubiquitin A-52 residue ribosomal protein fusion product 1) was one of the up-regulated genes observed in MM. In the PPI network, it was also recognized as a core gene, representing high connectivity of 22 interactions, the most significant of which are those with the other 6 hub genes. This gene is essential for the selective degradation of proteins by the ubiquitin-proteasome system (UPS). In addition, it has been found to play a central role in cellular functions, including cell cycle control, apoptosis, signaling, and transcriptional regulation. UBA52 has been shown to promote the development of cancer through multiple mechanisms, which may indicate its involvement in the development of MM. UBA52 has been the most novel target gene proposed as a driver of MM until there are enough studies on its involvement in myeloma. Further genomic studies are needed to explain the function of UBA52 in the process of the disease.

(Murine double minute protein) represented important connectivity with about 21 interactions in the network. The MDM2 gene is a pleiotropic protein located on the 12q13-14 chromosome, known to facilitate the ubiquitination of p53 required for its proteasome-mediated turnover. In MM, this gene is highly and constitutively expressed through processes such as chromosomal trisomy or gene amplification and functions both to drive cell proliferation and to enhance cell survival. It activates E2F-1, which promoted the transition from G1 to the S phase, and suppressed the function of wild-type p53 (wtp53), which can enhance cell survival. MDM2 binds constitutively to E2F-1, wtp53, and mtp53 (mutant-type p53) in all myeloma cells, as well as p21 in malignancy cells lacking p53. Then p21 was activated by wtp53 and potentiated the tumor suppressor function of wtp53 by strongly inhibiting cell cycle regulatory proteins such as cyclin E and CDK2. This indicated that MDM2 can improve the cell cycle progression of MM cells by downregulating cell cycle inhibitory proteins (wtp53 and p21) and activating E2F-1. Therefore, overexpression of MDM2 may contribute to the growth and survival of MM cells, indicating that therapeutic strategies targeting MDM2 have potential utility in MM.[16,74]

KEGG pathway enrichment analysis showed that the 7 key genes were significantly enriched in EBV infection (hsa05169), microRNAs in cancer (hsa05206), PI3K-Akt signaling pathway (hsa04151) and p53 signaling pathway (hsa04115). Many studies have been in concordance with our results, which demonstrated the existence of an association between EBV infection and MM development, in which a polymerase chain reaction was performed and revealed the presence of an elevated level of EBV DNA in patients with MM compared with healthy controls.[75,76] MicroRNAs (miRNAs) are small noncoding RNAs aberrantly expressed in solid and hematopoietic malignancies, where they play a central role as post-transcriptional regulators of gene expression. Recently, some studies have demonstrated the efficacy of miRNAs as specific and sensitive biomarkers for the classification, prognosis, and diagnosis of human cancer. Jones et al discovered 3 miRNAs in serum, miR-1308, miR-1246, and miR 720, which have the potential to be used as diagnostic biomarkers in myeloma. They demonstrated that the joint use of miR-1308 and miR-720 provided a powerful diagnostic tool for distinguishing normal healthy controls from patients with MGUS/myeloma. Furthermore, the combination of miR-1246 and miR-1308 can distinguish patients with MGUS from those with myeloma. The PI3K/Akt pathway has attracted considerable attention as a promising therapeutic target in MM. A study conducted on the PI3K/Akt/mTOR pathway (mammalian target of rapamycin) by Vijay Ramakrishnan and Shaji Kumar has shown that it played a critical role in the biology of diseases. This pathway was aberrantly activated in a large proportion of patients with MM by multiple mechanisms. It may also play a role in resistance to several existing therapies, making it a central pathway in the pathophysiology of MM. The study by Vikova et al showed after whole-exome sequencing on a large cohort of 30 human MM cell lines, representative of a large molecular heterogeneity of MM, and 8 control samples, that many canonical pathways known to be implicated in the proliferation and survival of MM cells have been mutated, including PI3K-AKT. In addition, 76% of myeloma cell lines had mutations in the p53 cell cycle pathway genes.

Conclusion

STRING, Cytoscape, and DAVID Software analyzed a set of 114 MM DEGs. TP53, MYC, CCND1, IL6, UBA52, EZH2, and MDM2 were identified as hub genes, which may play a key role in the pathogenesis of MM. They were implicated during different stages of disease progression. The 7 hub genes were identified as being significantly enriched in various pathways, particularly in EBV infection, microRNAs in cancer, the PI3K-Akt signaling pathway, and the p53 signaling pathway, which can be critical for the progression of myeloma. This predictive study can be a relevant tool for a better understanding of the evolution of MM in its complex microenvironment. Understanding the pathogenesis of MM could reveal new biomarkers of myeloma cells and promising potential therapeutic targets.

Click here for additional data file.

Supplemental material, sj-docx-1-bbi-10.1177_11779322221115545 for Multiple Myeloma: Bioinformatic Analysis for Identification of Key Genes and Pathways by Chaimaa Saadoune, Badreddine Nouadi, Hasna Hamdaoui, Fatima Chegdani and Faiza Bennis in Bioinformatics and Biology Insights