Upregulation of FHL1, SPNS3, and MPZL2 predicts poor prognosis in pediatric acute myeloid leukemia patients with FLT3-ITD mutation

Abstract

Chromosomal translocations and gene mutations are characteristics of the genomic profile of acute myeloid leukemia (AML). We aim to identify a gene signature associated with poor prognosis in AML patients with FLT3-ITD compared to AML patients with NPM1/CEBPA mutations. RNA-sequencing (RNA-Seq) count data were downloaded from the UCSC Xena browser. Samples were grouped by their mutation status into high and low-risk groups. Differential gene expression (DGE), machine learning (ML) and survival analyses were performed. A total of 471 differentially expressed genes (DEGs) were identified, of which 16 DEGs were used as features for the prediction of mutation status. An accuracy of 92% was obtained from the ML model. FHL1, SPNS3, and MPZL2 were found to be associated with overall survival in FLT3-ITD samples. FLT3-ITD mutation confers an indicative gene expression profile different from NPM1/CEBPA mutation, and the expression of FHL1, SPSN3, and MPZL2 can serve as prognostic indicators of unfavorable disease.

Keywords:

• Acute myeloid leukemia
• FLT3-ITD
• CEBPA
• NPM1
• prognostic indicators

Introduction

Acute myeloid leukemia (AML) makes up about 20% of acute leukemia in pediatric patients [1]. Common clinical symptoms of AML include leukocytosis, anemia, and thrombocytopenia. It is a very heterogeneous disease accounting for more than half of the deaths from leukemia [2,3]. Deciphering disease heterogeneity at the molecular level is crucial for accurate diagnosis, treatment and prognosis, and identifying possible gene therapeutic targets requires deciphering the genetic patterns underlying the etiology of the disease. AML clonal expansion results from abnormal genetic and epigenetic changes in hematopoietic stem and progenitor cells that cause changes or impairment of important physiological processes, such as self-renewal, proliferation, and differentiation [4,5]. Consequently, AML results in the insufficient generation of normal mature blood cells. In addition, AML is associated with multiple chromosomal translocations and mutations that are responsible for the disease pathology and influence disease prognosis. Three gene mutations are proven to be prognostically significant in AML namely, NPM1, CEBPA, and FLT3-ITD [6]. Mutations in FLT3 are associated with a higher rate of relapse and unfavorability [7], while mutations in NPM1 and CEBPA are associated with favorable survival prognoses [8,9].

Of great interest to many research studies are the genome-wide detection of differentially expressed genes (DEGs) from two or more conditions of interest [10]. Currently, one of the main techniques of choice used for gene expression profiling (GEP) is RNA-sequencing (RNA-Seq). GEP can be used to find many systematic differences between cancer and normal conditions thereby defining new clinically relevant disease subtypes [11]. The European LeukemiaNet (ELN-2017) provides guidelines to stratify AML patients into good, intermediate, and adverse risk groups based on cytogenetics and mutation status of some genes, including ASXL1, CEBPA, FLT3, NPM1, RUNX1, and TP53 [12].

Based on results of a preliminary data analysis of high- and low-risk AML samples from the GDC TARGET AML dataset, the PCA analysis showed that samples with NPM1, CEBPA, and FLT3-ITD mutations were clustered together (Figure 1), although they have different survival rates. In addition, samples with t(8;21) and inv(16) variations were well separated into two other independent clusters (Figure 1). In line with these results, a 36-gene expression signature enabled the accurate classification of AML samples with t(8;21) and inv(16) mutations [13]. However, segregation between AML samples with NPM1 and CEBPA using GEP is less accurate, revealing that the expression pattern between the two types of samples is likely similar [14]. Therefore, the fact that FLT3-ITD samples are grouped with NPM1 and CEBPA samples of favorable prognosis raises the question of why the survival of AML patients with FLT3-ITD is low? Thus, it is crucial to look for DEGs between AML FLT3-ITD on the one hand and AML NPM1/CEBPA on the other hand, which are potentially the cause of poor survival in FLT3-ITD samples.

