Screening and identification of key candidate genes and pathways in myelodysplastic syndrome by bioinformatic analysis

Introduction

Myelodysplastic syndrome (MDS) is a clonal hematopoietic stem cell (HSC) disease, mainly involving cytogenetic changes and/or genetic mutations, and there is also a widespread gene hypermethylation at advanced stage (Adès, Itzykson & Fenaux, 2014). The main features of MDS are myeloid cell cytopenias, morphologic dysplasia, ineffective hematopoiesis and a high risk of transformation to acute myeloid leukemia (AML) (Weinberg & Hasserjian, 2019). Although there are many new drugs for the treatment of MDS in recent years, about one-third of patients with MDS experience transformation to AML and the overall survival (OS) of MDS remains not ideal (Duchmann & Itzykson, 2019). There is accumulating evidence that MDS are associated with abnormal genetic mutations (Ogawa, 2019). Therefore, identifying more molecular biomarkers is critical for early diagnosis, treatment and prognosis of MDS.

In recent years, gene chip technology has developed rapidly and has been widely used in gene detection (Vogelstein et al., 2013). It has been confirmed that the pathogenesis and progression and heterogeneity of MDS are closely related to the genetic landscape (Maciejewski & Balasubramanian, 2017), such as tet methylcytosine dioxygenase 2 (TET2), isocitrate dehydrogenase1/2 (IDH1/2), additional sex comb-like 1 (ASXL1), enhancer of zeste homolog 2 (EZH2), NRAS proto-oncogene (NRAS), RUNX family transcription factor 1 (RUNX1) and GATA binding protein 2 (GATA2). TET2 is the most common mutations in MDS and it is reported that TET2 deficiency disturbed erythroid differentiation, contributing to ineffective erythropoiesis (Qu et al., 2018). IDH1/2 mutation was positively correlated with bone marrow blast counts and negatively correlated with absolute neutrophil counts. Patients with IDH1 mutations had shorter OS and progression-free survival (PFS), while IDH2 mutations had no effect on OS or PFS of MDS patients (Wang et al., 2017). The ASXL1 mutation can significantly and independently predict adverse outcomes of myeloid tumors (Mangaonkar et al., 2018). Aberrant EZH2 splicing and/or EZH2 mutations which impaired histone H3 lysine 27 tri-methylation may be associated with the pathogenesis of MDS (Shirahata-Adachi et al., 2017). NRAS was thought to be significantly associated with MDS, especially in the progression to leukemia (Leite et al., 2017; Makishima, 2019). RUNX1 mutation was observed in approximately 10% of patients with MDS. Both the repressor and activator of RUNX1 resulted in the tumor suppressor activity in MDS and AML (Behrens et al., 2016). GATA2 has been shown to be related to the progression of MDS, such as the progression from low-risk to high-risk or AML (Makishima, 2019; McReynolds et al., 2019). Although the discovery of these genes contributes to the research of MDS, the pathogenesis of MDS remains unclear. Re-integration and distribution of these abnormal genes provides the possibility to identify diseases associated with biological markers. The bioinformatics approach provides a good platform to correlate gene expression profiles with disease.

In this study, I downloaded two microarray datasets from the Gene Expression Omnibus (GEO) database to obtain differentially expressed genes (DEGs) from bone marrow CD34+ cells of MDS and healthy controls. Then, I performed Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and protein–protein interaction (PPI) network analyses on these DEGs. Finally, I screened the module genes and the hub genes and analyzed the relationship between gene expression and prognosis to help understand the underlying molecular mechanism of MDS.

Materials and Methods

Results

Identification of DEGs in MDS

I obtained gene expression profiles from the GSE19429 to GSE58831 datasets of CD34+ cells of bone marrow samples from MDS patients and healthy controls and analyzed DEGs using GEO2R. As shown in Figs. 1A and 1B, I identified 179 and 500 DEGs from GSE19429 to GSE58831, respectively. The overlap among the two datasets between MDS patients and healthy control contained 133 genes as shown in the Venn diagram (Fig. 1C).

GO and KEGG pathway enrichment analyses of DEGs

