Weighted gene coexpression network analysis identifies critical genes in different subtypes of acute myeloid leukaemia

Abstract

Acute myeloid leukaemia (AML) is a common blood cancer with rapid progression and a high death rate. The aim of this study was to illustrate the molecular mechanism and identify potential prognostic genes by performing weighted gene coexpression network analysis (WGCNA) of AML. WGCNA of AML was performed with R software based on GDC The Cancer Genome Atlas (TCGA) Acute Myeloid Leukaemia (LAML) RNA-seq data from TCGA database. The gene modules showing significant correlation with the M0-M7 subtypes and the NPMc mutation were analysed using Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses. Hub genes from the key modules were identified using the cytoHubba package, and the prognostic values of these hub genes were analysed with a Cox proportional hazards model based on TCGA clinical data. A total of 151 patients from the TCGA database with RNA-seq data were included in the present study. The weighted gene coexpression network contained 21 coexpression modules. Three modules (blue, yellow and purple) were highly correlated with AML subtypes and the NPMc mutation. In total, six key hub genes were identified from the AML subtype- and NPMc mutation-related modules: TLR8, SLC15A3, ADAP2, HOXA6 and HOXA10. Kaplan–Meier curve analysis showed that the five hub genes in the gene expression network are significantly associated with AML prognosis. WGCNA of AML performed in the present study revealed potential prognostic genes that not only have important clinical significance for prognosis but also provide important clues for identifying effective targets in AML therapeutic strategies.

Keywords:

• Acute myeloid leukaemia
• weighted gene coexpression network analysis
• hub gene
• bioinformatics analysis
• RNA-seq

Introduction

Leukaemia, as the most common hematologic disease worldwide, causes more than 300,000 people to die each year [1]. Although the outcomes have improved recently (with a 57% 5-year survival rate in America) [2], leukaemia is still associated with a poor prognosis in developing countries [3]. Generally, genetic [4–6] and environmental factors [7,8] cause leukaemia jointly. Smoking, ionizing radiation, prior chemotherapy, down syndrome and exposure to chemicals such as benzene are major risk factors for leukaemia [9–11]. Leukaemia comprises multiple types (including two acute and two chronic types) [12,13] with a different pathogenesis. Acute myeloid leukaemia (AML) shows rapid progression and causes a high death rate within weeks or months if not treated [14,15]. Records of AML by the American Cancer Society in 2016 showed that 19,950 new cases and 10,430 deaths occurred in the United States [16]. In AML, abnormal cell accumulation in the bone marrow and blood hinders the normal function of blood cells [17,18]. The dysfunction of myeloblasts leads to the uncontrolled growth of immature clones of cells and results in AML [19,20]. In addition, cytogenetic abnormalities are found in AML patients, indicating that these cytogenetic factors may have prognostic significance [21,22]. In addition to cytological changes, mutations in several genes, including CEBPA, NPM1, FLT3, MLL, IDH1 or IDH2, DNMT3A and ASXL1, are associated with the prognosis of AML [23–25].

Weighted gene coexpression network analysis (WGCNA) is a popular approach to build biological gene networks by analysing pairwise correlations between genes [26]. Compared to unweighted correlation networks, WGCNA uses a soft threshold that prevents information loss during processing [27]. To date, this approach has been employed to construct coexpression modules for genes [28], microRNAs [29], long noncoding RNAs (lncRNAs) [30] and even DNA methylations [31]. In addition to the construction of coexpressing modules, WGCNA is also used to find hub genes by determining genes and module eigengenes. By using WGCNA, Pan et al. first identified potential lncRNA modules as prognostic biomarkers in AML [32]. In another study, Luo et al. showed a relapse factor and mechanism in T cell acute leukaemia via WGCNA [33]. WGCNA of single-cell gene expression data revealed a core module between granulocyte/monocyte progenitors (GMPs) and leukaemia networks [34]. These studies opened a new horizon for illustrating complex biological systems in leukaemia pathogenesis.

To elucidate the dynamic changes and identify potential core modules during AML development, we analysed 151 tissue samples from The Cancer Genome Atlas (TCGA) database, representing the stepwise progression of AML from the M0–M7 subtypes. These transcriptome data with clinical trait information also facilitate trait-gene expression correlation studies. To the best of our knowledge, this is the pioneer report elucidating global gene expression changes in AML by WGCNA, although several studies have investigated the gene expression network using unweighted correlation network analyses.