Figure 1. Clustering of the high- and low-risk AML samples based on their cytogenetics and genetic aberrations using principal component analysis (PCA).

This study aims to identify a gene expression signature that can differentiate AML FLT3-ITD mutation from the AML NPM1 and AML CEBPA mutations, which are implicated in poor survival. In addition, we examined the predictive value of this signature in risk classification. These novel gene expression signatures may assist clinicians and pathologists in patient diagnosis and assessment, thereby ensuring more accurate and individualized treatment options.

Materials and methods

The workflow describing the steps and methods undertaken in this study is illustrated in Figure 2. It includes six essential steps: dataset retrieval, differential gene expression (DGE), PCA and clustering, machine learning (ML), Cox regression, and Kaplan–Meier’s (K–M) analyses.

Figure 2. Workflow depicting the steps and methods used for the identification of genetic signature associated with poor prognosis in AML samples with FLT3-ITD mutation. AML: acute myeloid leukemia; CEBPA: CCAAT enhancer binding protein alpha; DGE: differential gene expression; FLT3-ITD: FMS-like tyrosine kinase 3-internal tandem duplications; MLP: multi-layer perceptron; NPM1: nucleophosmin-1; PCA: principal component analysis; TARGET: therapeutically applicable research to generate effective treatment.

Datasets

The therapeutically applicable research to generate effective treatment (TARGET) project employed a multi-omic strategy to molecularly characterize hard-to-treat pediatric cancers, including AML. TARGET data are accessible through the TARGET data matrix and the UCSC Xena Browser, a web-based visual integration and exploration tool for multi-omic data. The UCSC Xena Browser is a high-performance visualization, exploration, and analysis tool for multi-omic data of large public repositories and private datasets [15].

The TARGET AML dataset in Xena consists of a total of 187 samples with their accompanying clinical data. XenaPython, a Python package implementing APIs to query and download data from Xena, was used to obtain gene expression in log2 transformed RNA-Seq counts format of the TARGET AML dataset (dataset ID: TARGET-AML.htseq_counts.tsv) from the GDC hub. The criteria used to query the dataset were the risk group (low and high) into which the samples were classified. The data samples were then filtered based on their cytogenetic abnormalities and genetic mutations. AML data samples with FLT3-ITD, NPM1, and CEBPA mutations were used for downstream analysis. Samples with multiple gene mutations were not considered, as multiple mutations may impact AML prognosis [16]. There is no precise risk classification of patients with FLT3-ITD and NPM1 double mutation. The National Comprehensive Cancer Network considers them to have a favorable or intermediate risk, however, an unfavorable prognosis was manifested by this category of patients in recent studies [17]. Similarly, CEBPA-mutated patients with WT1 mutations were reported to have an unfavorable risk [18]. Due to these uncertainties and changes in the risk classification of AML patients with multiple mutations, two samples, one with FLT3-ITD and NPM1 double mutation and the other with both CEBPA and WT1 mutations, were not considered.

Differential gene expression analysis

The gene expression values from the AML samples with FLT3-ITD, NPM1, and CEBPA were converted to raw counts using the mathematical expression R=$2^L$ – 1, where R is the raw count and L is the log2 normalized value. Filtering to remove low expressed genes was done using edgeR’s filterByExpr function [19]. Statistically, eliminating low expressed genes enables the mean–variance association in the data to be measured more accurately and reduces the number of computational checks conducted in downstream differentially expressed tests [10]. The DGE analysis between FLT3-ITD and NPM1/CEBPA samples was performed using the DESeq2 package in R [20]. DESeq2 uses shrinkage estimators for dispersion and fold change for comparative DGE estimation [20]. Genes that met the criteria of an adjusted p value <.01 were considered significant and therefore differentially expressed.

Machine learning