I used DAVID to perform and analyze the biological classification of DEGs, which included GO analysis and KEGG pathway enrichment analysis. GO analysis showed that the function of upregulated genes in DEGs were mainly enriched in type I interferon signaling pathway, defense response, response to external stimulus, immune effector, cell surface receptor signaling pathway, cytokine-mediated signaling pathway, regulation of cell action potential, sensory perception, immune response, Notch signaling pathway (Table 1). Downregulated genes in BP were mainly enriched in immune response, cell surface receptor signaling pathway, cell activation, cell differentiation, regulation of immune system, regulation of macromolecule metabolic, cell-cell adhesion, biosynthetic process, cell proliferation, cell death, specific DNA binding, apoptotic process, regulation of metabolic process, epithelium development, blood vessel morphogenesis, cell migration, muscle tissue development. For CC, downregulated genes were mainly enriched in side of membrane and extracellular space. For MF, enrichment of downregulated genes was primarily in transcription factor activity, receptor binding, signal transducer activity, steroid hormone receptor activity, hormone receptor binding, core promoter binding, macromolecular complex (Table 1). KEGG pathway analysis revealed that enrichment of DEGs was mostly in primary immunodeficiency, hematopoietic cell lineage, FOXO signaling pathway, hippo signaling pathway, transcriptional misregulation in cancer (Fig. 2).

Table 1:

Gene Ontology term enrichment analysis of DEGs in MDS.

Term	Description	Count	P-value
Upregulated
GO:0060337	Type I interferon signaling pathway	4	2.28E−04
GO:0098542	Defense response	5	4.77E−04
GO:0009605	Response to external stimulus	10	4.04E−03
GO:0002252	Immune effector	6	6.19E−03
GO:0007166	Cell surface receptor signaling pathway	11	6.70E−03
GO:0019221	Cytokine-mediated signaling pathway	5	1.16E−02
GO:0098911	Regulation of cell action potential	2	2.25E−02
GO:0050954	Sensory perception	3	2.69E−02
GO:0006955	Immune response	7	3.50E−02
GO:0045746	Notch signaling pathway	2	4.76E−02
Downregulated
GO:0006955	Immune response	27	9.26E−09
GO:0007166	Cell surface receptor signaling pathway	33	3.21E−07
GO:0000982	Transcription factor activity	11	4.51E−06
GO:0001775	Cell activation	16	2.84E−05
GO:0040011	Locomotion	20	1.06E−04
GO:0005102	Receptor binding	19	1.14E−04
GO:0045595	Cell differentiation	20	1.20E−04
GO:0002684	Regulation of immune system	15	1.99E−04
GO:0010605	Regulation of macromolecule metabolic	24	5.27E−04
GO:0022407	Cell-cell adhesion	9	5.71E−04
GO:0009891	Biosynthetic process	20	6.39E−04
GO:0004871	Signal transducer activity	19	1.12E−03
GO:0008283	Cell proliferation	20	1.26E−03
GO:0060548	Cell death	13	1.27E−03
GO:0098552	Side of membrane	9	1.54E−03
GO:0043565	Specific DNA binding	14	1.47E−03
GO:0006915	Apoptotic process	19	1.57E−03
GO:0009892	Regulation of metabolic process	24	1.63E−03
GO:0060429	Epithelium development	14	1.72E−03
GO:0048514	Blood vessel morphogenesis	9	2.64E−03
GO:0003707	Steroid hormone receptor activity	4	2.71E−03
GO:0016477	Cell migration	14	4.50E−03
GO:0060537	Muscle tissue development	7	8.93E−03
GO:0005615	Extracellular space	14	2.40E−03
GO:0051427	Hormone receptor binding	4	3.39E−02
GO:0001047	Core promoter binding	4	4.05E−02
GO:0044877	Macromolecular complex	12	4.16E−02

DOI: 10.7717/peerj.8162/table-1

PPI network establishment and module analysis

