Identification of Key Genes and Pathways Associated with RUNX1 Mutations in Acute Myeloid Leukemia Using Bioinformatics Analysis.

BACKGROUND RUNXl plays a key regulatory role in the process of hematopoiesis and is a common target for multiple chromosomal translocations in human acute leukemia. Mutations of RUNX1 gene can lead to acute leukemia and affect the prognosis of AML patients. We aimed to identify pivotal genes and pathways involved in RUNX1-mutated patients of with acute myeloid leukemia (AML) and to explore possible molecular markers for novel therapeutic targets of the disease. MATERIAL AND METHODS The RNA sequencing datasets of 151 cases of AML were obtained from the Cancer Genome Atlas database. Differentially expressed genes (DEGs) were identified using edgeR of the R platform. PPI (protein-protein interaction) network clustering modules were analyzed with ClusterONE, and the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses for modules were performed. RESULTS A total of 379 genes were identified as DEGs. The KEGG enrichment analysis of DEGs showed significantly enriched pathways in cancer, extracellular matrix (ECM)-receptor interaction pathway, and cyclic adenosine monophosphate (cAMP) signaling pathway. The top 10 genes ranked by degree were PRKACG, ANKRD7, RNFL7, ROPN11, TEX14, PRMT8, OTOA, CFAP99, NRXN1, and DMRT1, which were identified as hub genes from the protein-protein interaction network (PPI). Statistical analysis revealed that RUNX1-mutated patients with AML had a shorter median survival time (MST) with poor clinical outcome and an increased risk of death when compared with those without RUNX1 mutations. CONCLUSIONS DEGs and pathways identified in the present study will help understand the molecular mechanisms underlying RUNX1 mutations in AML and develop effective therapeutic strategies for RUNX1-mutation AML.

Background

Acute myeloid leukemia (AML) is the most common type of adult leukemia and its incidence rate is increasing every year [1]. The treatment has advanced greatly in recent years, resulting in an increased remission rate of 50% to 80%. However, there are still many patients who do not respond well to multiple-induction chemotherapy and the early recurrence rate is still high. The 5-year overall survival (OS) is still at a low level of 30% to 40%, while recurrent and refractory patients are less than 15% [1,2], and 60% to 80% of patients with complete remission (CR) will eventually experience recurrence. Such recurrence is more likely to cause resistance to chemotherapy and failed remission of leukemia stem cells (LSCs) [3,4].

A variety of gene mutations have been found in AML patients. Runt-related transcription factor 1 (RUNXl) is a member of the RUNX transcription factor family and plays an important role in the determination of cell lineage differentiation, normal hematopoietic cell formation, and stem cell proliferation [5,6]. Many studies have shown that RUNX1 is involved in the differentiation of myeloid cells, B cells, and T cells. Recent studies have found that RUNXl also promotes leukemia cell proliferation, suggesting that RUNXl plays a different role in different hematological malignancies [7-9]. As a transcription factor, RUNXl regulates the expression of many hematopoietic-related genes that control hematopoietic cell differentiation, apoptosis, and self-renewal. The mutations of RUNXl often lead to acute leukemia [6] and seriously affect disease outcome [10].

The aim of the present study was to find the key genes and pathways involved in RUNX1 mutations of AML.

Material and Methods

RNA sequencing data

Identification of differentially expressed genes (DEGs)

EdgeR was used according to the user’s guide for screening differential expression of genes at gene levels [14,15]. We identified differentially expressed genes (DEGs) with fold change (FC) ≥2, and defined P value and false discovery rate (FDR) cut-offs of <0.05 to be statistically significant. A heat map and volcano plot of the DEGs were drawn using the ggplot package in the R platform. In order to annotate input genes, classify gene functions, identify gene conversions, and carry out gene ontology (GO) term analysis, we use integrated Discovery (DAVID) v. 6.8 (; accessed as in August 16, 2017) [16]. To analyze the DEGs at the functional level, GO enrichment and KEGG pathway analysis were performed using the DAVID online tool. P value <0.05 was considered statistically significant.

Integration of protein–protein interaction (PPI)

The Search Tool for the Retrieval of Interacting Genes (STRING) database (; accessed August 16, 2017) is an online tool designed to evaluate the protein–protein interaction (PPI) information [17]. To evaluate the interactive relationships among DEGs, we mapped the DEGs to STRING in order to evaluate the interactive relationships among DEGs. Experimentally validated interactions with a combined score >0.4 were selected as significant. Using the PPI networks, module screening was performed using Molecular Complex Detection (MCODE) (scores >3 and nodes >4) in Cytoscape, an integrated bioinformatics platform [18].