Previously, thymidine kinase 1 [35] and bax [36] expression have been confirmed as prognostic and diagnostic indicator in AML patients. Although these indicators are demonstrated as potential indicators for AML, disease, diagnostic/prognostic indicators are still scarce, especially in different subtypes of AML. In the present study, we aimed to identify new biomarkers for AML early diagnosis and prognosis. Thus, by using a gene expression dataset from databases, we constructed a coexpression network with the WGCNA method and identified 21 modules in these genes. The correlation of AML progression and modules was estimated. The enrichment of canonical GO terms and KEGG pathways of the modules was analysed on three modules closely related to AML subtypes and the NPMc mutation which is a well-known mutation that drives AML development as an oncogene [37]. Then, we performed an upstream regulator analysis and identified hub genes in two modules. The expression changes and possible use of the diagnostic indicators of these hub genes for AML are discussed. The present results highlight the advantages of regulatory networks combined with a functional module analysis in elucidating the mechanism and determining prognostic biomarkers in AML.

Materials and methods

Data and processing

The transcriptome data of AML were downloaded from the TCGA database (https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Acute%20Myeloid%20Leukemia%20(LAML)&removeHub = https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443). A total of 151 samples which had key information of phenotypes were included in the present study (Supporting Information DATA S1 and DATA S2). The patients covered the M0–M7 subtypes, including acute myeloblastic leukaemia with minimal maturation, acute myeloblastic leukaemia with maturation, acute promyelocytic leukaemia, acute myelomonocytic leukaemia, acute myelomonocytic leukaemia with eosinophilia, acute monocytic leukaemia, acute erythroid leukaemia and acute megakaryoblastic leukaemia. The fragments per kilobase of transcript per million mapped reads (FPKM) of the genes were used to construct the network (Supporting Information DATA S1). First, the threshold of the top 25% variance of expression values of the transcripts were retained. Then, we used PCA analysis to check covariates and cofounders of these data.

Construction of the weighted gene coexpression network

WGCNA was applied to identify the coexpression network of all genes in the dataset. WGCNA is a newly developed approach to find key gene sets and hub genes in specific phenotypes of organisms by mathematical methods which is not based on previous biological findings. Using this method, new critical genes participating in AML and potential prognostic and diagnostic indicator for AML patients could be identified. In this study, we used the WGCNA package in R and selected the power of β = 9 as the soft-thresholding estimator to ensure a scale-free network. The parameters were: networkType = unsigned and the minModuleSize = 30. The highly correlated coexpression module is shown as a topological overlap matrix plot. We used a one-step function to generate the network and detect the consensus modules. The constructed modules were assigned different colours. To determine the coexpression similarity of the modules, the eigengene adjacency of the correlations between modules was calculated in a pairwise fashion.

GO and KEGG enrichment analyses of the coexpression modules