The PPI network of DEGs was established based on the STRING online database and Cytoscape software (15 upregulated and 103 downregulated) (Fig. 3A). There were 133 DEGs in the two datasets, however, 15 intersection genes were not found in the STRING. We obtained four significant modules based on the degree of importance by utilizing cluster analysis of the PPI network in Cytoscape MCODE. Module 1 contained 12 nodes and 53 edges (Fig. 3B); Module 2 contained 6 nodes and 15 edges (Fig. 3C); Module 3 contained 4 nodes and 5 edges (Fig. 3D); Module 4 contained 3 nodes and 3 edges (Fig. 3E). The enrichment analyses of genes involved in these modules were preformed using DAVID. As showed in Table 2, GO analysis indicated these genes in modules were mainly concentrated in immune response, immune cell differentiation, side of membrane, transcription factor, DNA binding, cell adhesion, hemopoiesis, cell activation, defense response, cell surface receptor, DNA metabolic, kidney development. As for KEGG pathway analyses, module genes mainly enriched in primary immunodeficiency, hematopoietic cell lineage, B cell receptor signaling pathway, FOXO signaling pathway.

Table 2:

Gene Ontology and KEGG pathway enrichment analysis of DEGs in the significant module.

Term	Description	Count	P-value
GO:0006955	Immune response	10	3.15E−06
GO:0030098	Immune cell differentiation	6	4.93E−05
GO:0098552	Side of membrane	5	5.71E−04
GO:0000982	Transcription factor	5	4.34E−04
GO:0003677	DNA binding	9	5.97E−04
GO:0016337	Cell adhesion	6	1.24E−03
GO:0030097	Hemopoiesis	6	2.09E−03
GO:0001775	Cell activation	6	2.34E−03
GO:0006952	Defense response	7	2.99E−03
GO:0007166	Cell surface receptor	8	1.22E−02
GO:0006259	DNA metabolic	5	1.83E−02
GO:0001822	Kidney development	3	3.90E−02
Hsa05340	Primary immunodeficiency	5	3.02E−07
Hsa04640	Hematopoietic cell lineage	5	1.36E−05
Hsa04662	B cell receptor signaling pathway	3	6.19E−03
Hsa04068	FOXO signaling pathway	3	2.09E−02

DOI: 10.7717/peerj.8162/table-2

Hub genes selection and biological analysis

In study, a total of four hub genes were identified by MCODE, as IFI44L, IRF4, IL7R, SKIL. I constructed a network of hub genes and their co-expression genes via cBioPortal, but did not find the co-expressed gene of IFI44L in cBioPortal (Fig. 4A). Then I analyzed the BPs of hub genes through BINGO as showed in Fig. 4B. Hierarchical clustering of hub genes was constructed using UCSC and blood derived cancers and healthy controls were roughly identified by hub genes (Fig. 4C).

Survival analysis

Kaplan–Meier curve was employed to predict the prognosis of the 25 identified module genes. First I had a simple understanding of these module genes, such as the full name, abbreviation and function of the genes (Table 3). Among the genes examined, MDS patients having BLNK, IRF4, IFITM1, IFIT1, ISG20, IFI44L alterations showed worse OS (Fig. 5). Although the difference between IFIT1 and IFI44L was not statistically significant, the OS of patients having IFIT1 and IFI44L alteration showed an obvious decline tendency.

Table 3:

Functional roles of 25 module genes.

