Bioinformatics analysis deciphering the transcriptomic signatures associated with signalling pathways and prognosis in the myelodysplastic syndromes

ABSTRACT

Background: Several studies scatteredly identified the myelodysplastic syndromes’ transcriptomic profiles (MDS). However, the exploration of transcriptional signatures, key signalling pathways, and their association with prognosis and diagnosis in the integrated multiple datasets remains lacking. Methods: We integrated the GSE4619, GSE19429, GSE30195, and GSE58831 microarray datasets of CD34 + cells for identifying the differentially expressed genes (DEGs) in the MDS. The series of bioinformatics methods are applied to identify the key hub genes, gene clusters, prognostic hub genes, and genes associated with diagnostic efficacy. Finally, we validated the expression differences of hub genes in the GSE114922 dataset. Results: We explored the DEGs related to gene ontology enrichment and KEGG pathways. We identified significant hub genes, including 168 upregulated hub genes (such as STAT1, IFIH1, EPRS, GRB2, RAC2, MAPK14, CASP1, and SPI1) and 52 downregulated hub genes (such as CREBBP, HIF1A, PIK3CA, EZH2, PIK3R1, MDM2, IRF4, CXCR4, PCNA, and CD19) in the MDS. In addition, we identified six significant molecular complex detection (MCODE)-derived upregulated gene clusters and one downregulated gene cluster, respectively. Moreover, we found that the higher expression level of MX2, GBP2, PXN, IFI44, FDXR, PLCB2, ASS1, ERCC4, PML, and RRAGD and the lower expression level of CD19, PAX5, TCF3, LEF1, NUSAP1, and TIMELESS hub genes are significantly correlated with shorter survival times of MDS patients. Furthermore, the area value under the ROC curve (AUC) of PXN, FDXR, PLCB2, PML, CD19, PAX5, and LEF1 prognostic genes are more than 0.80, indicating that these genes could be effectively used for the diagnostic efficacy of MDS patients. Conclusions: Identifying key hub genes and their association with the prognosis and diagnostic efficacy may provide substantial clues for the treatment and diagnosis of MDS patients.

KEYWORDS:

• Myelodysplastic syndrome
• Bioinformatics
• hub genes
• Survival times
• Diagnostic efficacy

1. Introduction

The most common malignant blood disorder is myelodysplastic syndromes (MDS) characterized by dysfunctional haematopoiesis [1]. The major molecular cause of MDS is genetic instability, which in turn transformed into acute myeloid leukaemia (AML) in a larger number of cases [1]. MDS is frequently occurring in older males and in other individuals who are earlier exposure to cytotoxic therapy [2]. Some previous studies demonstrated the genetic abnormalities and molecular pathogenesis of MDS patients [3–6]. It has been reported that the pathogenesis, initiation of disease, disease progression, and heterogeneity of MDS are crucially associated with the genetic abnormalities [3]. The deregulated transcriptomes have a greater contribution to the development of MDS [4]. The primary hallmark of MDS development is genes or chromosomal changes, copy-number alterations, and somatic mutations [6]. For example, Ying le demonstrated that the dysfunctional genes and signalling pathways are altered during the progression of MDS [5]. Aberrantly expressed transcriptomes are substantially associated with the clinical outcome of MDS patients [4]. Autophagy-associated genes predicted the poor prognosis of MDS patients [7]. Also, deregulated metabolic genes are correlated with poor prognosis of MDS patients [8]. These MDS-associated studies provided a clue that the deregulated and genetically altered transcriptomes are associated with the pathogenesis, development, progression, and clinical outcomes of the MDS patients.

Herein, a meta-analysis integrated the four microarray gene expression datasets of MDS to identify the differentially expressed genes (DEGs). Then, we explored the functional enrichment analysis, included biological processes, molecular functions, cellular components, and KEGG pathways that are significantly associated with the upregulated and downregulated genes, respectively. In addition, we constructed the protein–protein interaction (PPI) network of DEGs and identified the key hub genes based on the degree of interaction of DEGs. We investigated the expression validation of top hub genes in the independent datasets of MDS patients. Moreover, we explored the association of hub genes with the survival prognosis of MDS patients. Finally, we investigated the diagnostic efficacy of prognostic hub genes in MDS patients.

2. Methods and materials

2.1 Datasets

We searched the NCBI gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) using the keywords ‘myelodysplastic syndromes ‘, ‘bone marrow’, ‘CD34 + cells’ ‘ CD34 + cells in MDS’, ‘bone marrow in MDS’ and ‘transcriptomes in MDS’, and identified a five gene expression datasets of CD34 + cells that are associated with MDS of human. We downloaded and integrated the GSE4619 (n = 66) [9], GSE19429 (n = 200) [10], GSE30195 (n = 19) [11], and GSE58831 (n = 176) [4] to identify the DEGs in the MDS. The platform of these datasets is Affymetrix human genome U133 plus 2.0 array. The integrated dataset (combining four datasets) has 461 samples, including 412 MDS samples and 49 healthy controls. Then, we used the GSE114922 (n = 90) [12] for validating the expression differences of key genes between the MDS and control samples. The platform of the GSE114922 is GPL20301 (Illumina HiSeq 4000 (Homo sapiens)). All data were transformed into log2-base normalization. For analyzing the survival differences between the two groups of patients, we downloaded the clinical data of GSE58831[4].

2.2 Identification of differentially expressed genes (DEGs)