To determine if the DEGs identified could serve as risk classification biomarkers in AML, a multi-layer perceptron (MLP) classifier was constructed using the pediatric DEGs and tested on adult AML samples, as explained below. MLP is a feed-forward artificial neural network (ANN) algorithm. It is implemented in Python under the Scikit-learn library [21] and imported using the MLPClassifer class. Prior to the classification task, a linear support vector machine and repeated stratified 10-fold cross-validation were used for recursive feature elimination (RFE). The RFE was applied to the DEGs to obtain the essential features for optimum model performance. The selected genes were then used as features for classification. Also, principal component analysis (PCA) and clustering analysis were applied to the expression values of these selected genes to assess their ability to discriminate the samples of the two AML groups. The PCA plot and clustering heatmap were generated using the ggplot2 and pheatmap R packages, respectively.

The features of the training set were extracted from the log2-transformed normalized counts of the pediatric AML samples. The test set was constructed from the adult AML dataset (dataset ID: TCGA-LAML.htseq_counts.tsv) obtained from the GDC hub in the UCSC Xena Browser. Construction of the test set followed the same process used to construct the training set, i.e. consider AML samples with FLT3-ITD mutation and AML samples with NPM1/CEBPA mutation. The features (genes) not present in the test dataset were removed from the training set and were not used in the ML classification. In addition, the gene expression values of the training and test sets were scaled using Scikit-learn’s StandardScaler function. The MLP model was then applied to predict the classification of the samples in the test set. The evaluation metric for the model performance was the accuracy of classification.

Cox’s regression analysis and Kaplan–Meier’s estimates

To determine among the genes from ML those that best correlated with patient's survival, we used a Cox regression model based on the Lasso algorithm of the glmnet R package [22–24]. The model assigns each gene a regression coefficient value. Genes with a zero coefficient were eliminated, having no effect on survival. Prognostic genes with positive coefficients suggest that their upregulation signifies low survival in FLT3-ITD patients. For each gene with positive coefficient, a score value is calculated for each patient as a product of the expression value of the gene and its corresponding coefficient obtained from the Cox regression model. A median value was inferred from the patient scores. Each score was then compared to the median, and patients were assigned a status value of 1 or 0 depending on whether the score was above or below the median. According to patient status information, K–M estimates were then calculated for overall survival (OS). K–M curves were generated using the ggsurvplot function from the survminer R package.

Results

Differential gene expression analysis

Querying UCSC Xena Browser for TARGET AML samples using risk group as search criteria returned 12 high-risk samples and 72 low-risk samples. The distribution and clustering of samples by cytogenetic and mutational aberrations are shown in Table 1 and Figure 1, respectively. Samples with FLT3-ITD, NPM1, and CEBPA mutations were selected for downstream analysis, comprising of 10 samples of FLT3-ITD mutation and 13 NPM1/CEBPA samples. All were clustered together in the exploratory analysis (Figure 1). Filtering of the gene expression data eliminated 38,980 low expressed genes and retained 21,503 genes for downstream analysis. The DGE analysis using DESeq2 identified 471 DEGs between the AML FLT3-ITD and AML NPM1/CEBPA groups. Of these, 208 genes were upregulated, while 263 genes were downregulated.

Table 1. Type of mutation, number of samples and risk category of TARGET high- and low-risk samples.

Type of mutation	Number of samples	Risk classification
CEBPA	9	Favorable
NPM1	6	Favorable
Del(5q)	1	Unfavorable
FLT3-ITD	11	Unfavorable
inv(16)	32	Favorable
t(8; 21)	25	Favorable

Machine learning