Statistical analysis

We conducted the statistical analyses with SPSS version 20.0 and R 3.3.0.

Using the Cox proportional hazards regression model, hazard ratios (HRs) and their 95% confidence intervals (CIs) were calculated. FDR in edgeR was adjusted for multiple testing with the Benjamini-Hochberg procedure to control FDR [19-21]. A value of P<0.05 was considered as statistically significant.

Results

Gene expression dataset

Information for 151 patients with adult AML and a corresponding bone marrow RNA-Seq dataset were obtained from the TCGA database. These data were obtained from the cBioPortal for Cancer Genomics website.

GO and KEGG

Altogether, we presented 379 DEGs (Figure 1) for further GO and KEGG pathway analyses using DAVID (Tables 1, 2). Heat maps of potential RUNX1-mutation-related DEGs are shown in Figure 2. The KEGG enrichment analysis of DEGs showed that these DEGs were significantly enriched in cancer pathways, as well as in the extracellular matrix-receptor interaction pathway and the cAMP signaling pathway.

Figure 1

Volcano plot of the differentially expressed genes. Red: upregulation; green: downregulation; black: non-differentially expressed genes. P-adj, adjusted P value.

Table 1

GO term analysis of RUNX1 mutation-DEGs associated with AML.

Category	GO ID	Term	Count	%	P-value
GOTERM_BP_DIRECT	GO: 0007155	Cell adhesion	20	3.86	0.001790
GOTERM_BP_DIRECT	GO: 0001525	Angiogenesis	12	2.32	0.004496
GOTERM_BP_DIRECT	GO: 0051897	Positive regulation of protein kinase B signaling	7	1.35	0.006020
GOTERM_BP_DIRECT	GO: 0035360	Positive regulation of peroxisome proliferator activated receptor signaling pathway	3	0.58	0.007486
GOTERM_BP_DIRECT	GO: 0043066	Negative regulation of apoptotic process	17	3.28	0.017451
GOTERM_BP_DIRECT	GO: 0071356	Cellular response to tumor necrosis factor	7	1.35	0.020903
GOTERM_BP_DIRECT	GO: 0050679	Positive regulation of epithelial cell proliferation	5	0.97	0.029629
GOTERM_BP_DIRECT	GO: 2000147	Positive regulation of cell motility	3	0.58	0.029648
GOTERM_BP_DIRECT	GO: 0050930	Induction of positive chemotaxis	3	0.58	0.033775
GOTERM_BP_DIRECT	GO: 0072604	Interleukin-6 secretion	2	0.39	0.038175
GOTERM_BP_DIRECT	GO: 0007568	Aging	8	1.55	0.043500
GOTERM_BP_DIRECT	GO: 0071300	Cellular response to retinoic acid	5	0.97	0.047991
GOTERM_BP_DIRECT	GO: 0005912	Adherens junction	6	1.16	0.002782
GOTERM_BP_DIRECT	GO: 0008083	Growth factor activity	8	1.55	0.031109
GOTERM_BP_DIRECT	GO: 0005112	Notch binding	3	0.58	0.042876

GO – gene ontology; DEGs – differentially expressed genes; AML – acute myeloid leukemia.

Table 2

KEGG pathway analysis of RUNX1mutation-DEGs associated with AML.

Pathway ID	Name	Count	%	P-value	Genes
hsa05200	Pathways in cancer	20	3.87	5.14E-05	PTGER1, PTGER3, LPAR4, RUNX1T1, FGF10, GNG11, MECOM, GLI2, GLI3, MMP2, CTNNA3, PRKACG, LAMB4, WNT7B, PTK2, LAMA3, LAMC3, NTRK1, NOS2, WNT8A
hsa04080	Neuroactive ligand-receptor interaction	12	2.32	0.009536961	LEP, CRHR2, HCRTR1, PTGER1, GABRR3, GABRE, GLRB, PTGER3, PRSS2, LPAR4, CHRNA6, CTSG
hsa05146	Amoebiasis	7	1.35	0.01088093	PRKACG, LAMB4, PTK2, LAMA3, LAMC3, NOS2, CTSG
hsa04512	ECM-receptor interaction	6	1.16	0.018487893	LAMB4, LAMA3, LAMC3, COL6A5, SV2B, TNN
hsa05032	Morphine addiction	6	1.15	0.022018528	PRKACG, GABRR3, GABRE, GNG11, PDE4C, KCNJ3
hsa04024	cAMP signaling pathway	9	1.74	0.023512773	PRKACG, PTGER3, NPY, ATP1B4, RYR2, PDE4C, CREB3L3, TNNI3, GLI3
hsa05205	Proteoglycans in cancer	9	1.73	0.024791345	PRKACG, WNT7B, PTK2, MRAS, IGF2, FLNC, MMP2, WNT8A, TWIST1
hsa04726	Serotonergic synapse	6	1.16	0.046110307	PRKACG, ALOX15B, MAOA, GNG11, KCNJ3, PLA2G4D