We used the web tool NetworkAnalyst software [13] to identify the DEGs between MDS and normal samples. The multiple probes of a single gene were averaged and normalized by log2 normalization. We reduced the batch effects of multiple datasets by using the ComBat method [14]. A total of 19,547 common genes were found after integrating the four datasets. We employed the R package ‘limma’ to identify the significant DEGs between MDS and normal samples [15]. We identified the DEGs with a threshold absolute value of combined effect size (ES) >0.585 and adjusted P-value≤0.05.

2.3 Gene ontology and gene-set enrichment analysis

We performed the DEGs’ functional and gene-set enrichment analysis by using the GSEA [16]. We used the GSEA tool for identifying the gene ontology terms and KEGG pathways for the significant upregulated and downregulated DEGs in the MDS. The biological processes, cellular components, molecular functions, and KEGG [17] pathways that are significantly associated with the upregulated DEGs and the downregulated DEGs were identified, respectively. The FDR-value < 0.05, calculated by the Benjamini – Hochberg method [18], was considered significant when selecting the functionally enriched terms.

2.4 Identifying hub genes and gene clusters from the protein-protein interaction (PPI) network of DEGs

We constructed the PPI network using the STRING database to understand better the regulatory relationship among the identified DEGs [19]. To determine the rank of hub genes, we used Cytoscape plug-in tool cytoHubba [20]. Hub genes were identified based on the degree of interactions with their neighbour genes. To construct the PPI network of DEGs, we selected the minimum required interaction score of 0.40. Hub genes were defined as a gene that was connected to a minimum of 20 other DEGs in the original PPI. We visualize the PPI networks by utilizing the Cytoscape 3.6.1 software A Cytoscape plug-in molecular complex detection (MCODE) was used to detect the significant modules from the original PPI network [21]. We identified the significant modules based on the MCODE score and node number. The threshold of the MCODE analysis was MCODE score > 5.0, Node Score Cut-off: 0.2, Haircut: true, K-Core: 2, and maximum depth from Seed: 100.

2.5 Validation of Hub gene expression levels in the independent dataset

The expression levels differences of top hub genes were further validated using the GSE114922 (n = 90) [12] of MDS patients. The GSE114922 is RNA-sequencing data that was performed on bone marrow CD34 + hematopoeitic stem and progenitor cells from patients with MDS (n = 82) and healthy controls (n = 8). We used the R package ‘limma’ to identify the significant DEGs between MDS and healthy samples [15]. The expression differences of hub genes with P < 0.05 were considered statistically significant between the two groups.

2.6 Survival analysis of Hub genes

We utilized the R package ‘survival’ to investigate the survival prognosis of MDS patients. We compared the overall survival (OS) of MDS patients classified based on the expression levels of all hub genes (High expression group > mean > Low expression group). Kaplan-Meier survival curves were used to show the survival differences, and the log-rank test was utilized to evaluate the significance of survival differences. The prognostic roles of screened hub genes in the MDS were analyzed using the expression and clinical data of GSE58831[4]. Cox P< 0.05 was taken as the significant between two groups of MDS patients.

2.7 Diagnostic efficacy evaluation for the prognostic hub genes in the MDS

To evaluate the diagnostic values of the prognostic key upregulated and downregulated hub genes, the receiver operating characteristic (ROC) curve was plotted and the area under the ROC curve (AUC) was calculated using the ‘pROC’ R package [22], indicating the capability to distinguish MDS and healthy control. The diagnostic roles of screened hub genes in the MDS were analyzed using the expression data of GSE58831[4]. The more excellent AUC value of individual genes indicated the differences between MDS and normal samples, and the key gene with AUC > 0.5 in the MDS datasets was defined as a diagnostic efficiency of the gene [23].

3. Results

3.1 Identification of DEGs in MDS with the integration of four gene expression data

In the integrated datasets, we explored 2404 differentially expressed genes (DEGs), which included 1619 upregulated (Supplementary Table S1) and 785 downregulated (Supplementary Table S2) genes in the MDS datasets when compared with healthy control. The top (the highest combined effect size (ES)) significantly upregulated 25 genes (such as DENND2D, SMIM20, KCNE3, TGDS, ALDH3B1, SNX3, PLEKHA8P1, CDC16, E2F6, RPS27L, COA6, LYL1, BNIP3, and CHCHD10) are displaying in Table 1. In addition, the top 25 downregulated genes (such as OR7A5, GPR176, IRF4, ARPP21, IGHV5-78, LRIG1, AKAP12, DUSP26, MME, RAG1, SH2D4B, and LDLRAD4) are displayed in Table 2.

Table 1. The top upregulated DEGs in MDS with relative to the healthy control.