No.	Gene symbol	Full name	Function
1	RAG2	Recombination activating 2	RAG2 is involved in B and T cell development and the expression of RAG2 is high in most AMLpatients
2	EBF1	Transcription factor, early B cell factor 1	EBF1 promotes of bone marrow CD34+ cell proliferation and migration,but suppresses the apoptosis of cells
3	BLNK	B cell linker	BLNK regulates the development, maturation and function of B cells
4	CD79B	CD79b molecule	CD79B is necessary for expression and function of the B-cell antigen receptor
5	VPREB1	V-set pre-B cell surrogate light chain 1	VPREB1 is involved in early B cell differentiation and VPREB1 deletion is often detected in relapsed leukemia
6	BLK	BLK proto-oncogene, Src family tyrosine kinase	BLK is involved in cell proliferation and differentiation and has a role in B-cell receptor signaling and B-cell development
7	IRF4	Interferon regulatory factor 4	IRF4 has oncogenic implications in a variety of hematological tumors including multiple myeloma, leukemia and lymphoma
8	IGLL1	Immunoglobulin lambda like polypeptide 1	IGLL1 is involved in cell proliferation and IGLL1 mutation can lead to B cell deficiency and agammaglobulinemia
9	DNTT	DNA nucleotidy-lexotransferase	DNTT acts as a catalyst to add nucleotides during the maturation of B cells and T cells
10	CD19	CD19 molecule	CD19 is a surface-specific marker of B lymphocytes and CD19 decreases the threshold of antigen receptor-dependent
11	PAX5	Paired box 5	PAX5 is involved in cell differentiation and neuro developmentand spermatogenesis
12	RAG1	Recombination activating 1	High expression of RAG1 in B-ALL patients promotes leukemic clonal evolution and high expression of RAG1 is associated with high proliferative markers
13	IFITM1	Interferon induced transmembrane protein 1	IFITM1 promotes proliferation, migration, and invasion of lung cancer
14	IFIT1	Interferon induced protein with tetratricopeptide repeats 1	IFIT1 and IFIT3 promote oral squamous cell carcinoma metastasis and contribute to theanti-tumor effect of gefitinib
15	IFIT3	Interferon induced protein with tetratricopeptide repeats 3	IFIT1 and IFIT3 promote oral squamous cell carcinoma metastasis and contribute to the anti-tumor effect of gefitinib
16	ISG20	Interferon stimulated exonuclease gene 20	ISG20 promotes local tumor immunity and has poor survival in human glioma
17	IFI44L	Interferon induced protein 44 like	High expression of IFI44L is associated with poor prognosis of osteosarcoma
18	IFI27	Interferon alpha inducible protein 27	IFI27 promotes cell proliferation and invasion in oral squamous cell carcinoma
19	CD24	CD24 molecule	CD24 modulates growth and differentiation signals to mature granulocytes and B cells
20	CXCR4	C-X-C motif chemokine receptor 4	High expression of CXCR4 has shorter OS time and shorter relapse-free survival time
21	IL7R	Interleukin 7 receptor	IL7R is overexpressed in adult acute lymphoblastic leukemia and is associatedwith JAK/STAT pathway mutations
22	MME	Membrane metalloendo-peptidase	MME inhibits tumor metastasis of esophageal squamous cell carcinoma and interrupt tumor cell adhesion
23	SKIL	SKI like proto-oncogene	SKIL has regulatory role in cell division or differentiation
24	SMAD7	SMAD family member 7	High expression of SMAD7 at diagnosis predicts poor prognosis in acute myeloid leukemia
25	SMAD1	SMAD family member 1	SMAD1 is involved in a range of biological activities including cell growth, apoptosis, morphogenesis, development and immune responses

DOI: 10.7717/peerj.8162/table-3

Discussion

MDS is a type of stem cell-derived bone marrow clonal disease with a relatively uneven presentation profile. The incidence of MDS increases with age, the onset of MDS before the age of 40 years is rare. However, when the age is older than 70, the incidence rate exceeds 50 per 100,000 people. In recent years, due to the gradual increase of the population aging problem, the number of patients with MDS has also increased year by year. The exact cause of MDS is uncertain and its relatively risk factors include chemotherapy or radiation therapy and tobacco or poisonous solvents or agricultural chemicals (Tefferi & Vardiman, 2009). In recent years, molecular biology and immune regulation have become increasingly important in MDS (Hosono, 2019; Liu et al., 2019). MDS, also known as pre-leukemia, is characterized by a high risk of transformation to AML. This may be one of the causes of poor prognosis in MDS, so identifying potential markers for MDS diagnosis and prognosis is extremely urgent.