The RFE selected 22 genes out of 471 DEGs as the most critical features for classifying samples according to their mutational status (Table 2). These 22 DEGs accurately discriminated and clustered the TARGET samples into their respective mutation groups (Figure 3). However, some of these selected genes were not found in the adult test dataset and were not included in the pediatric training set used to create the MLP model. Sixteen DEGs (features) were thus used for the ML classification, including FHL1, VWF, MPZL2, SPNS3, LINC00515, TCEA3, DCN, MACC1, ADD2, CCDC152, KCNA6, ADAMTS3, KCTD15, KCNMB4, HESX1, and IL3RA. The training set was constructed based on the gene expression values of the TARGET samples from the FLT3-ITD group with 10 samples and the NPM1/CEBPA group with 13 samples. The test set was 25 samples, including five FLT3-ITD samples and 20 NPM1/CEBPA samples. Evaluation of the MLP classifier on the test set yielded 92% accuracy of samples correctly assigned to their corresponding mutational status. All FLT3-ITD samples were correctly predicted, while 18 of the 20 NPM1 and CEBPA samples were correctly predicted.

Figure 3. Clustering of the AML samples with FLT3-ITD and NPM1/CEBPA mutations using the 22 DEGs selected by RFE using (A) principal component analysis (PCA) and (B) hierarchical clustering using pheatmap R package.

Table 2. List of the 22 DEGs selected by RFE genes with * was used in the ML model.

Ensmbl_Id	Symbol	log2FoldChange	padj
ENSG00000250696.4	LOC105377267	–7.31571544058504	5.57786607011193E-16
ENSG00000022267.15	FHL1*	3.3904649981173	1.74302159426198E-11
ENSG00000110799.12	VWF*	5.57344738786003	7.63061866000091E-11
ENSG00000182557.6	SPNS3*	2.34741592806669	1.73370534268721E-08
ENSG00000149573.7	MPZL2*	4.11373735379172	2.24605771346956E-07
ENSG00000182183.13	SHISAL2A	–2.80284553156498	3.31597387364724E-07
ENSG00000011465.15	DCN*	–4.93071112156197	3.74939407438753E-07
ENSG00000183742.11	MACC1*	3.11711209042932	4.62361610488009E-06
ENSG00000156140.7	ADAMTS3*	–6.20242014449732	5.12366228681635E-06
ENSG00000075340.21	ADD2*	–3.56747714605576	6.92781902701875E-06
ENSG00000135643.4	KCNMB4*	–2.96743485475082	1.20611585085678E-05
ENSG00000153885.13	KCTD15*	–2.73614037690847	2.01801229683281E-05
ENSG00000247774.5	PCED1B-AS1	–1.90010539713958	0.000355270349374
ENSG00000224420.3	ADM5	–1.99554499874056	0.000408776507021
ENSG00000204219.8	TCEA3*	–1.99082117539315	0.000651667363783
ENSG00000163666.7	HESX1*	1.59282492017001	0.001841360690536
ENSG00000185291.9	IL3RA*	1.36386204250051	0.002217439812688
ENSG00000178075.18	GRAMD1C	1.93515803839245	0.003109807383799
ENSG00000130035.5	KCNA6*	–2.40840491942698	0.003146881886037
ENSG00000260583.1	LINC00515*	2.48087500485109	0.003149939931692
ENSG00000198865.8	CCDC152*	2.70119814011179	0.003498603901287
ENSG00000138678.9	GPAT3	2.22603438762192	0.004203371885436

Cox’s regression analysis and Kaplan–Meier’s estimates

The Cox regression model associated five of the 22 DEGs obtained from RFE with patient survival in the two AML groups. The model created using the LASSO algorithm assigned non-zero, positive or negative, coefficients to the five genes (Table 3). The genes with a positive coefficient (FHL1, SPNS3, and MPZL2) had high expression in the FLT3-ITD group with over fivefold increase, while those with a negative coefficient (KCNMB4 and ADD2) had high expression in the NPM1/CEBPA group with over sevenfold increase. Kaplan–Meier’s estimates for OS based on patient statuses of each gene with a positive coefficient were derived and presented in Figure 4. Expression of the three genes in the FLT3-ITD group correlated with low survival (p<.0001). Patients with a score above the median were predominantly FLT3-ITD, and those below mainly were NPM1/CEBPA.