Entrez ID	Symbol	Name	Combined ES	Adjusted P-Value
79961	DENND2D	DENN domain containing 2D	1.47	0
389203	SMIM20	small integral membrane protein 20	1.42	0
10008	KCNE3	potassium voltage-gated channel subfamily E regulatory subunit 3	1.41	0
23483	TGDS	TDP-glucose 4,6-dehydratase	1.36	3.1E-07
221	ALDH3B1	aldehyde dehydrogenase 3 family member B1	1.34	0
8724	SNX3	sorting nexin 3	1.32	0
51054	PLEKHA8P1	pleckstrin homology domain containing A8 pseudogene 1	1.31	9.44E-14
8881	CDC16	cell division cycle 16	1.31	9.44E-14
1876	E2F6	E2F transcription factor 6	1.30	9.44E-14
51065	RPS27L	ribosomal protein S27 like	1.29	1.63E-07
388753	COA6	cytochrome c oxidase assembly factor 6	1.28	3.34E-13
4066	LYL1	LYL1 basic helix-loop-helix family member	1.28	3.34E-13
664	BNIP3	BCL2 interacting protein 3	1.28	1.4E-12
400916	CHCHD10	coiled-coil-helix-coiled-coil-helix domain containing 10	1.27	3.34E-13
197259	MLKL	mixed lineage kinase domain like pseudokinase	1.27	3.95E-13
84984	CEP19	centrosomal protein 19	1.26	1.22E-07
10538	BATF	basic leucine zipper ATF-like transcription factor	1.26	6.74E-13
10007	GNPDA1	glucosamine-6-phosphate deaminase 1	1.24	1.58E-12
5321	PLA2G4A	phospholipase A2 group IVA	1.24	5.67E-10
7187	TRAF3	TNF receptor associated factor 3	1.24	1.7E-12
54858	PGPEP1	pyroglutamyl-peptidase I	1.23	2.1E-12
51001	MTERF3	mitochondrial transcription termination factor 3	1.22	2.33E-05
4013	VWA5A	von Willebrand factor A domain containing 5A	1.22	2.81E-12
9937	DCLRE1A	DNA cross-link repair 1A [Source:HGNC Symbol;Acc:HGNC:17660]	1.22	3.86E-12
9342	SNAP29	synaptosome associated protein 29	1.22	3.35E-12

Table 2. The top downregulated DEGs in MDS with relative to the healthy control.

Entrez ID	Name	Name	Combined ES	Adjusted P-Value
26659	OR7A5	olfactory receptor family 7 subfamily A member 5	−2.78	7.39E-12
11245	GPR176	G protein-coupled receptor 176	−2.36	0
3662	IRF4	interferon regulatory factor 4	−2.34	0
10777	ARPP21	cAMP regulated phosphoprotein 21	−2.30	0
28387	IGHV5-78	immunoglobulin heavy variable 5–78 (pseudogene)	−2.08	0
26018	LRIG1	leucine rich repeats and immunoglobulin like domains 1	−2.03	0
9590	AKAP12	A-kinase anchoring protein 12	−2.01	0
78986	DUSP26	dual specificity phosphatase 26	−2.00	0
4311	MME	membrane metalloendopeptidase	−1.94	0
5896	RAG1	recombination activating 1	−1.92	0
387694	SH2D4B	SH2 domain containing 4B	−1.91	0
753	LDLRAD4	low density lipoprotein receptor class A domain containing 4	−1.90	0
5897	RAG2	recombination activating 2	−1.87	0
4092	SMAD7	SMAD family member 7	−1.86	0
974	CD79B	CD79b molecule	−1.86	2.97E-07
4852	NPY	neuropeptide Y	−1.85	0
8974	P4HA2	prolyl 4-hydroxylase subunit alpha 2	−1.84	8.19E-11
57824	HMHB1	histocompatibility minor HB-1	−1.83	1.48E-09
1.01E+08	LINC01013	long intergenic non-protein coding RNA 1013	−1.83	0
29802	VPREB3	V-set pre-B cell surrogate light chain 3	−1.82	0
7441	VPREB1	V-set pre-B cell surrogate light chain 1	−1.79	0
7466	WFS1	wolframin ER transmembrane glycoprotein	−1.77	0
94241	TP53INP1	tumor protein p53 inducible nuclear protein 1	−1.73	0
640	BLK	BLK proto-oncogene, Src family tyrosine kinase	−1.71	6.91E-10
114757	CYGB	cytoglobin	−1.71	0.00E+00

3.2 Bone marrow of CD34+ cells-derived transcriptomes are associated with gene ontology and the enrichment of KEGG pathways

We used the GSEA tool to identify the enriched gene ontology terms and KEGG pathways that were associated with identified DEGs (Figures 1 and 2). The identified gene ontology included biological processes, cellular components, and molecular functions (Figure 1A and 1B). The upregulated genes are significantly associated with the enrichment of biological processes (Such as small molecule metabolic process, cellular macromolecule localization, establishment of protein localization, protein containing complex subunit organization, intracellular transport, regulation of intracellular signal transduction, immune effector process, lipid metabolic process, and positive regulation of biosynthetic process) (Supplementary Table S3), cellular components (mitochondrion, envelope, vacuole, mitochondrial envelope, endosome, organelle subcompartment, Golgi apparatus, mitochondrial matrix, Golgi apparatus subcompartment, endoplasmic reticulum) (Supplementary Table S4), and molecular functions (identical protein binding, enzyme binding, transcription regulator activity, DNA binding transcription factor activity, transition metal ion binding, protein dimerization activity, and oxidoreductase activity) (Supplementary Table S5). Besides, downregulated genes are significantly associated with the enrichment of biological processes (Supplementary Table S6), cellular components (Supplementary Table S7), and molecular functions (Supplementary Table S8). Figure 1A and Figure 1B illustrated the top ten upregulated and downregulated GO, respectively.

Figure 1. Significantly enriched GOs that are associated with CD34 + cells-derived DEGs in the MDS. A. Significantly enriched GOs that are associated with CD34 + cells-derived upregulated genes in the MDS. B. Significantly enriched GOs that are associated with CD34 + cells-derived Downregulated genes in the MDS.

Figure 2. Significantly enriched KEGG pathways that associated with CD34 + cells-derived DEGs in the MDS. A. The top 25 enriched KEGG pathways are considerably associated with upregulated genes. B. The KEGG pathways that significantly associated with downregulated genes.