Recently, gene chips have been widely used in research and the corresponding data have been stored in public databases. A number of studies have been conducted based on microarray data profiles to elucidate the mechanisms of disease. But if there is only one data, errors may occur. Therefore, combining bioinformatics methods with expression profiling techniques can partially reduce bias. In this study, I analyzed the expression of genes in two mRNA microarray datasets based on CD34+ cells of bone marrow samples from MDS patients and healthy controls. A total of 133 DEGs were identified, including 15 up-regulated genes and 103 down-regulated genes, however, 15 DEGs were not found in STRING. GO and KEGG pathway enrichment analysis was performed by DAVID online software to explore the biological functions of DEGs. GO enrichment analysis indicated that the identified DEGs were mainly enriched in the immune process. Cancer cells involved in dynamic interference of immune response, evading immune destruction is one of the main features of cancer, involving the participation of a variety of molecular events (Grivennikov, Greten & Karin, 2010; Hanahan & Weinberg, 2011). In a mouse model, IFN-γ-induced apoptosis gene expression promotes destruction of HSCs and progenitor cells (Chen et al., 2015). Interferons are proteins secreted by human immune cells in response to infections (mainly viral). Most of the genes I have identified are immune genes and/or related to the interferon pathway. Therefore, I speculate that the pathogenesis of MDS may be related to immune factors, similar to the study by Silzle et al. (2019) which showed that lymphopenia has a significantly reduced OS in low-risk MDS patients, these data subtly support the role of the immune system in MDS.

In the current study, I selected IFI44L, IRF4, IL7R, SKIL4 as hub genes with the highest degrees. IFI44L, IRF4 and IL7R are immune genes and/or related to the interferon pathway. IFI44L, a type I interferon-stimulated gene, has been shown to be involved in human carcinogenesis and tumor progression (Gao et al., 2013). Pellagatti et al. (2018, 2010) showed that IFI44L was upregulated by more than twofold in at least 60% of the MDS patients and indicated that IFI44L is one of the significant survival predictors in MDS. A recent study found that down-regulation of IRF4 is a necessary and independent effector of immunomodulatory drugs in primary effusion lymphoma (Patil, Manzano & Gottwein, 2018). IRF4 has oncogenic implications in a variety of hematological tumors including multiple myeloma, leukemia and lymphoma (Agnarelli, Chevassut & Mancini, 2018; Cherian et al., 2018; So et al., 2014; Zhang et al., 2013). IL7R is over expressed in adult acute lymphoblastic leukemia (ALL) and is associated with JAK/STAT pathway mutations (Gianfelici et al., 2019). In contrast, IL7R is down-regulated in MDS and this differential expression may be contributed to identifying MDS subtypes and leukemia (Pellagatti et al., 2004). SKIL, also known as SNON, promotes cell invasion of immortalized human mammary epithelial cells. These differential expression could be either cause of the ailment or an effect (or an indirect effect).

The PPI network of DEGs was established and four modules were considered significant in this network. A total of 25 genes were identified in the module models and the functions of these genes were analyzed and Kaplan–Meier curves were employed to predict prognosis. Moreover, the BLNK, IRF4, IFITM1 and ISG20 alteration are significantly associated with worse OS. While the IFIT1 and IFI44L alteration were not statistically significant, they showed an obvious downward tendency. IFI44L and IRF4 have been introduced previously. BLNK (B Cell Linker Protein), which is a selective target of repression, can result in the development of ALL by differentiation block of PAX5-PML protein (Imoto et al., 2016). The results of Nakayama et al. (2009) are consistent with this and somatic cell loss and mutation of BLNK lead to the generation of pre-B-cell leukemia. IFIT1 and IFITM1 are up-regulated in the CD34+ cells of most patients with MDS. IFIT1 and IFITM1 are interferon-stimulated genes (ISG) and up-regulation of IFIT1 is a potential diagnostic marker for MDS (Pellagatti et al., 2006). The ISG20 gene was first introduced by Gongora et al. (1997) as a promyelocytic leukemia nuclear body (PML-NB)-associated protein. The expression of ISG20 is reduced in IDH mutant tumors and ISG20 can promote local tumor immunity and may serve as a potential therapeutic target of human glioma (Gao et al., 2019). This is very interesting: all these identified genes are somehow related to interferons and/or immune response. However, due to a shortfall of the underlying data acquisition technology, false positive results may occur, which is also the biggest flaw of this paper. Next, I plan to collect samples of MDS patients and verify the conclusions of this paper through experiments.

Conclusion

In conclusion, this study aimed to perform a comprehensive bioinformatics analysis of datasets for MDS patients and healthy controls. Identify DEGs that may be involved in MDS progression and prognosis. Interestingly, most of prognostic genes and hub genes are interferons and immune-related genes. However, further studies are warranted to elucidate the biological function of these genes in MDS.

Supplemental Information