KEGG – Kyoto Encyclopedia of Genes and Genomes; DEGs – differentially expressed genes; cAMP – cyclic adenosine monophosphate.

Figure 2

Heat map of the RUNX1 mutation -related differentially expressed genes. Red: upregulation; green: downregulation.

Networks and module analysis

We constructed the protein–protein interactome networks and identified some RUNX1 mutation-associated hub genes (Figure 3). The top 10 genes ranked by degree were identified as hub genes, including PRKACG, ANKRD7, RNFL7, ROPN11, TEX14, PRMT8, OTOA, CFAP99, NRXN1, and DMRT1. PRKACG had the highest degree of nodes. Modules of genes in PPI networks were identified by the MCODE plugin in Cytoscape. The top 3 notable modules were chosen for bioinformatics analysis (Figure 4).

Figure 3

The protein–protein interactome networks and hub genes.

Figure 4

Top 3 modules from the PPI interaction networks. (A) module 1; (B) module 2; (C) module 3; (D) the enriched GO term of module 1 module 2 and module 3.

Survival analysis

In the present study, we observed that RUNX1 mRNA expression was different between RUNX1 mutation and wild-type patients’ bone marrow tissues (Figure 5A). OS analysis indicated that RUNX1 mutant AML patients had a shorter MST than those without RUNX1 mutations (335 days vs. 608 days, log-rank P=0.012). Univariate Cox proportional hazards regression analysis showed that RUNX1 mutations were significantly associated with a poor clinical outcome and an increased risk of death (log-rank P=0.012, HR=2.145, 95%CI=1.156–3.982, Figure 5).

Figure 5

The comparison of mRNA expression and survival between AML patients with RUNX1 mutation and wild type. (A) The mRNA expression of the RUNX1 gene in AML patients’ bone marrow tissue between RUNX1 with mutations and the wild type; (B) Kaplan-Meier survival curves for AML patients stratified by RUNX1 mutation. AML – acute myeloid leukemia; OS – overall survival.

Discussion

It is well accepted that genetic variants can be used as independent prognostic biomarkers for AML due to their potential effect on efficacy of chemotherapy. RUNX1 mutations are the dangerous element in AML. RUNX1 might be a novel biomarker for early diagnosis and a therapeutic target for treatment of AML. Therefore, further investigation is essential for better understanding of the biological roles of RUNX1 mutations in AML. In this study, the KEGG enrichment analysis showed that DEGs were notably abundant in pathways of cancer, ECM-receptor interaction, and cAMP signal. While the cancer pathway was expected and is in line with a previous study [22], it is unclear how the ECM-receptor interaction pathway and the cAMP signaling pathway affect the pathogenesis and prognosis of AML. ECM-receptor interaction pathway enrichment was reported in esophageal squamous cell carcinoma [23]. Li et al. found multiple DEGs and miRNAs were potential biomarkers for the prognosis of epithelial ovarian cancer (EOC) metastasis through the ECM-receptor interaction pathway [24] and affected tumor prognosis [25]. XuS et al. found that prostaglandin E2 (PGE2) receptor EP4 signaling activates the cyclic (c)-AMP signaling pathway, which may be closely related to the occurrence of prostate cancer [26]. It is reported by Kumar N1 that the cAMP signaling pathway was associated with non-small cell lung cancer [27]. Thus, it is possible that RUNX1 mutations may affect development and prognosis of AML through the ECM-receptor interaction pathway and cAMP signaling pathway.