Moreover, we identified the significantly enriched KEGG pathways (Figure 2) for upregulated and downregulated genes, respectively. The upregulated pathways are mainly associated with immune regulation (RIG-I-like receptor signaling pathway, adipocytokine signaling pathway, chemokine signaling pathway, Fc epsilon RI signaling pathway, Toll-like receptor signaling pathway, Fc gamma R-mediated phagocytosis, natural killer cell mediated cytotoxicity, B cell receptor signaling pathway, T cell receptor signaling pathway, epithelial cell signaling in Helicobacter pylori infection, leukocyte transendothelial migration, NOD-like receptor signaling pathway, cytokine-cytokine receptor interaction), metabolisms (lysine degradation, propanoate metabolism, valine, leucine and isoleucine degradation, fatty acid metabolism, beta-Alanine metabolism, butanoate metabolism, Fructose and mannose metabolism, glycerolipid metabolism, glycerophospholipid metabolism, inositol phosphate metabolism, limonene and pinene degradation, amino sugar and nucleotide sugar metabolism, alanine, aspartate and glutamate metabolism, glycosaminoglycan biosynthesis – chondroitin sulfate, sulfur metabolism), cellular signaling and development (lysosome, MAPK signaling pathway, cytosolic DNA-sensing pathway, peroxisome, apoptosis, endocytosis, VEGF signaling pathway, Jak-STAT signaling pathway, homologous recombination, Wnt signaling pathway, phosphatidylinositol signaling system, oxidative phosphorylation, tight junction, nucleotide excision repair, ErbB signaling pathway, basal transcription factors), cancers (acute myeloid leukemia, pathways in cancer, and pancreatic cancer) (Figure 2A and Supplementary Table S9). Besides, we identified significantly enriched KEGG pathways that are associated with downregulated genes in the MDS. Some of these downregulated pathways are B cell receptor signaling pathway, primary immunodeficiency, cell cycle, circadian rhythm – mammal, cytokine-cytokine receptor interaction, p53 signaling pathway, neuroactive ligand–receptor interaction, cell adhesion molecules (CAMs), hematopoietic cell lineage, renin-angiotensin system, Oocyte meiosis, and apoptosis (Figure 2B and Supplementary Table S10).

3.3 STRING-based PPI analysis identifying key hub genes and modules in the MDS

We constructed the PPI of all significant upregulated and downregulated genes to identify the key hub genes and significant gene modules in the MDS. The STRING-based PPI information was inputted into the Cytoscape for further analysis. We revealed that the 1446 upregulated genes are involved with the PPI network (Supplementary Table S11) and 632 downregulated genes are (Supplementary Table S12) associated with the PPI network. Then, we investigated the rank of all upregulated and downregulated genes based on the Cytoscape plug-in CytoHubba tool. We identified the168 upregulated hub genes (minimum degree of interaction is 20 with other neighbours DEGs) (Supplementary Table S11) and 52 downregulated hub genes (minimum degree of interaction is 20 with other neighbours DEGs) (Supplementary Table S12) in the MDS. The PPI networks of upregulated and downregulated hub genes were illustrated in Figure 3A and B, respectively. The top ten upregulated hub genes are STAT1, STAT3, IFIH1, EPRS, GRB2, RAC2, MAPK14, CASP1, SPI1, and TBK1. The identified top ten downregulated hub genes are CREBBP, HIF1A, PIK3CA, EZH2, PIK3R1, MDM2, IRF4, CXCR4, PCNA, and CD19. The top upregulated and downregulated hub genes and their regulatory status, names, the rank of hub genes, degree of interaction, combined ES, and adjusted P values are shown in Table 3.

Figure 3. Identification of the hub genes that are associated with PPI in the MDS. A. The interaction of 168 upregulated hub genes identified from the PPI in MDS. B. The interaction of 52 downregulated hub genes identified from the PPI in MDS. The color intensity of nodes indicates the interconnected degree of hub genes. The more red color indicated the higher the degree of interactions.

Table 3. The top ten upregulated and top ten downregulated hub genes were extracted from the original PPI. The hub genes were identified by using the Cytoscape plug-in CytoHubba tool.

Regulatory status	Rank	Symbol	Name	Degree of interaction	Combined ES	Adjusted P value
Upregulated	1	STAT1	signal transducer and activator of transcription 1	80	0.90872	3.17E-07
	2	STAT3	signal transducer and activator of transcription 3	71	0.93007	1.62E-04
	3	IFIH1	interferon induced with helicase C domain 1	63	0.9983	5.39E-06
	4	EPRS	glutamyl-prolyl-tRNA synthetase 1	57	0.64802	3.63E-04
	5	GRB2	growth factor receptor bound protein 2	56	0.6316	5.26E-04
	6	RAC2	Rac family small GTPase 2	51	0.60996	8.48E-04
	7	MAPK14	mitogen-activated protein kinase 14	50	0.75339	2.71E-05
	7	CASP1	caspase 1	50	0.71999	2.76E-02
	9	SPI1	Spi-1 proto-oncogene	49	0.71557	1.35E-03
	9	TBK1	TANK binding kinase 1	49	0.64902	3.53E-04