Figure 4. Kaplan–Meier’s plots of overall survival of patients with FLT3-ITD and NPM1/CEBPA mutations for survival scores above the median (continuous line), corresponding to high expression, and below the median (dotted line), corresponding to low expression, for the genes with positive Cox coefficients (A) FHL1, (B) SPNS3, and (C) MPZL2.

Table 3. The list of the DEGs related to overall survival and their corresponding coefficient values from the Cox regression model.

Gene name	Coefficient
FHL1	0.077996384008932
SPNS3	0.050799306065861
MPZL2	0.08555275154584
ADD2	–0.000713971400328
KCNMB4	–0.03112122883124

Discussion

Proper risk stratification of AML patients at diagnosis is essential for making optimal therapeutic decisions. The results of the DGE analysis demonstrated that the gene expression profile of FLT3-ITD AML is different from that of the other mutation groups mentioned above. FLT3-ITD mutation is associated with unfavorable AML disease outcome [7]. Similar prognostic driver genes have been cataloged in several cancers. For example, in neuroblastoma, another pediatric cancer, MYCN amplification is a predictor of poor prognosis [25]. Some of the DEGs identified here have been associated with AML and several other cancers [26–28]. The downstream analysis identified DEGs associated with survival (Table 3) and possible responsibility for the poor survival in AML FLT3-ITD as opposed to AML NPM1/CEBPA, which clustered together in the exploratory analysis (Figure 1).

DEGs in the expression signature

Some of the genes in the identified signature have been reported to have tumor-suppressive roles and downregulated in several cancers. TCEA3, a gene with apoptosis-promoting functions [29], was downregulated in the FLT3-ITD group in our study. It has been reported downregulated as well in rhabdomyosarcoma [30], gastric cancer [31], and ovarian cancer [32]. DCN, another downregulated gene in the FLT3-ITD group, is well known for its oncosuppressive function [33]. DCN expression is decreased in different cancer types [34–36]. In contrast, some upregulated genes have tumor-promoting roles. VWF, upregulated in the FLT3-ITD group, is a glycoprotein involved in hemostasis whose expression has been associated with lymph node metastasis in prostate cancer patients [37]. It was also demonstrated to be related to metastatic activities in glioma and osteosarcoma cells [38].

MACC1 regulates the HGF-MET pathway, which is involved in cellular activities like growth, motility, angiogenesis, invasiveness, epithelial–mesenchymal transition, and metastasis [39,40]. As expected, its high expression has therefore been reported in some cancers such as retinoblastoma [41], gastric cancer [42], and colon cancer [28]. MACC1 is an independent prognostic indicator for metastasis in colon cancer [28]. LINC00515 has also been reported upregulated in glioma [43] and multiple myeloma [44]. IL3RA expression was also correlated with FLT3-ITD-mutated AML, although it is also expressed in some NPM1-mutated AML cases [45]. High expression of IL3RA also seemed to be associated with worse clinical outcomes in AML [46]. Arai et al. [47] reported IL3RA expression associated with chemotherapy response failure and poor survival in de novo AML patients. It, therefore, suggests that FLT3-ITD mutation causes the downregulation of protective tumor-suppressing genes while upregulating genes with oncogenic function.

Machine learning

The high classification accuracy of 92% obtained indicates that the 16-gene expression signature can be used as prognostic markers of FLT3-ITD mutation and unfavorable disease outcome in AML and could be potential therapeutic targets. This high accuracy was achieved despite using an adult AML dataset as the testing set. Differences in molecular mechanisms between pediatric and adult AML have been reported [3,48,49]. Our classification results suggest that pediatric and adult AML have likely similar disease mechanisms for the mutations investigated in our study. Different gene expression signatures for prognosis in AML have been proposed before. Using expression microarray, Bullinger et al. [50] found a 20-gene expression signature predictive of FLT3-ITD mutation status using prediction analysis of microarray method [51] and Kaplan–Meier’s analysis. They obtained a prediction accuracy of 81% for identifying cases with FLT3-ITD mutation. Our study used an RNA-Seq dataset for analysis, performed DGE analysis, and applied machine-learning methods to find and validate a 16-gene expression signature predictive of FLT3-ITD mutation with 92% higher accuracy achieved. RNA-Seq has become the primary technology used for GEP [10]. Our expression signature of a lesser number of genes should be implementable and potentially cost-effective.

