Prognostic significance of ferroptosis-related genes and their methylation in AML

ABSTRACT

Background

Ferroptosis involves in the development and therapeutic response of various types of tumors. This study aims to explore ferroptosis-related prognostic genes that could further accurately stratify AML patients.

Methods

We investigated the prognosis significance of ferroptosis-related genes in AML by Univariate and multivariate Cox proportional hazards regression analyses. With the methylation data of TCGA samples, we looked for methylation sites associated with prognostic genes and compared the correlation between methylation and mRNA expression. R software and ‘edgeR’ packages were used to identify the DEGs between the high-and-low-risk groups divided by the FRPGs prognosis model and then run GO enrichment, KEGG pathway, and PPI network.

Results

We found a prognostic risk model that included AKR1C2 and SOCS1 predicted outcomes in AML patients. Methylation analysis showed that AKR1C2 and SOCS1 are negatively regulated by their methylation, leading to their low expression in AML patients. Besides, both decreased SOCS1 expression and hypermethylation predicted favorable OS and PFS in AML patients. Finally, this prognostic risk model exhibited a close correlation with several clinical features, especially with age (P=0.005), cytogenetic type (P=0.031), risk_cytogenetic (P=0.001), and risk_molecular (P<0.001). Functional enrichment analysis showed that DEGs are most enriched in the regulation of cell death and the PI3K-Akt signaling pathway.

Conclusion

AKR1C2 and SOCS1 are promising biomarkers for predicting prognosis in patients with AML.

KEYWORDS:

• Ferroptosis
• AML
• AKR1C2
• SOCS1
• methylation

Introduction

Acute myeloid leukemia (AML) is a group of molecularly heterogeneous hematologic malignancies. Its initiation mainly attributes to changes in pivotal cellular regulatory factors, such as transcription and epigenetic factors [1]. Currently, cytogenetic outcomes and molecular abnormalities at diagnosis consider being the most important prognostic factors for AML and highly predictive of complete response rate(CR), disease-free survival(DFS), overall survival(OS), and recurrence risk [2,3].

Current clinical guidelines classify AML into three groups (favorable, intermediate, and adverse risk) according to cytogenetic and molecular risks of AML patients. The intermediate-risk group is the most subgroups, accounting for approximately 45% of adult AML patients. AML patients with normal karyotypes (NK-AML) and patients without adverse or favorable molecular abnormalities are in this subgroup [2]. The ideal therapy strategy for these patients endures unclear, and treatment outcomes are heterogeneous. Several studies have found molecular alterations to further refine risk stratification (such as KIT, AF1q, MN1, ERG, MLL, and BAALC) in this subgroup of patients [4–6]. It is necessary to further perfect the risk stratification of AML patients and provide more accurate and appropriate therapy for the patients with intermediate-risk.

Epigenetic modification abnormalities appear in a large proportion of AML patients [7]. As a critical regulator of gene expression, DNA methylation study may also help to the more precise grouping of AML patients. DNA methylation is a key epigenetic mechanism involved in transcriptional regulation and cell development and function [8]. Hypomethylation rules in introns and repeats elements, while hypermethylation usually happens in so-called CpG islands and is correlated with gene promoter regions. Numerous studies have established a link between abnormal DNA methylation and gene silencing in disorders [9–11]. These suggest the potential importance of DNA methylation as a prognostic and risk marker in AML.

Ferroptosis is an iron-dependent form of regulatory cell death caused by excessive lipid peroxidation, associates with the development and therapeutic response of various types of tumors [12,13]. In this study, we first selected the ferroptosis-related genes (FRGs) of AML in the TCGA database, explored the ferroptosis-related prognostic genes (FRPGs) and their DNA methylation in AML. Then, we validated the prognostic-related genes in patients at intermediate risk classified by the NCCN AML classification according to their molecular and cytogenetic risk status. Finally, according to the prognostic model, patients were divided into low- and high-risk groups. We calculated the differentially expressed genes (DEGs) between the two groups and examined the biological processes of DEGs through gene enrichment analysis to study the functional mechanism of FRPGs in AML.

Materials and methods

Acquisition of ferroptosis-related genes in AML