Downregulated	1	CREBBP	CREB binding protein	55	−0.6001	1.06E-03
	2	HIF1A	hypoxia inducible factor 1 subunit alpha	50	−0.7376	4.44E-05
	2	PIK3CA	"phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha	50	−0.7585	6.69E-05
	4	EZH2	enhancer of zeste 2 polycomb repressive complex 2 subunit	49	−0.9086	3.05E-07
	5	PIK3R1	phosphoinositide-3-kinase regulatory subunit 1	48	−0.7313	4.85E-05
	6	MDM2	MDM2 proto-oncogene	43	−0.8032	6.97E-06
	7	IRF4	interferon regulatory factor 4	41	−2.335	0.00E+00
	8	CXCR4	C-X-C motif chemokine receptor 4	40	−1.298	9.44E-14
	8	PCNA	proliferating cell nuclear antigen	40	−0.8258	3.65E-06
	10	CD19	CD19 molecule	39	−1.5513	1.86E-08

Furthermore, we extracted the significant gene modules from the original PPI networks by using Cytoscape plug-in MCODE tool. We found the 6 essential gene modules/clusters (MCODE score > 5.0) for the upregulated genes in the MDS (Figure 4) (Supplementary Table S13). The 42 genes are interconnected in cluster 1, including XAF1, OAS1, IFIT1, MX2, IFIH1, GBP2, STAT1, IRF9, IFIT3, and IFITM1 are the top interconnected genes (Figure 4A and Supplementary Table S13). The enrichment analysis revealed the 3 KEGG pathways (RIG-I-like receptor signaling pathway, Toll-like receptor signaling pathway, and Jak-STAT signaling pathway) that are associated with the genes of cluster 1 (Table 4). In cluster 2, some of the mitochondrial ribosomal proteins (MRPS6, MRPL46, MRPS10, MRPS30, and MRPL40) are interconnected in the network (Figure 4B and Supplementary Table S13). The 16 genes, including RXRA, MED8, NCOA1, MED18, FAM120B, MED31, MED20, and MED21 are interconnected in cluster 3 (Figure 4C and Supplementary Table S13). The enrichment analysis identified the 4 KEGG pathways (apoptosis, natural killer cell mediated cytotoxicity, cytokine-cytokine receptor interaction, and p53 signaling pathway) that are significantly correlated with cluster 3 (Table 4). Cluster 4 has 36 interconnected genes, including IKBKB, CASP1, TBK1, TNFRSF1A, TRAF2, FADD, TICAM1, TRADD, MLKL, and CASP10 (Figure 4D and Supplementary Table S13). The significantly enriched 23 KEGG pathways that are associated with cluster 4 are shown in Table 4. The interconnected nodes of cluster 5 (Figure 4E and Supplementary Table S13) and cluster 6 (Figure 4F and Supplementary Table S13) are also shown in Figure 4. The 36 KEGG pathways (such as peroxisome, oxidative phosphorylation, endocytosis, melanogenesis, B cell receptor signaling pathway, Fc epsilon RI signaling pathway, chemokine signaling pathway, Fc gamma R-mediated phagocytosis, T cell receptor signaling pathway, natural killer cell mediated cytotoxicity, acute myeloid leukemia, NOD-like receptor signaling pathway, focal adhesion, PPAR signaling pathway) are significantly enriched in cluster 6 (Table 4).

Figure 4. MCODE-derived significant Clusters of upregulated DEGs in the MDS

Table 4. The enriched KEGG pathways in the significant clusters of upregulated genes in the MDS. The pathways are identified using the GSEA tool.

Name of cluster	KEGG pathways	Genes in Overlap	p-value	FDR value
Cluster 1	RIG-I-like receptor signaling pathway	3	5.72E-05	1.06E-02
	Toll-like receptor signaling pathway	3	1.68E-04	1.56E-02
	Jak-STAT signaling pathway	3	5.73E-04	3.55E-02
Cluster 3	Apoptosis	3	5.33E-06	9.91E-04
	Natural killer cell mediated cytotoxicity	3	2.08E-05	1.94E-03
	Cytokine-cytokine receptor interaction	3	1.48E-04	9.15E-03
	p53 signaling pathway	2	3.31E-04	1.54E-02
Cluster 4	SNARE interactions in vesicular transport	8	8.40E-18	1.56E-15
	RIG-I-like receptor signaling pathway	7	3.11E-13	2.89E-11
	Epithelial cell signaling in Helicobacter pylori infection	6	3.44E-11	2.13E-09
	Apoptosis	6	1.57E-10	7.29E-09
	Adipocytokine signaling pathway	5	3.95E-09	1.47E-07
	Toll-like receptor signaling pathway	5	3.33E-08	1.03E-06
	Vibrio cholerae infection	4	1.64E-07	4.11E-06
	Cytosolic DNA-sensing pathway	4	1.77E-07	4.11E-06
	Oxidative phosphorylation	4	5.96E-06	1.23E-04
	Pathways in cancer	5	1.01E-05	1.89E-04
	Acute myeloid leukemia	3	1.85E-05	3.13E-04
	NOD-like receptor signaling pathway	3	2.39E-05	3.51E-04
	Chemokine signaling pathway	4	2.45E-05	3.51E-04
	B cell receptor signaling pathway	3	4.23E-05	5.61E-04
	MAPK signaling pathway	4	9.38E-05	1.16E-03
	Amyotrophic lateral sclerosis (ALS)	2	1.04E-03	1.21E-02
	Pancreatic cancer	2	1.80E-03	1.97E-02
	Chronic myeloid leukemia	2	1.96E-03	2.02E-02
	Fc epsilon RI signaling pathway	2	2.29E-03	2.24E-02
	Small cell lung cancer	2	2.58E-03	2.40E-02
	Prostate cancer	2	2.89E-03	2.56E-02
	T cell receptor signaling pathway	2	4.22E-03	3.57E-02
	Neurotrophin signaling pathway	2	5.70E-03	4.61E-02