Similarly, Zhu et al. [52] found a six immune-related gene signature to predict prognosis in AML using Cox and LASSO regression analysis. However, their classifier is based on immune-related genes, and the genes identified may not provide a larger picture of the mechanisms at play in the disease pathogenesis. As a limitation, the sample size of our test dataset was small. Therefore, validation in much larger sample cohorts may still be essential. However, we expect that further studies of the genes identified in our study will shed light on the underlying pathogenesis of AML with FLT3-ITD mutation.

Cox’s regression analysis and Kaplan–Meier’s estimates

By reducing the list of genes obtained from the ML analysis, the Cox regression model selected five genes related to patient OS including three genes with positive Cox coefficients. The K–M curves (Figure 4) show the ability of these three genes to distinguish between shorter and longer OS, i.e. FLT3-ITD and NPM1/CEBPA groups, respectively (p<.0001). The results of the survival analysis (Figure 4) validated the poor prognosis in the FLT3-ITD group. Since a patient score above the median is associated with low OS, higher expression of FHL1, SPNS3, and MPZL2 with positive Cox coefficients is a prognostic indicator of poor outcome in AML with FLT3-ITD mutation. The roles of the genes with negative Cox coefficients, ADD2 and KCNMB4, in the survival of AML patients with NPM1/CEBPA require further investigation.

FHL1 belongs to a family of genes that play a role in focal adhesion and differentiation. High expression of FHL1 has been described as a powerful prognostic indicator of worse survival and poor outcome, independent of existing genetic factors for the prognosis of AML. Its overexpression was further found to be related to chemotherapy-resistance and relapse in AML, and FHL1 knockdown enhanced the sensitivity of AML cells response to treatment [53]. SPNS3 is a transmembrane transporter whose overexpression has been associated with poor prognosis in AML patients receiving chemotherapy or allogeneic hematopoietic stem cell transplantation [54]. Huang et al. [54] suggested that over-expression of SPNS3 may regulate and control proliferation and differentiation of AML by autophagy. MPZL2 expressed in the lymphoid organ, thymus, and other epithelial structures maintains stemness in glioblastoma [55]. Its increased expression has been reported in AML and hepatocellular carcinoma, which is associated with poor prognosis and recurrence [56,57].

Notably, the high expression of these genes in FLT3-ITD patients promotes poor outcomes, and the implication of FHL1 and SPNS3 in chemotherapy resistance could be the reason for treatment failure, which may therefore explain the low survival in this category of patients. Functional studies such as knockdown experiments to decipher the precise roles of these genes in FLT3-ITD AML would be illuminative.

Conclusions

In this study, we differentiated AML FLT3-ITD and AML NPM1/CEBPA using GEP. We also identified and validated a 16-gene expression signature for risk classification in AML, which has diagnostic and prognostic value. The upregulation of FHL1, SPNS3, and MPZL2 was found to be associated with poor survival in AML FLT3-ITD. These genes could therefore be potential therapeutic targets in pediatric and adult AML.

Disclosure statement

Authors have no conflict of interest to declare.

Additional information

Funding

The work reported herein was made possible through funding by the South African Medical Research Council (SAMRC) through its Division of Research Capacity Development under the Mid-Career Scientist Programme from the South African National Treasury. The South African Research Chairs Initiative of the Department of Science and Innovation and the National Research Foundation of South Africa (Grant ID 64751) also supported this work; the Poliomyelitis Research Foundation; South African National Research Foundation.