Furthermore, we constructed the protein–protein interactome networks and identified some hub genes associated with RUNX1 mutations. The PRKACG gene, which had the highest degree of nodes, is a member of the protein kinase super-family. PRKACG is known to promote tumor progression because it affects the prognosis of colorectal cancer and leads to poor clinical outcomes [28]. Matos et al. [29] showed that RNF17/TDRD4 (cancer antigen) is a known cancer-associated hub gene of the individual GRNs markers, which are associated with cancer-related diagnostic and prognostic features. RNF17 is a potential liver cancer CT antigen and is associated with unfavorable prognosis [30]. ROPN1L variants are significantly associated with high risk of breast cancer in women [31]. Moreover, Tian et al. reported that ROPN1L was enriched in lung squamous cell carcinoma [32]. Lowe et al. [33] showed that ROPN1L gene was overexpressed in pancreatic cancer. However, there is no correlation between this gene and AML susceptibility and prognosis. Again, our results inspire further research on this gene. Testis-expressed gene 14 (Tex14) is a gene encoding a preferentially expressed protein kinase and promotes the survival of breast cancer [34]. DMRT1 gene is a candidate regulator of sexual development in vertebrates and encodes conserved transcription factor essential for gonadal function. Moreover, DMRT1 promotes testicular germ cell tumors (TGCT) progression [35]. PRMT8 is a type I PRMT that catalyzes the MMA ω-NG-monomethyl arginine and aDMA (aDMA (ω-NG, NG-asymmetric, dimethylarginine) modifications. Yanhong et al. [36] reported that PRMT8 is involved in gliomagenesis molecular pathogenesis. Simandi et al. showed that PRMT8 might also be relevant in the development of human brain malignancy [37]. Otoancorin (OTOA) was reported to be associated with ovarian and pancreatic cancer [38]. NRXN1, also known as KIAA0578, belongs to the neurexin family. Lee et al. found that NRXN1 may result in ALL leukemic cells harboring a novel 3-way translocation t(2;19;11) (p12;p13.3;q23), thereby affecting development of the disease [39].

However, no study has reported the association of CFAP99 or DMRT1 with AML or other cancers. Thus, it is important to further explore the role of CFAP99 and DMRT1 in the pathogenesis and prognosis of AML.

KEGG pathway enrichment showed the top 3 module genes. Genes enrichment of module 1 was prkacg and ankrd7. The genes of module 2 were most abundant in azu1, elane, mpo, ctsg, and rnase3. Furthermore, 3 module genes were abundant in apob and pon1. Apolipoprotein B (ApoB) belongs to the lipoprotein family, which plays important roles in lipid metabolism. In diabetes and metabolic syndrome, serum APOB levels are abnormal [40,41], which can promote the development of cancer [42,43]. Moreover, recent research revealed that high APOB levels can lead to increased risk of lung cancer and colorectal cancer [44]. Serum paraoxonase-1 (PON1) is a 45 kDa glycoprotein [45]. Goncalves et al. [46] found that PON1 rs854560 (L55M) may be related to increased incidence of childhood leukemia, while Çebi et al. [47] reported that PON1 activity can affect the pathology of AML disease.

Statistical analysis indicated that mutated leukemia patients had a shorter MST and poor clinical outcome and increased mortality compared to those without RUNX1 mutation. Our results are consistent with previous reports that RUNX1 mutation adversely affects overall survival (OS) and prognosis of AML patients [48-50]. Germline et al. [51] found that RUNX1 mutations can lead to increased risk of developing AML. The mutations of RUNX1 often result in poor clinical outcomes in AML patients [52]. Gaidzik et al. [49] reported that for AML patients carrying RUNX1-mutation, the rates of event-free survival (EFS), relapse-free survival (RFS), and overall survival (OS) were lower after intensive chemotherapy when compared with wild-type patients. Other studies also showed RUNX1 mutations were related to chemoresistance and increased chance of becoming refractory [53-55]. Tang et al. analyzed RUNX1 mutant AML patients treated with intensive chemotherapy and found that the complete remission rate was lower [54]. The mutations of RUNX1 may lead to increased mortality rate of acute myeloid leukemia. Stengel et al. recently reported that RUNX1-mutant AML showed a distinct pattern of genetic abnormalities and an adverse prognosis [56].

Conclusions

Our results indicate that mutated RUNX1 promotes poor OS in AML patients, which is consistent with previous reports. Results of bioinformatic analysis in the study also reveal that genes such as PRKACG, ANKRD7, RNFL7, ROPN11, TEX14, PRMT8, OTOA, CFAP99, NRXN1, and DMRT1, as well as pathways such as cancer pathways, ECM-receptor interaction pathway, and cAMP signaling pathway, may play a crucial role underlying the effect of RUNX1 mutations in AML prognosis. Our findings provide a novel biomarker for early diagnosis and a therapeutic target for personized treatment of AML. This study had certain limitations: the sample size was small, survival data were incomplete, and survival analysis was only focused on a single factor. Therefore, further studies in molecular pathogenesis and large-scale clinical tumor specimen validation are still needed.