Cluster 6	Aminoacyl-tRNA biosynthesis	6	2.89E-10	2.07E-08
	Peroxisome	7	3.10E-10	2.07E-08
	Parkinson's disease	8	3.34E-10	2.07E-08
	Wnt signaling pathway	7	3.20E-08	1.49E-06
	Oxidative phosphorylation	6	3.62E-07	1.35E-05
	Alzheimer's disease	6	1.39E-06	4.30E-05
	Endocytosis	6	2.29E-06	6.09E-05
	Basal cell carcinoma	4	5.52E-06	1.14E-04
	Pathways in cancer	7	5.52E-06	1.14E-04
	Huntington's disease	5	4.11E-05	7.64E-04
	Melanogenesis	4	6.14E-05	1.04E-03
	Lysine degradation	3	1.09E-04	1.68E-03
	Renal cell carcinoma	3	4.32E-04	6.18E-03
	Chronic myeloid leukemia	3	4.89E-04	6.49E-03
	B cell receptor signaling pathway	3	5.29E-04	6.56E-03
	Fc epsilon RI signaling pathway	3	6.16E-04	7.16E-03
	Chemokine signaling pathway	4	6.74E-04	7.37E-03
	beta-Alanine metabolism	2	9.65E-04	9.44E-03
	Biosynthesis of unsaturated fatty acids	2	9.65E-04	9.44E-03
	Fc gamma R-mediated phagocytosis	3	1.08E-03	1.01E-02
	T cell receptor signaling pathway	3	1.52E-03	1.35E-02
	Citrate cycle (TCA cycle)	2	1.92E-03	1.62E-02
	Propanoate metabolism	2	2.17E-03	1.76E-02
	Natural killer cell mediated cytotoxicity	3	2.99E-03	2.32E-02
	Fatty acid metabolism	2	3.50E-03	2.60E-02
	Jak-STAT signaling pathway	3	4.23E-03	3.01E-02
	Notch signaling pathway	2	4.36E-03	3.01E-02
	Endometrial cancer	2	5.32E-03	3.53E-02
	Non-small cell lung cancer	2	5.72E-03	3.67E-02
	Acute myeloid leukemia	2	6.36E-03	3.94E-02
	Colorectal cancer	2	7.48E-03	4.35E-02
	NOD-like receptor signaling pathway	2	7.48E-03	4.35E-02
	Glioma	2	8.20E-03	4.61E-02
	Focal adhesion	3	8.42E-03	4.61E-02
	PPAR signaling pathway	2	9.19E-03	4.88E-02
	Pancreatic cancer	2	9.45E-03	4.88E-02

In addition, we investigated the significant gene clusters for the downregulated genes. We revealed that one cluster (Supplementary Table S13), with an MCODE score of 13.935, is significantly derived from the original PPI of downregulated genes. The 32 genes are interconnected in the cluster 1 network (Figure 5A and Supplementary Table S13). The highly connected nodes in cluster 1 are EZH2, HELLS, MKI67, CD19, PAX5, RRM2, TCF3, EBF1, CD79A, VPREB1, KIF11, NUSAP1, E2F8, and IL7R. We found that the 3 pathways, included primary immunodeficiency, B cell receptor signaling pathway, and hematopoietic cell lineage are significantly enriched in cluster 1 (Figure 5B).

Figure 5. MCODE-derived significant Cluster of downregulated genes in the MDS (A) and its functional enrichment analysis (B)

3.4 Expression validation of top hub genes in the independent dataset

To determine the aberrant expression of MDS-associated hub genes in RNA-sequence data, we investigated the expression level of top ten upregulated and top ten downregulated hub genes in the independent dataset GSE114922. Interestingly, we found that the upregulated STAT1, IFIH1, RAC2, and CASP1 are significantly upregulated in the GSE114922 when compared to the MDS patients with controls (Figure 6). Besides, the 9 hub genes out of the top ten downregulated hub genes (CREBBP, HIF1A, PIK3CA, EZH2, PIK3R1, MDM2, IRF4, CXCR4, and CD19), are consistently downregulated in the GSE114922 (Figure 6). This result indicates that the aberrant expression of these key hub genes is associated with the pathogenesis of MDS.

Figure 6. Consistent deregulation of hub genes in the independent dataset.

3.5 Hub genes are associated with poor survival prognosis in the MDS

To find the correlation of hub genes with the survival time of MDS patients, we investigated the survival prognosis of 168 upregulated and 52 downregulated hub genes based on the mean expression of individual genes. We found that upregulated hub genes’ expression levels included MX2, GBP2, PXN, IFI44, FDXR, PLCB2, ASS1, ERCC4, PML, and RRAGD are significantly correlated with shorter survival times of MDS patients (Figure 7A). On the other hand, the lower expression level of downregulated hub genes, including CD19, PAX5, TCF3, LEF1, NUSAP1, and TIMELESS are significantly associated with the shorter survival time of MDS patients (Figure 7B). Since, these hub genes are involved in the PPI network with other neighbour genes, indicating that the prognostic hub genes are associated with the clinical outcome of MDS patients through the PPI-based signaling in the cellular environment. Interestingly, the functional enrichment analysis revealed that the sixteen hub genes are related to the enrichment of acute myeloid leukemia, melanogenesis, Wnt signaling pathway, chemokine signaling pathway, and pathways in cancer (P < 0.01). Altogether, our analysis revealed that the expression level of CD34 + cells-derived hub-genes is related to MDS pathogenesis and the shorter survival time of MDS patients.