The blue, yellow and purple modules that showed a strong correlation with AML subtypes and the NPMc mutation were subjected to GO and KEGG enrichment analyses by DAVID, an online gene functional classification tool (https://david.ncifcrf.gov/). The GO terms and KEGG pathways were confirmed as enrichment when p < 0.05. A heatmap showing the p value of the GO and KEGG enrichment analyses was drawn.

Identification of hub genes

The correlation weight value from each gene in the blue, yellow and tan modules was generated with the WGCNA package. The gene correlations were visualized using Cytoscape (v3.7.1). Using the cytoHubba package, we identified the hub genes that were highly interconnected with nodes in a module with the maximal clique centrality method. The top genes with correlation weight value > 0.1 were defined as hub genes.

Kaplan-Meier survival analysis

Available data and tools from the SurvExpress website (http://bioinformatica.mty.itesm.mx:8080/Biomatec/SurvivaX.jsp) were employed to evaluate the prognostic values for AML. Data were obtained from the TCGA. Two clinical datasets were used: AML TCGA (containing data from 168 patients) and LAML - TCGA - AML (containing data from 149 patients). The hub genes (TLR8, SLC15A3, ADAP2, HOXA6 and HOXA10) were used as input. The risk score obtained from Cox fitting was used to generate risk groups based on the median risk score. Kaplan-Meier plots were drawn, log-rank tests were performed to analyse differences between risk groups, and hazard ratios were estimated.

Results and discussion

Data processing

AML is a common type of leukaemia in adults that threatens human health worldwide. The main objective of the present study was to unveil the molecular mechanism of the subtypes of AML and identify potential candidate gene groups associated with prognosis by constructing a gene coexpression network using WGCNA and data from the TCGA. We identified significantly changed modules that were highly correlated with morphology (M0–M7 subtypes), and the key modules were annotated with GO terms and KEGG pathway enrichment analyses. Subsequently, the hub genes in these modules were identified, and their prognostic values in AML were evaluated.

In total, 151 samples with RNA-seq data were downloaded from the TCGA database. The expression data were used for WGCNA. The clinical characteristics are shown in Table 1. The age of the patients ranged from 20 to 90 years. The 151 patients covered the M0-M7 subtypes, and the NPMc mutation was recorded in 151 patients (33 NPMc-positive and 117 NPMc-negative; Table 1). First, the data downloaded from GEO followed normalization and PCA analysis was performed. The PCA plots showed that the variation in different groups had no significant difference. Thus, WGCNA analysis from the tested samples was reasonable (Figure 1).

Figure 1. PCA analysis of the M0–M7 subtypes using the FPKM of the data for the present study.

Table 1. Clinical characteristics for all patients.

Characteristic (n)	Total
Age
20–30	14
31–40	19
41–50	24
51–60	31
61–70	36
71–80	24
81–90	3
Sex
Male	83
Female	68
Subtypes
M0	15
M1	35
M2	38
M3	15
M4	29
M5	15
M6	2
M7	1
Not Classified	1
NPMc mutant
NPMc Positive	33
NPMc Negative	117
Not Classified	1

Construction of the weighted gene coexpression network from AML samples

The soft-thresholding power based on the criterion of approximate scale-free topology was defined as 9. The value of 9 was the lowest power for the curve that made the scale-free topology fit index flat upon reaching a value above 0.9 (Figure 2(A)) and showed moderate mean connectivity (Figure 2(B)). The weighted gene coexpression network was constructed from the filtered probes, and 21 modules contained 14,536 genes, 511 microRNAs and 74 lincRNAs were identified (Figure 3(A); Supporting Information DATA S3). The interactions of the 21 coexpression modules were analysed with eigengene networks (Figure 3(B)). The transcripts number ranged from 35 to 2920 for each module. Generally, the modules were clustered into two main clades, including a small clade with only 2 modules and a large clade with 19 other modules. The eigengene adjacency heatmap showed the correlations between the modules (Figure 3(C)).

Figure 2. Network topology constructed with the WGCNA package from R to evaluate the soft-thresholding power. The numbers in the figure show the corresponding soft-thresholding powers. The scale-free topology model indicates that the soft-thresholding power of 9 is optimal.

Figure 3. WGCNA of AML. (A) Gene modules identified by WGCNA are shown by clustering the dissimilarity based on the consensus topological overlap using different colours. Each colour represents a module containing a group of highly correlated genes. (B) The dendrogram of consensus module eigengenes based on the consensus correlation. The red line shows the merging threshold. The expression profiles of the module eigengenes below the threshold could be merged due to their similarity. (C) Heatmap of the adjacencies of modules. The adjacency ranged from 0 to 1 (blue to red).

Correlations between the modules and phenotypes

Multiple modules were found to be correlated with the M0–M7 subtypes of AML. In total, six modules had a high correlation (p < 0.05) with the subtypes of AML (morphology code). The highest positively correlated module was the blue module (R² = 0.53), and the highest negatively correlated module was the yellow module (R² = –0.62). In addition to the subtypes of AML, we also observed NPMc, a common mutation in AML. The most correlated module was the tan module, with R² = 0.68 and p = 2e – 21 (Figure 4). These modules were correlated with AML progression and thus focused on the following analyses.

Figure 4. Correlations between consensus module eigengenes and the phenotypes of AML. Rows and columns represent the consensus module and phenotypes, respectively. In the cells, the numbers indicate the correlations of the corresponding module eigengenes and the phenotypes. The numbers in brackets show the p values. The correlation ranged from –1 to 1 (green to red).

Functional enrichment analysis

The highest correlated modules with AML subtypes and the NPMc mutation, including the blue, yellow and purple modules, were used for the functional enrichment analysis. The most enriched GO terms in the purple module were cellular processes, anatomical structure development and the regulation of metabolic processes. In the blue module, the genes were enriched in the immune response, the regulation of immune system processes and regulation of the response to stimuli. The regulation of cellular processes, the regulation of metabolic processes and biosynthetic processes were the most enriched GO terms in the yellow module (Figure 5(A)). Regarding the KEGG pathways, malaria was enriched in the purple module. The most enriched pathways in the blue module were osteoclast differentiation, tuberculosis, lysosome and other pathways. In the yellow module, the most enriched pathways were lysosome, vascular smooth muscle contraction, the cGMP-PKG signalling pathway, and other pathways (Figure 5(B)).

Figure 5. GO and KEGG enrichment analyses of key modules correlated with AML subtypes and the NPMc mutation. (A) The enrichment analysis of GO terms in the tan, blue and yellow modules shown as a heatmap. The colours represent p values ranging from 0.0 to 0.05 (red to white). (B) The enrichment analysis of KEGG pathways in the tan, blue and yellow modules shown as a heatmap. The colours represent p values ranging from 0.0 to 0.05 (red to white).

The most correlated modules with morphology were the blue and yellow modules. The blue module was enriched in immune-related functions, such as the immune response, the regulation of immune system processes and regulation of the response to stimuli, which were changed significantly in all AML subtypes, and the changes in this module were exacerbated with progression. Recently, the tumour microimmune environment has received increased attention with the development of immunology and immunotherapy in cancer [38–40]. Tumour-infiltrating lymphocytes [38] have been demonstrated to participate in several types of cancer, including AML. The immune genes that are also sensitive to chemotherapy and radiotherapy have been regarded as effective prognostic predictors, such as CCL5, CXCL9, CXCL13, CD80, CD21, CD8A and IGKC [41]. Thus, the identification of the immune-related module provides a full glance of the roles of immune genes in the progression of AML.

The yellow and purple modules showed enrichment in functional genes associated with cellular and metabolic processes. It is well known that cellular processes, including proliferation, are altered at the very early stage of AML. When abnormal cell proliferation occurs with metabolic dysfunction in tissues, oncogenesis occurs. It seems that immune changes and cellular and metabolic dysfunction drive AML progression. This evidence suggests that altered energy metabolism is a crucial change in AML, consistent with recent studies [42–44]. Glycolytic flux provides the dominant energy source for cancer cells due to defective mitochondrial respiration [45,46]. We found that in AML, metabolism, especially the regulation of lipolysis in adipocytes, the cGMP-PKG signalling pathway and the purine metabolism pathway, is associated with the development of AML subtypes. These findings suggest that multiple genes constitute modules together that promote AML progression.

Hub gene analysis

To further identify hub genes in the modules related to AML subtypes and the NPMc mutation, we first qualified their relevance using the MMC method based on the eigenvalues of the network modules. We identified the blue and purple modules which were positively correlated with AML subtypes and NPMc. The blue module had the highest gene significance for TLR8, SLC15A3 and ADAP2 (Figure 6(A)). Meanwhile, we found that two microRNAs including miR-373 and miR-4802 were included in these gene networks, suggesting epigenetic regulation by these hub genes. These three hub genes were correlated with 235, 170 and 180 genes, respectively. In the purple module, HOXA6 and HOXA10 were the top two hub genes, which were correlated with 15 and 61 genes, respectively (Figure 6(B); Supporting Information DATA S4).

Figure 6. Protein–protein interaction networks of genes in the blue (A) and purple (B) modules constructed with Cytoscape v3.7.1 software. The red dots represent hub genes identified by cytoHubba, and blue dots represent correlated genes.

Kaplan-Meier survival analysis

Kaplan-Meier survival analysis was used to determine the relationship between the gene expression levels of hub genes and the survival rate of AML patients to show the prognostic values of the genes. We used three genes from the five hub genes, namely, SLC15A3, HOXA6 and HOXA10, as multigene biomarkers for clinical outcomes. The expression levels of the three hub genes were significantly high (p < 0.05) in the high-risk group (Figures 7 and 8). These results indicate that the hub genes are able to effectively distinguish the risk for AML prognosis.

Figure 7. Prognostic analyses of hub genes (TLR8, SLC15A3, ADAP2, HOXA6 and HOXA10) in the AML TCGA dataset based on SurvExpress (n = 168). Kaplan–Meier survival analysis of AML based on the expression of the hub genes. ‘+’ on the right upper figure represents censored samples. The concordance index (CI) and p value of the log-rank test are shown in the figure. Red and green lines indicate the high- and low-risk groups, respectively. The numbers below the x-axis show the death events of individuals over time.

Figure 8. Prognostic analyses of hub genes (TLR8, SLC15A3, ADAP2, HOXA6 and HOXA10) in the LAML - TCGA - AML dataset based on SurvExpress (n = 149). Kaplan–Meier survival analysis of AML based on the expression of the hub genes. ‘+’ in the right upper panel represents censored samples. The CI and p value of the log-rank test are shown in the figure. Red and green lines indicate the high- and low-risk groups, respectively. The numbers below the x-axis show the death events of individuals over time.

In the module, multiple genes that showed a correlation were coregulated in cancer cells. Among these genes, the hub genes that are highly interconnected with nodes in a module have been recognized as functionally significant. The hub genes from the immune-related blue module contained TLR8, SLC15A3 and ADAP2. Thus, TLR8 agonists could be a potential therapeutic strategy for AML. The present results also support the role of TLR8 in promoting AML [47]. High SLC15A3 expression has been shown in a previous study of AML [48]. Moreover, the methylation levels of SLC15A3 were different between LNM-positive and -negative samples, which could be used to predict LNM in ESCC [49]. In the present study, we also demonstrated that these two genes were correlated with miRNAs which may regulate AML in an epigenetic mechanism. In accordance with these results, the involvement of SLC15A3 in AML progression is not surprising. In cancer, ADAP2 binds to β-tubulin and COPRS, which mediates the stability of microtubules and activates oncogenes [50]. However, there are no results regarding ADAP2 function in AML. In the present study, we first found that ADAP2 is a hub gene in an immune-related module that demonstrates the significant functions of these genes in AML. In the tan module, we identified HOXA6 and HOXA10 as hub genes. In addition, we found that the expression levels of HOXA6 and HOXA10 were significantly high in AML. This result is similar to findings in lung adenocarcinomas [51], breast cancer [52] and prostate cancer [53]. However, no such evidence has been reported in AML. These three genes from the Hox family are upregulated in AML. Taken together, these data suggest that the Hox genes are activated with AML progression.

Several studies have suggested that the expression levels of key genes could be considered prognostic indicators. In the novel gene expression network, we evaluated six hub genes as potential prognostic genes, and the results showed that these genes that were highly expressed in AML could be used for prognostic and diagnostic purposes. Cox analysis using the Cox proportional hazards model confirmed these gene sets in two series of AML patients from the TCGA database. The expression levels of FIS1, SPI1, PDCD7 and Ang2 as prognostic indicators have been demonstrated in 100 newly diagnosed cytogenetically normal (CN-AML) patients using the real-time quantitative polymerase chain reaction (qPCR) method [54]. Compared to the qPCR approach, high-throughput sequencing data provide more systematic information. Nevertheless, limitations exist in the present study. We did not validate the levels of gene expression in tissue samples to confirm our findings. Further study is needed to verify our study findings and unravel the mechanism of the hub genes in AML.

Conclusions

The present study provides numerous potential prognostic genes, and the functional modules are presented systematically. A weighted gene coexpression network was constructed of AML with important clinical significance for early diagnosis and prognosis in AML. Of note, TLR8, SLC15A3, ADAP2, HOXA6 and HOXA10 are potential indicators for early diagnosis and prognosis. This pipeline for identifying new biomarkers could be used in other types of disease. Further biological experiments are needed to evaluate the functions and potential prognostic values of these hub genes.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability

The gene expression profiles and phenotype datasets of this study are provided as the Supporting Information Data Sheets S1, S2.

Additional information

Funding

This work was supported by the Science Foundation of Health Commission of Shanxi Province [Grant No. 20162039],  Science Foundation of Education Department of Shanxi Province [Grant No. 2019L686] and Science Foundation of Innovation Team of Changzhi Medical College [Grant No. CX201504].