First, we downloaded the clinical data, transcriptional data, along with methylation profiles of patients with AML from the TCGA database via the Cbioportal website (https://www.cbioportal.org/study/summary?id=laml_tcga_pub). The inclusion criteria were patients with complete transcriptional, methylation, and clinical data. 154 AML samples were retrieved for subsequent analysis. Then, we downloaded the FRGs data from the FerrDb website (http://www.zhounan.org/ferrdb/operations/download.html). FerrDb is the database of ferroptosis regulators and markers and ferroptosis-disease associations [14]. The genes in AML of the TCGA database were intersected with FRGs to obtain FRGs in AML.

Construction of a prognostic signature

We performed the univariate Cox regression analysis of FRGs in AML by using the ‘survival’ package in R software to obtain the FRPGs (P<0.05). Then, we searched an online website for gene expression profiling interactive analysis (GEPIA) to investigate the DEGs of the FRPGs in AML samples and healthy control samples. Moreover, the DEGs of prognostic significance were further put into multivariate Cox regression analysis to obtain the FRPGs for constructing the prognosis model and generated the regression coefficients of these crucial genes. To reflect the effect of the prognostic model on the survival status of AML patients, we calculated the risk score of each patient by multiplying the normalized gene expression of each FRPG with its corresponding regression coefficient and divided patients into low-and- high-risk groups according to the mean risk score. We further conducted the Kaplan-Meier curve analysis to estimate the relationship between the risk score, overall survival, and disease-free survival. We also performed time-dependent receiver operating characteristic (ROC) curve analysis to evaluate the 3-year and 5-year predictive value of the os status and calculated the area under the curve (AUC) of the ROC curve to examine the classifier performance. Besides, we analyzed the relationship between the different clinical parameters and the risk score based on FRPGs. To explore whether these prognostic genes can use as supplements to the current AML molecular risk stratification, we further applied the prognosis model in patients at intermediate risk classified by the NCCN AML classification according to their molecular and cytogenetic risk status.

DNA methylation analysis of prognostic genes

With the methylation data of TCGA samples, we looked for methylation sites associated with prognostic genes and compared the correlation between methylation and mRNA expression. MethSurv (https://biit.cs.ut.ee/methsurv/) contains the clinical and methylation data on AML in the TCGA dataset and provides online analysis [15], was used to draw the heat maps of methylation sites and to do survival analysis of methylation sites. We also investigated the detailed association of FRPGs and their methylation with a range of clinical features. Graph of methylation and gene expression, survival curves, and density figures of the methylation sites were plotted.

GO and KEGG functional enrichment analysis

In the part of construction of a prognostic signature, the TCGA cohorts were divided into high-and-low-risk groups by multivariate regression analysis. Then, the ‘edgeR’ package of R software was used to compare the gene expression profiles of patients in high-and-low-risk groups and the DEGs between these two groups were obtained. Then, analyzed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on the DEGs (P-value < 0.05, |log2FC| ≥ 1.5) between the high-risk and low-risk groups using the ‘clusterProfiler’ R package.

Protein–protein interaction network construction and module analysis

The protein–protein interaction (PPI) network of DEGs was constructed using the gene/protein retrieval tool (https://string-db.org/), and the results were visualized by Cytoscape software (version 3.8.0). Cytoscape MCODE plug-in provided access to select the significant modular of PPI network [16,17]. The patterns default parameters are as follows: degree cutoff=3, node score cutoff=0.2, k-core=2, and max. depth=100. For genes in the hub modules, we use ClueGO + CluePedia for functional enrichment analysis.

Statistical analysis

The OS and DFS between different groupswere compared by Kaplan-Meier analysis with the Nonparametric test. The student's t-testwas used to compare gene expression between high-and-low groups. Pearson correlation and Spearman's rho wereused to analyzing the correlation between methylation and mRNA expression. A chi-square test wasused to investigate the detailed association of prognostic-related genes, AKR1C2 expression, SOCS1 expression, and their methylation with a set of clinical features. All statistical analyses were performed with R software and SPSS software, and statistical significance was set at probability values of p < 0.05.

Results

Ferroptosis-related gene extraction

In this study, we used three FRGs data sets provided by the FerrDb database, including driver genes that promote ferroptosis, suppressor genes that prevent ferroptosis, and marker genes that indicate the occurrence of ferroptosis. Among them, four genes were listed as drivers and inhibitors simultaneously and were excluded from our study. Excluding the multi-annotated genes, 130 overlapping ferroptosis-related genes from the TCGA dataset and FRGs were selected for further analysis. The Wien diagram isshown in Figure 1.

Figure 1. Wien diagram of FRGs and genes of AML in TCGA dataset

Construction of a prognostic signature

Populations of 154 patients from the TCGA dataset were identified and included in this study. The univariate Cox proportional hazards regression analysis showed that 20 genes had significant prognostic relevance (p < 0.05) (Table 1). GEPIA identified 6 genes that were differentially expressed between AML samples and normal samples (Figure 2A-2F). Multivariate Cox proportional hazards regression analysis was further performed on the 6 genes, which screened AKR1C2 and SOCS1. The risk score for predicting overall survival was calculated as follows: Risk score=0.160* AKR1C2 + 0.383* SOCS1 (Figure 3A). According to the risk score, patients were divided into low- and high-risk groups. Survival analysis showed that low-risk patients had longer OS and DFS than high-risk patients (median OS: 33.56 months vs. 14.04 months, p<0.0001; median DFS:26.79 months vs. 9.55 months, p<0.0001) (Figure 3B、3C). The AUC of 3-year and 5-year os survival time roc curve analysis were 0.73 and 0.76, respectively (Figure 3D).

Figure 2. The expression of the FRPGs between AML patients and their normal samples. A. AIFM2 B. AKR1C2 C. ENPP2 D. LPCAT3 E. MT1G. F. SOCS1. Red represents AML samples; Blue represents healthy control samples.

Figure 3. Survival prognosis model on the 2 FRPGs. A. Forest map of multivariate cox regression. B. Survival analysis showed that the patients in the high-risk group had worse overall survival than those in low-risk group in TCGA cohort. C. DFS analysis showed that the patients in the high-risk group had shorter DFS than those in low-risk group in TCGA cohort D. ROC analysis was performed to calculate the most optimal cutoff value to divide the AML patients into high risk and low risk group.

Table 1. Significant prognostic genes obtained by univariate Cox proportional hazards regression analysis.

Gene Symbol	HR	z	pvalue
SOCS1	1.476447	4.991562	5.99E−07
MYB	0.671575	−3.22989	0.001238
G6PD	1.45728	2.948866	0.003189
AKR1C2	1.172824	2.941784	0.003263
LAMP2	0.463622	−2.9167	0.003538
PIK3CA	0.625872	−2.62536	0.008656
MT1G	1.126474	2.619178	0.008814
SRC	1.199648	2.618093	0.008842
NCOA4	0.624779	−2.59129	0.009562
OTUB1	1.535339	2.393224	0.016701
MAPK8	0.663913	−2.3275	0.019939
YY1AP1	2.053237	2.303706	0.021239
HRAS	1.424803	2.291709	0.021922
ACVR1B	0.766157	−2.18144	0.029151
AIFM2	1.20518	2.14205	0.032189
BACH1	0.688397	−2.1401	0.032346
HSF1	1.494968	2.130842	0.033102
TMBIM4	0.663498	−2.10892	0.034951
LPCAT3	1.505839	2.0978	0.035923
ENPP2	1.101019	2.082965	0.037254

We then used age as a control variable to compare the risk score based on the 2 FRPGs with different clinical factors, such as BM blast percentage, PB blast percentage, cytogenetic risk, molecular risk, WBC count, and sex on survival. The results showed that risk scores based on 2 genes, cytogenetic stratification, and molecular stratification were more common than age factors and increased the occurrence of survival. The occurrence of BM blast percentage, PB blast percentage, WBC count, and sex were lower than that of age and did not affect the occurrence of survival (Figure 4).

Figure 4. Cox Forest of the 2 FRPGs with different clinical features.

Validation of the prognostic value of the 2 FRPGs in AML patients with intermediate risk

The prognostic model was further applied in patients at intermediate risk classified by the NCCN AML classification according to their molecular and cytogenetic risk status. We also found that among these intermediate-risk patients, those low-risk patients stratified according to FRPGs had longer overall survival than high-risk patients (median OS of patients stratified according to cytogenetic risk: 33.87 vs. 16.07 months, p<0.0007; median OS of patients stratified according to molecular risk: 33.73 vs. 16.69 months, p<0.0014, Figure 5A, 5B). A similar and statistically significant trend of DFS was also observed (median DFS of patients stratified according to cytogenetic risk: 26.53 vs. 10.64 months, p<0.0004; median DFS of patients stratified according to molecular risk: 25.94 vs. 10.73 months, p<0.0001, Figure 5C,5D).

Figure 5. Survival prognosis model on the 2 FRPGs in patients with intermediate risk. A. OS analysis of patients stratified according to cytogenetic risk. B. OS analysis of patients stratified according to molecular risk C. DFS analysis of patients stratified according to cytogenetic risk. D. DFS analysis of patients stratified according to molecular risk.

DNA methylation analysis of FRPGs

We explored the methylation status of AKR1C2 and SOCS1 genes and found that AKR1C2 has 4 methylation sites and SOCS1 has 24 methylation sites. Figure 6A,6B clearly shows the distribution of CpG sites of each gene. Then Pearson correlation and Spearman's rho were used to analyzing the correlation between methylation and mRNA expression. We observed that two methylation sites were negatively correlated with gene expression. The relationship between methylation and mRNA expression is shown in Figure 7A,7B. We then selected cg18841653 and cg03014241 CpG sites to investigate the prognostic values of these significant DNA CpG sites in patients with AML. Kaplan-Meier plots demonstrated that a high level of methylation of the cg03014241 CpG site was correlated with more favorable OS (Figure 8A) among patients with AML (median OS: 25.24 vs. 22.93 months, P=0.001). Then, the AML patients in TCGA were dichotomized into low or high subgroups when divided by the median value of SOCS1 mRNA expression. In contrast to the effect of methylation on prognosis, survival analysis showed that a high level of mRNA expression was correlated with worse OS (Figure 8B) among patients with AML (median OS: 22.93 vs. 25.24 months, p<0.041).

Figure 6. The heatmap of CpG sites of AKR1C2 and SOCS1 A. The distribution of CpG sites of AKR1C2. B. The distribution of CpG sites of SOCS1.

Figure 7. A. The relationship between methylation and AKR1C2 mRNA expression. B.The relationship between methylation and SOCS1 mRNA expression.

Figure 8. A. Kaplan-Meier plots of the cg03014241 CpG site. B. Kaplan-Meier plots of the SOCS1 expression.

We used the chi-square test to investigate the detailed association of FRPGs, AKR1C2, SOCS1 expression, and their methylation with a range of clinical features. As shown in Table 2, the FRPGs were correlated with age (P=0.005), cytogenetic type (P=0.031), risk_cytogenetic (P=0.001), and risk_molecular (P<0.001). The expression of SOCS1 was correlated with age (P=0.002) and FAB type (P=0.045). The methylation of SOCS1 cg03014241 CpG site was correlated with age (P=0.039), WBC count (P=0.023), PB_blast (P<0.001), and cytogenetic type (P=0.033). The FRPGs (P<0.001), AKR1C2 (P<0.001), and SOCS1 (P<0.001) expression, as well as SOCS1 methylation levels (cg03014241) (P<0.001) both affect the DFS time and OS time. Similarly, the expression of AKR1C2 (P<0.0001), the expression of SOCS1 (P<0.0001), and the SOCS1 methylation levels (cg03014241) (P=0.001) also affect the risk of FRPGs.

Table 2. Correlation between AKR1C2、SOCS1 mRNA expression/methylation、FRPGs、and clinical features in TCGA database.

Clinical features		RISK		P	AKR1C2		P	SOCS1		P	cg03014241		P	cg18841653		P
		high	low		high	low		high	low		high	low		high	low
N		76	78		77	77		77	77		77	77		77	77
SEX	Female	32 (42.1)	39 (50.0)	0.412	37 (48.1)	34 (44.2)	0.746	33 (42.9)	38 (49.4)	0.518	38 (49.4)	33 (42.9)	0.518	34 (44.2)	37 (48.1)	0.746
	Male	44 (57.9)	39 (50.0)		40 (51.9)	43 (55.8)		44 (57.1)	39 (50.6)		39 (50.6)	44 (57.1)		43 (55.8)	40 (51.9)
AGE		59.50 (15.26)	52.26 (16.14)	0.005	57.56 (15.38)	54.10 (16.66)	0.183	59.87 (14.97)	51.79 (16.22)	0.002	53.17 (16.04)	58.49 (15.76)	0.039	56.22 (16.28)	55.44 (15.96)	0.765
WBC		33.34 (36.90)	45.56 (55.31)	0.11	38.38 (40.77)	40.68 (53.42)	0.764	38.34 (47.53)	40.72 (47.51)	0.757	48.19 (52.51)	30.86 (40.12)	0.023	40.43 (45.49)	38.63 (49.48)	0.814
BM_BLAST		69.88 (18.81)	67.79 (18.93)	0.494	69.53 (18.93)	68.12 (18.84)	0.643	69.22 (19.76)	68.43 (18.00)	0.795	70.31 (17.80)	67.34 (19.83)	0.329	69.18 (19.07)	68.47 (18.73)	0.815
PB_BLAST		36.96 (32.58)	45.45 (31.65)	0.104	37.96 (32.29)	44.55 (32.16)	0.208	37.52 (32.99)	45.00 (31.32)	0.153	51.66 (30.23)	30.95 (31.10)	<0.001	40.32 (33.00)	42.14 (31.76)	0.728
FAB	M0	6 (7.9)	10 (12.8)	0.161	7 (9.1)	9 (11.7)	0.358	7 (9.1)	9 (11.7)	0.045	7 (9.1)	9 (11.7)	0.401	10 (13.0)	6 (7.8)	0.423
	M1	23 (30.3)	20 (25.6)		21 (27.3)	22 (28.6)		26 (33.8)	17 (22.1)		22 (28.6)	21 (27.3)		24 (31.2)	19 (24.7)
	M2	16 (21.1)	20 (25.6)		15 (19.5)	21 (27.3)		15 (19.5)	21 (27.3)		22 (28.6)	14 (18.2)		14 (18.2)	22 (28.6)
	M4	13 (17.1)	21 (26.9)		16 (20.8)	18 (23.4)		11 (14.3)	23 (29.9)		17 (22.1)	17 (22.1)		15 (19.5)	19 (24.7)
	M5	12 (15.8)	6 (7.7)		12 (15.6)	6 (7.8)		12 (15.6)	6 (7.8)		8 (10.4)	10 (13.0)		11 (14.3)	7 (9.1)
	M6	2 (2.6)	0 (0.0)		2 (2.6)	0 (0.0)		2 (2.6)	0 (0.0)		0 (0.0)	2 (2.6)		1 (1.3)	1 (1.3)
	M7	3 (3.9)	0 (0.0)		2 (2.6)	1 (1.3)		3 (3.9)	0 (0.0)		0 (0.0)	3 (3.9)		2 (2.6)	1 (1.3)
	Nc	1 (1.3)	1 (1.3)		2 (2.6)	0 (0.0)		1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		0 (0.0)	2 (2.6)
CYTOGENETIC TYPE	BCR-ABL1	2 (2.6)	1 (1.3)	0.031	2 (2.6)	1 (1.3)	0.679	2 (2.6)	1 (1.3)	0.358	1 (1.3)	2 (2.6)	0.033	2 (2.6)	1 (1.3)	0.32
	CBFB-MYH11	0 (0.0)	10 (12.8)		3 (3.9)	7 (9.1)		3 (3.9)	7 (9.1)		7 (9.1)	3 (3.9)		5 (6.5)	5 (6.5)
	Complex Cytogenetics	15 (19.7)	7 (9.0)		10 (13.0)	12 (15.6)		15 (19.5)	7 (9.1)		6 (7.8)	16 (20.8)		14 (18.2)	8 (10.4)
	Intermediate Cytogenetics	10 (13.2)	9 (11.5)		8 (10.4)	11 (14.3)		9 (11.7)	10 (13.0)		7 (9.1)	12 (15.6)		10 (13.0)	9 (11.7)
	MLL translocation, poor risk	1 (1.3)	2 (2.6)		1 (1.3)	2 (2.6)		2 (2.6)	1 (1.3)		3 (3.9)	0 (0.0)		1 (1.3)	2 (2.6)
	MLL translocation, t (9;11)	1 (1.3)	0 (0.0)		1 (1.3)	0 (0.0)		1 (1.3)	0 (0.0)		0 (0.0)	1 (1.3)		0 (0.0)	1 (1.3)
	N.D.	1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		2 (2.6)	0 (0.0)
	Normal Karyotype	40 (52.6)	38 (48.7)		44 (57.1)	34 (44.2)		39 (50.6)	39 (50.6)		45 (58.4)	33 (42.9)		36 (46.8)	42 (54.5)
	Poor Risk Cytogenetics	5 (6.6)	4 (5.1)		5 (6.5)	4 (5.2)		4 (5.2)	5 (6.5)		6 (7.8)	3 (3.9)		6 (7.8)	3 (3.9)
	RUNX1-RUNX1T1	1 (1.3)	6 (7.7)		2 (2.6)	5 (6.5)		1 (1.3)	6 (7.8)		1 (1.3)	6 (7.8)		1 (1.3)	6 (7.8)
RISK_CYTO	Good	1 (1.3)	16 (20.5)	0.001	5 (6.5)	12 (15.6)	0.313	4 (5.2)	13 (16.9)	0.073	8 (10.4)	9 (11.7)	0.777	6 (7.8)	11 (14.3)	0.11
	Intermediate	51 (67.1)	47 (60.3)		53 (68.8)	45 (58.4)		49 (63.6)	49 (63.6)		52 (67.5)	46 (59.7)		46 (59.7)	52 (67.5)
	N.D.	1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		2 (2.6)	0 (0.0)
	Poor	23 (30.3)	14 (17.9)		18 (23.4)	19 (24.7)		23 (29.9)	14 (18.2)		16 (20.8)	21 (27.3)		23 (29.9)	14 (18.2)
RISK_MOLECULAR	Good	1 (1.3)	16 (20.5)	0.001	5 (6.5)	12 (15.6)	0.347	4 (5.2)	13 (16.9)	0.087	8 (10.4)	9 (11.7)	0.798	6 (7.8)	11 (14.3)	0.192
	Intermediate	46 (60.5)	44 (56.4)		48 (62.3)	42 (54.5)		45 (58.4)	45 (58.4)		48 (62.3)	42 (54.5)		43 (55.8)	47 (61.0)
	N.D.	1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		1 (1.3)	1 (1.3)		2 (2.6)	0 (0.0)
	Poor	28 (36.8)	17 (21.8)		23 (29.9)	22 (28.6)		27 (35.1)	18 (23.4)		20 (26.0)	25 (32.5)		26 (33.8)	19 (24.7)
DFS_MONTHS		9.67 (10.62)	27.12 (26.82)	<0.001	9.67 (10.62)	27.12 (26.82)	<0.001	9.67 (10.62)	27.12 (26.82)	<0.001	9.67 (10.62)	27.12 (26.82)	<0.001	14.99 (19.68)	22.02 (24.13)	0.05
OS_MONTHS		14.25 (15.48)	33.97 (27.39)	<0.001	14.25 (15.48)	33.97 (27.39)	<0.001	14.25 (15.48)	33.97 (27.39)	<0.001	14.25 (15.48)	33.97 (27.39)	<0.001	20.90 (23.76)	27.58 (24.62)	0.089
Risk	high	——	——	——	50 (64.9)	26 (33.8)	<0.001	68 (88.3)	8 (10.4)	<0.001	27 (35.1)	49 (63.6)	0.001	41 (53.2)	35 (45.5)	0.42
	low	——	——		27 (35.1)	51 (66.2)		9 (11.7)	69 (89.6)		50 (64.9)	28 (36.4)		36 (46.8)	42 (54.5)
AKR1C2	high	50 (65.8)	27 (34.6)	<0.001	——	——	——	42 (54.5)	35 (45.5)	0.334	39 (50.6)	38 (49.4)	1	35 (45.5)	42 (54.5)	0.334
	low	26 (34.2)	51 (65.4)		——	——		35 (45.5)	42 (54.5)		38 (49.4)	39 (50.6)		42 (54.5)	35 (45.5)
SOCS1	high	68 (89.5)	9 (11.5)	<0.001	42 (54.5)	35 (45.5)	0.334	——	——	——	25 (32.5)	52 (67.5)	<0.001	42 (54.5)	35 (45.5)	0.334
	low	8 (10.5)	69 (88.5)		35 (45.5)	42 (54.5)		——	——		52 (67.5)	25 (32.5)		35 (45.5)	42 (54.5)
cg03014241	high	27 (35.5)	50 (64.1)	0.001	39 (50.6)	38 (49.4)	1	25 (32.5)	52 (67.5)	<0.001	——	——	——	36 (46.8)	41 (53.2)	0.519
	low	49 (64.5)	28 (35.9)		38 (49.4)	39 (50.6)		52 (67.5)	25 (32.5)		——	——		41 (53.2)	36 (46.8)
cg18841653	high	41 (53.9)	36 (46.2)	0.42	35 (45.5)	42 (54.5)	0.334	42 (54.5)	35 (45.5)	0.334	36 (46.8)	41 (53.2)	0.519	——	——	——
	low	35 (46.1)	42 (53.8)		42 (54.5)	35 (45.5)		35 (45.5)	42 (54.5)		41 (53.2)	36 (46.8)		——	——

GO and KEGG functional enrichment analysis

240 DEGs were found between high-riskand low-risk groups. The volcanic map of DEGs is shown in Figure 9A. The results of functional and pathway enrichment analysis are shown in Figure 9B, 9C. Gene Ontology (GO) analysis of DEGs showed that regulation of cell population proliferation, cell death, protein phosphorylation, cell migration, NMDA glutamate receptor clustering, and engulfment of the apoptotic cell, were mostly enriched. Among the functional and pathway enrichment analysis, ECM-receptor interaction, protein digestion, and absorption, PI3K-Akt signaling pathway, enzyme-linked receptor protein signaling pathway, transmembrane receptor protein tyrosine kinase signaling pathway, cGMP-mediated signaling, positive regulation of peroxisome proliferator-activated receptor signaling pathway, and cAMP-mediated signaling pathway were enriched most.

Figure 9. Functional and Pathway Enrichment Analysis of DEGs A. Volcanic maps of all genes. Red spot, the expression is up-regulated; Blue spot, the expression is down-regulated; Gray spot, no significantly dysregulated; The first 10 significantly differentially upregulated genes were labeled. B.The GO biological process enrichment analysis. C. The KEGG enrichment analysis.

PPI network and analysis on clusters

STRING mapped 240 DEGs into PPI network containing 154 nodes and 287 edges (Figure 10A). Then MCODE was used to discover clusters in the network. 3 clusters were selected for GO term and KEGG pathway enrichment analysis. Among them, cluster 1 contained 8 nodes and 28 edges, with the highest score (Figure 10B), cluster 2 contained 6 nodes and 14 edges (Figure 10C), cluster3 contained 5 nodes and 10 edges (Figure 10D). We then performed a functional analysis for the 3 important clusters using ClueGO. The results showed that chloride channel complex, negative regulation of chemotaxis, and extracellular matrix structural constituent conferring tensile strength were most enriched (Figure 10E).

Figure 10. A. The PPI analysis of DEGs at STRING (https://string-db.org). B、C、D. Cluster analysis of the PPI network (B. Cluster1; C. Cluster2; D. Cluster3). E. GO and KEGG analyses on the 3 hub modules.

Discussion

In this study, we screened out the FRGs in AML based on TCGA and FerrDb datasets. We found the expression of AKR1C2 and SOCS1 were significantly different in AML samples and healthy control samples and were negatively related to the OS and DFS of the AML patients. We used age as a control variable to compare the risk score based on these 2 genes with different clinical factors and found that as the same as cytogenetic stratification and molecular stratification, risk based on these 2 genes was more common than age factors and increased the occurrence of survival. Moreover, we also verified the prognostic role of AKR1C2 and SOCS1 expression in patients with intermediate molecular and cytogenetic risk status in the NCCN AML classification. The results all emphasized the promising prognostic value of AKR1C2 and SOCS1 expression in AML patients. We compared the correlation between methylation and mRNA expression of AKR1C2 and SOCS1 genes and identified methylation sites that were negatively correlated with mRNA expression of these two genes. We also explored the prognostic role of methylation sites.

Aldo-keto reductases family 1 member C2 (AKR1C2), a member of the aldo/keto reductase superfamily. Previous studies have shown that its expression is selectively deleted in prostate and breast cancer and leads to clonal expansion of tumor cells. In addition, some studies have shown that its expression is negatively correlated with apoptosis-inducing factor (AIF), and is involved in the metastasis of liver cancer and the occurrence of NSCLC drug resistance [18–20]. Suppressor of cytokine signaling 1 (SOCS1) is identified as a tumor suppressor gene. Silencing of SOCS1 expression due to promoter methylation has been found in several malignancies such as pancreatic cancer, breast cancer, lymphoma, and leukemia [21–25]. There is a study mentioned that SOCS1 methylation leads to SOCS1 gene silencing, which promotes AML cell growth and proliferation by inhibiting the downstream JAK2/STAT signaling pathway [26]. However, in our study, low expression of SOCS1 suggested a better prognosis. The specific role of SOCS1 methylation in AML still needs to be further investigated. Besides, no studies have been conducted on the prognostic effect of AKR1C2 and SOCS1 on AML. In this study, we systematically investigated the relationship between AKR1C2 and SOCS1 expression and survival using the TCGA database and found that low AKR1C2 and SOCS1 expression was highly correlated with more favorable OS and DFS in AML patients. Multivariate Cox regression further confirmed this conclusion, suggesting that low expression of AKR1C2 and SOCS1 is a strong prognostic factor for OS and DFS in AML patients. Finally, we further evaluated the overall prognostic value of ARK1C2 and SOCS1 in patients with a molecular and cytogenetic risk status of intermediate risk in the NCCN AML classification, proving that decreased expressions of ARK1C2 and SOCS1 are indeed associated with better OS and DFS in AML patients. Moreover, the expressions of ARK1C2 and SOCS1 are possible to further refine the risk stratification in clinic and to help accurately stratify patients with intermediate-risk. In conclusion, our analysis highlights that ARK1C2 and SOCS1 are promising biomarkers for predicting prognosis in patients with AML.

Increasing evidence has proved that abnormal DNA methylation plays a crucial role in the occurrence and development of cancers [27,28]. In our analysis, we first determined whether the methylation status of AKR1C2 and SOCS1 would affect their mRNA expression through Pearson correlation and Spearman's rho. There was a significant negative correlation between cg18841653 and cg03014241 CpG sites methylation and ARK1C2 and SOCS1 mRNA expression in AML. This negative correlation can well explain the low expression of AKR1C2 and SOCS1 in AML. We then investigated the prognostic significance of two selected CpG sites in AKR1C2 and SOCS1 DNA methylation and found that SOCS1 hypermethylation was strongly associated with better OS and PFS in AML patients. In summary, AKR1C2 and SOCS1 are negatively regulated by their methylation, and SOCS1 methylation status may be a strong indicator of good OS and PFS.

Ferroptosis is a new type of iron-dependent programed cell death. Its main mechanism is catalyzing the highly expressed unsaturated fatty acids on the cell membrane under the action of ferric divalent or esteroxylase, resulting in lipid peroxidation, thus inducing cell death [12]. Ferroptosis involves in a variety of diseases, including inhibition of tumor cell growth [29–31]. In our study, we divided patients into low and high-risk groups based on the expression of FRPGs and identified DEGs of the two groups. Functional enrichment analysis showed that DEGs enriched in regulation of cell death, unsaturated fatty acid metabolic process, PI3K-Akt signaling pathway, ECM-receptor interaction, regulation of cAMP-mediated signaling, and cGMP-mediated signaling pathway. Significantly, the metabolism of unsaturated fatty acids is associated with ferroptosis. PI3 K signaling pathway and ECM-receptor interaction are associated with cell differentiation, proliferation, and apoptosis [32,33]. The cAMP and cGMP-mediated signaling pathways are involved in regulating a variety of cellular physiological processes. However, although we have explored the biological processes involved in FRPGs in AML through enrichment analysis, how the expression of FRPGs and their methylation affect the prognosis of AML still needs to be further explored in biomedical experiments. Nevertheless, the role of FRPGs in the stratification of prognosis in AML is now clear. The more precise prognostic stratification, especially in patients with intermediate-risk is encouraging, and FRPGs are prognosis biomarkers worthy of further exploration.

Availability of data and materials

The clinical data, transcriptional data, along with methylation profiles of patients with AML are available in the TCGA database via the Cbioportal website (https://www.cbioportal.org/study/summary?id=laml_tcga_pub). The FRGs data are available in the FerrDb website (http://www.zhounan.org/ferrdb/index.html).

Acknowledgments

We acknowledge TCGA and FerrDb database for providing their platforms and contributors for uploading their meaningful datasets.

Disclosure statement

The authors declare that there is no conflict of interest regarding the publication of this paper.

Additional information

Funding

This work was supported by the [Key Medical Projects of Jiangsu Province #1] under Grant [number BE2017747] and [Jiangsu Social Development Project#2] under Grant [number BE2018711].