Figure 7. Identifying prognostic hub genes in the MDS. A. Upregulated 10 hub genes are poorly correlated with the survival prognosis of patients’ B. Downregulated 6 hub genes are poorly associated with a survival prognosis of patients.

3.6 Prognostic Hub genes are associated with the diagnostic value in the MDS

We speculate that the prognostic sixteen hub genes (MX2, GBP2, PXN, IFI44, FDXR, PLCB2, ASS1, ERCC4, PML, RRAGD, CD19, PAX5, TCF3, LEF1, NUSAP1, and TIMELESS) have diagnostic value in MDS patients. We used the GSE58831 dataset (n = 176) to validate our hypothesis, and the results showed that the ROC curve of the expression levels of MX2, GBP2, PXN, IFI44, FDXR, PLCB2, ASS1, ERCC4, PML, RRAGD, CD19, PAX5, TCF3, LEF1, NUSAP1, and TIMELESS hub genes showed excellent diagnostic value for MDS patients (n = 159) when compared with normal samples (n = 17) (Figure 8). The cut-off points, AUC, 95% CI, sensitivity, specificity, PPV, and NPV of these16 genes are illustrated in Table 5. The value of area under the ROC curve (AUC) of PXN, FDXR, PLCB2, PML, CD19, PAX5, and LEF1 genes are more than 0.80, indicating that these genes could be effectively used for diagnostic efficacy in the MDS.

Figure 8. Diagnostic value of prognostic hub genes in the MDS. The MDS patients considering the cases and healthy sample is control.

Table 5. The 16 prognostic genes with cut-off points, AUC, 95% CI, sensitivity, specificity, PPV, and NPV.

Gene	Cut off value	AUC	95% CI	Sensitivity	Specificity	PPV	NPV
MX2	6.937	0.729	62.006-85.089	59.748	94.118	98.958	20
GBP2	9.12	0.62	52.78-72.34	37.107	100	100	14.53
PXN	3.507	0.822	72.545-92.752	88.679	64.706	95.918	37.931
IFI44	7.171	0.761	64.116-88.899	53.459	94.118	98.837	17.778
FDXR	6.937	0.802	69.295-93.117	59.748	94.118	98.958	20
PLCB2	4.325	0.881	82.659-93.442	77.987	94.118	99.2	31.373
ASS1	4.744	0.733	61.812-84.136	61.006	76.471	96.04	17.333
ERCC4	4.188	0.682	54.188-80.44	49.057	82.353	96.296	14.737
PML	3.739	0.8	69.365-91.789	86.164	64.706	95.804	33.333
RRAGD	3.739	0.71	59.235-84.679	86.164	64.706	95.804	33.333
CD19	3.06	0.855	69.516-91.12	84.277	82.353	97.81	35.897
PAX5	4.91	0.859	81.346-92.164	76.73	100	100	31.481
TCF3	5.42	0.725	58.97-86.794	80.503	64.706	95.522	26.19
LEF1	3.872	0.825	76.558-90.145	64.78	100	100	23.288
NUSAP1	9.221	0.652	56.446-75.037	47.799	94.118	98.701	16.162
TIMELESS	5.581	0.681	57.887-80.7	51.572	82.353	96.471	15.385

AUC is area under the ROC curve, PPV is positive predictive value, and NPV is negative predictive value.

4. Discussions

Genome-wide expression profiling is the fundamental step to identify the differentially expressed genes in disease conditions [24]. We combine the four microarray gene expression profiling datasets into a single integrated dataset for investigating the differential expression of genes in the MDS relative to the healthy controls. We revealed that the 1619 upregulated (Supplementary Table S1) and 785 downregulated genes (Supplementary Table S2) in the MDS patients by a meta-analysis (Combined ES > 0.585, Adjusted P < 0.05). Previous studies consistently explored that CD34 + cells-derived transcriptomes are associated with MDS and acute myeloid leukemia. For example, the expression level of the ALDH3B1 gene is elevated in the MDS [10]. E2F6, another upregulated gene, is associated with CHK proteins in cell-cycle checkpoint control and consistently deregulated in the MDS [10]. Deregulated expression of IRF4 is associated with acute myeloid leukemia and myelodysplastic syndrome susceptibility [25]. AKAP12, a downregulated gene with tumour suppressor properties, is related to the myeloid malignancies [26].

It was previously reported that the functional enrichment analysis, including gene ontology and the enrichment of pathways, are substantially associated with significant DEGs in the MDS [5,8]. We found that the upregulated and downregulated DEGs are significantly related to various biological processes, cellular components, and molecular functions, respectively (Figure 1). Since deregulated pathways are critically associated with the molecular pathogenesis of MDS [27], we identified the KEGG pathways that are significantly associated with identified DEGs. The upregulated pathways are mainly involved in immune system deregulation, metabolic abnormalities, and aberrant cellular signalling and developmental process (Figure 2A and Supplementary Table S9). In addition, the downregulated genes are aberrantly involved with the deregulation of immunity, cancer, and cellular signalling in the MDS (Figure 2A and Supplementary Table S10). A. Pellagatti et al. demonstrated that the Wnt pathways, immunodeficiency, apoptosis, and chemokine signaling are the deregulated pathways in the MDS [10]. Ying le et al. reported that hematopoietic cell lineage is deregulated in the MDS [5]. B cell receptor signaling pathway is associated with shared gene signature between MDS and AML [28]. Altogether, these results indicated that the aberrantly expressed genes are involved with the deregulation of signaling pathways and key biological functions in the MDS.

Various human diseases are originated from the result of abnormal protein–protein interactions [29]. We constructed and identified the key interacting genes from the PPI network of DEGs (Figure 3). We found that the STAT1, STAT3, IFIH1, EPRS, GRB2, RAC2, MAPK14, CASP1, and SPI1 were the top upregulated hub genes (Supplementary Table S11) and CREBBP, HIF1A, PIK3CA, EZH2, PIK3R1, MDM2, IRF4, CXCR4, PCNA, and CD19 (Supplementary Table S12) were the top downregulated hub genes in the MDS. Digdem Yoyen-Ermis et al. demonstrated that the STAT3 pathway is activated in the AML and MDS [30]. In stem cells of MDS patients, the master myeloid differentiation – driving transcription factor SPI1 was expressed at higher levels [31]. In our expression validation analysis, we checked the expression differences of the top 10 upregulated and top 10 downregulated hub genes in the independent dataset (GSE114922). Interestingly, we revealed that the thirteen genes are consistently deregulated in the validation cohort (Figure 6), indicating that the expression of these genes is consistently deregulated in the multiple datasets. Some previous studies reported the consistent results. For example, the elevated IFN-γ/STAT1 signaling is associated with MDS progression and treatment response [32]. Rac2 may play a crucial role in the initial development of myelodysplastic syndrome [33]. The expression of CREBBP is associated with hematopoietic malignancies, including myelodysplastic syndrome (MDS) [34]. Variants in IRF4 related to acute myeloid leukemia and MDS susceptibility [25]. Anke C. Spoo et al. suggested that CXCR4 should be incorporated into the risk assessment of AML patients [35]. The expression level of CD19 is downregulated and one of the crucial genes in the MDS [36]. Since genes with similar characteristics have tended to interact with each other through PPI networks in the diseases/disorders, clustering analysis provides the information to predict human disease-related gene clusters/modules [37]. We investigated the significant gene clusters (MCODE score > 5.0) in the MDS, and we revealed that these gene clusters are substantially related to the enrichment of KEGG pathways in the MDS (Figures 4 and 5 and Table 4), indicating that these gene clusters are crucially associated with the pathogenesis of MDS.

Furthermore, we analyzed the association of hub genes with the survival prognosis of MDS patients. The upregulation of MX2, GBP2, PXN, IFI44, FDXR, PLCB2, ASS1, ERCC4, PML, and RRAGD hub genes are significantly correlated with shorter survival times of MDS patients (Figure 7A) and the downregulation of CD19, PAX5, TCF3, LEF1, NUSAP1, and TIMELESS (Figure 7B) are significantly correlated with the overall survival prognosis of MDS patients. It was stated that the overexpression of PXN is associated with poor prognosis in AML patients [38]. M. A. ElBaiomy et al. concluded that the higher expression level of LEF1 is a favourable prognostic marker in the AML patient risk [39]. Finally, we found that these survival-associated hub genes have diagnostic efficacy in the MDS when compared with healthy control (Figure 8). Altogether, this study provides the clue that these deregulated sixteen genes (MX2, GBP2, PXN, IFI44, FDXR, PLCB2, ASS1, ERCC4, PML, RRAGD, CD19, PAX5, TCF3, LEF1, NUSAP1, and TIMELESS) may be associated with the prognostic and diagnostic value of MDS patients. The aberrant expression level of CD19 is indicated as the marker of monocytes lineage in acute myelogenous leukemia [40]. In hematolymphoid, it was demonstrated that the lower expression of PAX-5 helps to establish B-cell lineage and increases diagnostic yield [41]. In the myelodysplastic syndromes, downregulation of the LEF1 is associated with disease progression [42].

The major limitation of our study is that the deregulated key genes and networks identified by bioinformatics analysis have not been validated by experimental research in the MDS patients. Another drawback of our study is that MDS is a group of heterogonous diseases, the different classification subgroup and the different prognosis subgroup should have other hub genes. However, we did not identify the hub genes in the different subtypes of MDS patients. Thus, although our findings could provide potential biomarkers for the MDS prognosis and diagnosis, as well as therapeutic targets, further experimental and clinical validation would be warranted to transform these findings into the clinical application of MDS and different subtypes of MDS patients.

5. Conclusions

This bioinformatics analysis provided compelling evidence about the deregulated CD34 + cells-derived hub transcriptomes that are associated with an inferior prognosis and diagnostic efficacy in the MDS patients. The association of these prognostic-hub genes should be experimentally validated in different subtypes of MDS patients. In summary, identifying key hub genes and their association with the prognostic and diagnostic efficacy may provide substantial clues for the diagnosis and treatment of MDS patients.

Authors’ contributions

All authors have contributed to the research concept, design, analysis, and interpretation of the results. All authors read, reviewed, and approved the document's final version.

Disclosure statement

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

Data availability

The GEO datasets (GSE4619 (n = 66), GSE19429 (n = 200), GSE30195 (n = 19), GSE58831 (n = 176), GSE114922) used in this study were downloaded from The National Biotechnology Information Center Gene (NCBI-Gene) database (https://www.ncbi.nlm.nih.gov/gene).

Additional information

Funding

This work was supported by National Natural Science Foundation of China: [Grant Number 81560027].