Ferroptosis-based molecular prognostic model for adrenocortical carcinoma based on least absolute shrinkage and selection operator regression.

BACKGROUND: This study aimed to find ferroptosis-related genes linked to clinical outcomes of adrenocortical carcinoma (ACC) and assess the prognostic value of the model.
METHODS: We downloaded the mRNA sequencing data and patient clinical data of 78 ACC patients from the TCGA data portal. Candidate ferroptosis-related genes were screened by univariate regression analysis, machine-learning least absolute shrinkage, and selection operator (LASSO). A ferroptosis-related gene-based prognostic model was constructed. The effectiveness of the prediction model was accessed by KM and ROC analysis. External validation was done using the GSE19750 cohort. A nomogram was generated. The prognostic accuracy was measured and compared with conventional staging systems (TNM stage). Functional analysis was conducted to identify biological characterization of survival-associated ferroptosis-related genes.
RESULTS: Seventy genes were identified as survival-associated ferroptosis-related genes. The prognostic model was constructed with 17 ferroptosis-related genes including STMN1, RRM2, HELLS, FANCD2, AURKA, GABARAPL2, SLC7A11, KRAS, ACSL4, MAPK3, HMGB1, CXCL2, ATG7, DDIT4, NOX1, PLIN4, and STEAP3. A RiskScore was calculated for each patient. KM curve indicated good prognostic performance. The AUC of the ROC curve for predicting 1-, 3-, and 5- year(s) survival time was 0.975, 0.913, and 0.915 respectively. The nomogram prognostic evaluation model showed better predictive ability than conventional staging systems.
CONCLUSION: We constructed a prognosis model of ACC based on ferroptosis-related genes with better predictive value than the conventional staging system. These efforts provided candidate targets for revealing the molecular basis of ACC, as well as novel targets for drug development.

INTRODUCTION

Programmed cell death has been shown to be a significant type of cell death. It acts as a natural barrier to prevent cells from developing into cancers. , Dysregulation of programmed cell death signaling pathways is emerging as a key factor in tumorigenesis.  The most thoroughly studied aspect of programmed cell death is apoptosis. Research has revealed new mechanisms of programmed cell death, one of which is ferroptosis. The concept of ferroptosis was first proposed by Stockwell et al. in 2012, and it is a non‐apoptotic programmed cell death process. Recent studies have focused on the role of ferroptosis in the progression, invasion, migration, and cell death of multiple types of cancers. , , For most anti‐cancer drugs, activation of programmed cell death pathways to kill tumor cells is a vital anti‐tumor mechanism. Due to the acquired and intrinsic resistance of tumor cells to apoptosis, the therapeutic efficacy of inducing apoptosis in tumor is limited.  Therefore, the use of other forms of non‐apoptotic cell death to clear tumor cells and control the proliferation of drug‐resistant cell clones provides a new therapeutic possibility. The potential of targeting ferroptosis in cancer treatment has generated high expectations. , ,

Adrenocortical carcinoma (ACC) is an isolated malignant tumor, which has attracted more and more attention since the end of the last century. It is a rare and highly aggressive malignant disease and can occur at any age. Localized tumors can be cured by surgery. Even if the tumor has been completely removed, however, recurrence is common. Unlike other tumors, treatment options after ACC recurrence are limited. , ,  The prognosis remains poor. Most studies have shown that the median survival time of ACC patients is about 12 months. It has been thought that changes in the Wnt / β‐Catenin and IGF‐2 signaling pathways lead to ACC, but recent studies have shown that these changes are not sufficient to cause the occurrence of malignant adrenal tumors. ,  Therefore, the mechanism of the development and occurrence of ACC remains incompletely understood, and numerous genes and their functions remain to be discovered and explained. , ACC shares some genetic profiles that are associated with promising therapeutic responsiveness in other cancers.  With the development of precision medicine, we have the opportunity to identify genes that are related to clinical outcomes and novel molecular targets for new drugs. A genomics‐guided clinical care approach offers the potential for prolonging life expectancy and also improving the quality of life for ACC patients.

In this study, we aimed to find candidates ferroptosis genes, which were related to clinical outcomes of ACC. We constructed a prognosis model of ACC based on ferroptosis‐related genes and then clarified the prognostic value of ferroptosis genes in ACC. These efforts may contribute to the development of better treatment strategies in the future.

METHODS

Data acquisition

We downloaded the RNA‐sequencing data and clinical data for 78 ACC patients from the TCGA data portal (https://tcga‐data.nci.nih.gov/tcga/dataAccessMatrix.htm). Regulator genes and marker genes for ferroptosis (ferroptosis‐related genes) were downloaded from the FerrDb database, and articles were downloaded from the PubMed database.

Candidate gene screening and validation, prediction model establishment

Two steps were involved in the candidate gene screening. First, we performed univariate regression analysis of every ferroptosis‐related gene and overall survival. Genes with p‐values < 0.05 were included in the next step. Univariate Cox regression was carried out using the “survival” R package. Then, machine‐learning least absolute shrinkage and selection operator (LASSO) were used to select independent risk factors that affected outcomes. LASSO Cox regression was implemented using the “glmnet” R package. Correlation coefficients at lambda.min were chosen for the final model, and cross‐validation was used to tune and optimize the LASSO penalty terms. K‐fold cross‐validation (k = 5) was used to train and test the model.

After candidate genes were selected at lambda.min, a prognostic model was then constructed using the formula below. RiskScore was then calculated for each patient.

Assessing the effectiveness of prediction models

We grouped the patients into high‐ and low‐risk groups based on the median riskScore. The KM curve for these data was used to compare the prognosis between high‐risk and low‐risk groups according to the riskScore. Receiver operating characteristic (ROC) curves and areas under the curve (AUCs) were calculated with the “survivalROC” and “survminer” R packages to demonstrate the predictive ability of riskScore for 1‐, 3‐, and 5‐year OS. A flow diagram of this trial is shown in Figure 1.

FIGURE 1

Flowchart of the experiment

External validation was done using the GSE19750 cohort. Data were downloaded from the GEO database. The riskScore was calculated using the formula mentioned above. The clinical data were also downloaded. We determined the ROC curve and the Kaplan–Meier curve to test the predictive value of the prognostic model.

We generated nomogram by combining the riskScore value and clinic‐pathological factors to predict survival probability at 1, 3, and 5 years. This is a quantitative and intuitive method to assess the association between variables and survival. We then measured the prognostic accuracy by calculating the Harrell's concordance index (C‐index). The larger the C‐index, the more accurate the prognostic prediction proved to be.  We compared the prediction model with conventional staging systems using the C‐index. We assessed calibration by comparing observed and predicted survival probabilities using the KM method and applied bootstraps with 100 replicates Nomogram was undertaken using the “rms” R package.

Functional analysis

We used Gene Ontology analysis (GO) to identify characteristic biological attributes of survival‐associated ferroptosis genes and performed Kyoto Encyclopedia of Genes and Genomes pathway (KEGG) enrichment analysis to identify functional attributes. GO and KEGG analysis was done using the following R packages: “DOSE” “org. Hs.eg.db”, “clusterProfiler” and “pathview”. For visualization of the data, the “ggplot2” package was used.

RESULTS

The RNA‐sequencing data and clinical data of 78 ACC patients were downloaded from TCGA database. Two patients were excluded from the analysis due to missing clinical information. Of those who were qualified for inclusion, 48 were female and 28 were male. The average overall survival time was 3.39 ± 2.69 years. Two hundred fifty‐nine ferroptosis genes were downloaded from the FerrDb database and Pubmed database (123 marker genes, 109 suppressor genes, and 150 driver genes).

First, we performed univariate regression analysis of every ferroptosis‐related gene and overall survival. Seventy genes were identified as survival‐associated ferroptosis‐related genes with p < 0.05. Figure 2A shows the HR level of each survival‐associated ferroptosis‐related genes.

FIGURE 2

Parameter selection. (A) Forest map of the univariate regression analysis. The horizontal axis represents the Hazard ratio (HR). The horizontal ordinate represents each gene with a p‐value < 0.05 in univariate regression analysis. (B) and (C) Tuning parameter selection using LASSO with k‐fold cross‐validation (k = 5)

Next, LASSO Cox regression was implemented for these 70 genes. Correlation coefficients at lambda.min were chosen for the final model (Figure 2B, C, optimal lambda.min =0.078). After fivefold cross‐validation, 17 genes were included in the final model. The Coef level for each gene is shown in Table 1. RiskScore was also calculated for each patient (Table 2).

TABLE 1

Seventeen genes included in the model and its corresponding Coef

Gene	Coef
STMN1	0.006855766
RRM2	0.003733332
HELLS	0.017375996
FANCD2	0.00161208
AURKA	0.007796465
GABARAPL2	−0.0054616
SLC7A11	0.016787224
KRAS	0.014846229
ACSL4	−0.021912674
MAPK3	−0.008147927
HMGB1	0.013098853
CXCL2	0.006211908
ATG7	−0.005985336
DDIT4	0.00576449
NOX1	−0.007679209
PLIN4	0.000928894
STEAP3	0.002633784

TABLE 2

RiskScore and clinical stage for each patient

	OS	Event	RiskScore	Risk	T	N	M	Stage
TCGA.OR.A5J2	1677	1	26.65612951	High	t3	n0	m1	Stage iv
TCGA.OR.A5J3	1942	0	−82.55082586	Low	t3	n0	m0	Stage iii
TCGA.OR.A5J5	365	1	220.2545248	High	t4	n0	m0	Stage iii
TCGA.OR.A5J6	2428	0	−60.07161748	Low	t2	n0	m0	Stage ii
TCGA.OR.A5J7	490	1	127.9296784	High	t3	n0	m0	Stage iii
TCGA.OR.A5J8	579	1	181.7346398	High	t3	n0	m0	Stage iii
TCGA.OR.A5J9	1183	0	53.85394279	High	t2	n0	m0	Stage ii
TCGA.OR.A5JA	922	1	20.83685542	High	t4	n0	m1	Stage iv
TCGA.OR.A5JB	551	1	249.1943157	High	t4	n0	m1	Stage iv
TCGA.OR.A5JD	2782	0	−87.26295608	Low	t2	n0	m0	Stage ii
TCGA.OR.A5JE	2105	1	37.11735056	High	t1	n0	m0	Stage i
TCGA.OR.A5JF	1259	0	0.159417811	Low	t2	n0	m0	Stage ii
TCGA.OR.A5JG	541	1	50.49723047	High	t4	n1	m1	Stage iv
TCGA.OR.A5JI	1424	0	−84.30563802	Low	t1	n0	m0	Stage i
TCGA.OR.A5JJ	309	0	79.28220068	High	t4	n1	m1	Stage iv
TCGA.OR.A5JK	1255	0	−13.39745347	Low	t4	n0	m1	Stage iv
TCGA.OR.A5JL	670	0	−124.2939039	Low	t1	n0	m0	Stage i
TCGA.OR.A5JM	562	1	46.32584829	High	t4	n0	m1	Stage iv
TCGA.OR.A5JO	889	0	30.20226513	High	t1	n0	m0	Stage i
TCGA.OR.A5JP	149	0	94.40884423	High	t2	n0	m0	Stage ii
TCGA.OR.A5JQ	674	0	−77.31526038	Low	t2	n0	m0	Stage ii
TCGA.OR.A5JR	3688	0	−130.3421379	Low	t1	n0	m0	Stage i
TCGA.OR.A5JS	383	0	29.70127434	High	t2	n0	m0	Stage ii
TCGA.OR.A5JT	488	0	−61.35190065	Low	t2	n0	m0	Stage ii
TCGA.OR.A5JV	1541	0	−95.23738127	Low	t2	n0	m0	Stage ii
TCGA.OR.A5JW	1924	0	8.031815994	Low	t2	n0	m0	Stage ii
TCGA.OR.A5JX	950	0	98.19251924	High	t3	n0	m0	Stage iii
TCGA.OR.A5JY	552	1	63.85014088	High	t4	n1	m1	Stage iv
TCGA.OR.A5JZ	211	0	−76.63893854	Low	t2	n0	m0	Stage ii
TCGA.OR.A5K0	1029	0	30.6211023	High	t2	n0	m0	Stage ii
TCGA.OR.A5K1	2723	0	−41.34566511	Low	t2	n0	m0	Stage ii
TCGA.OR.A5K2	994	1	94.27576045	High	t4	n0	m0	Stage iii
TCGA.OR.A5K3	2842	0	−69.84016404	Low	t2	n0	m0	Stage ii
TCGA.OR.A5K4	528	0	−52.32932887	Low	t4	n0	m0	Stage iii
TCGA.OR.A5K5	253	0	27.93880869	High	t3	n0	m0	Stage iii
TCGA.OR.A5K6	1130	0	1.900433368	Low	t2	n0	m0	Stage ii
TCGA.OR.A5K8	504	0	40.08479735	High	t2	n0	m0	Stage ii
TCGA.OR.A5K9	344	1	100.4672018	High	t2	n0	m0	Stage ii
TCGA.OR.A5KO	1414	0	20.61064176	Low	t4	n0	m1	Stage iv
TCGA.OR.A5KT	2673	0	10.91838006	Low	t1	n0	m0	Stage i
TCGA.OR.A5KU	4673	0	18.75358502	Low	t2	n0	m0	Stage ii
TCGA.OR.A5KV	3659	0	41.91500616	High	t2	n1	m0	Stage iii
TCGA.OR.A5KW	1525	0	−23.1429177	Low	t2	n1	m0	Stage iii
TCGA.OR.A5KX	1091	0	115.6707949	High	t2	n1	m0	Stage iii
TCGA.OR.A5KY	391	1	130.5067414	High	t4	n1	m1	Stage iv
TCGA.OR.A5KZ	125	1	218.2568428	High	t2	n0	m0	Stage ii
TCGA.OR.A5L3	3897	0	−10.86225829	Low	t1	n0	m0	Stage i
TCGA.OR.A5L4	724	0	−262.5860277	Low	t4	n0	m0	Stage iii
TCGA.OR.A5L5	840	0	−51.65798978	Low	t1	n0	m0	Stage i
TCGA.OR.A5L6	628	0	33.8583362	High	t2	n0	m0	Stage ii
TCGA.OR.A5L8	555	0	29.68040353	High	t2	n0	m0	Stage ii
TCGA.OR.A5L9	645	0	−38.37036871	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LA	487	0	−75.732104	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LB	1204	1	80.98199258	High	t4	n0	m1	Stage iv
TCGA.OR.A5LC	159	1	198.3267946	High	t4	n0	m1	Stage iv
TCGA.OR.A5LD	1197	1	68.83890176	High	t4	n0	m0	Stage iii
TCGA.OR.A5LE	662	1	64.07674666	High	t2	n0	m0	Stage ii
TCGA.OR.A5LG	1589	0	27.02398016	High	t3	n0	m0	Stage iii
TCGA.OR.A5LH	2385	1	−58.9701129	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LJ	1105	1	49.85160852	High	t2	n1	m1	Stage iv
TCGA.OR.A5LK	2222	0	−44.29088098	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LL	1613	1	24.1344422	High	t2	n0	m0	Stage ii
TCGA.OR.A5LM	1858	0	−14.74304601	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LN	1916	0	−93.96683194	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LO	1949	0	147.1147267	High	t2	n0	m0	Stage ii
TCGA.OR.A5LP	1583	0	−175.3190167	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LR	639	0	−87.18232022	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LS	882	0	−8.787969716	Low	t2	n0	m0	Stage ii
TCGA.OR.A5LT	365	0	−31.61654757	Low	t3	n0	m0	Stage iii
TCGA.OU.A5PI	709	0	12.48413456	Low	t2	n1	m1	Stage iv
TCGA.P6.A5OF	207	1	227.6880918	High	t4	n0	m0	Stage iii
TCGA.P6.A5OG	383	1	119.2763823	High	t4	n0	m1	Stage iv
TCGA.PA.A5YG	470	0	−99.68999465	Low	t2	n0	m0	Stage ii
TCGA.PK.A5H8	3240	0	−72.91560771	Low	t2	n0	m0	Stage ii
TCGA.PK.A5H9	307	0	−85.04336783	Low	t2	n0	m0	Stage ii
TCGA.PK.A5HA	830	0	−72.43647323	Low	t1	n0	m0	Stage i

OS, overall survival in days. Events indicate survival status. 1 represents patient was dead. 0 represents patient was alive. The patients were classified into low‐risk group and high‐risk group according to the median value of the risk scores.

Figure 3 showed that patients with poorer prognosis had lower riskScores. For patients who died during the follow‐up, the average riskScore was −25.51 (SD = 74.47), while patients who survived follow‐up had an average riskScore of 80.84 (SD = 101.83). It is clear that these groups were significantly different with regard to RiskScores (p = 1.21E‐07, Figure 3).

FIGURE 3

RiskScores of patients with different survival statuses during follow‐up. 0 representing death and 1 representing survival

The median riskScore was 19.68 for all patients. Patients were grouped into high‐ and low‐risk groups based on their riskScores. The high‐risk group (riskScore >19.68) had 37 patients, and the low‐risk group (riskScore ≤19.68) had 39. KM curve showed that the high‐risk group had poorer prognoses (p < 0.0001, Figure 4A). Then, we determined the time‐dependent ROC curve to find the prognostic performance of riskScore for survival prediction. The AUC of the ROC curve for predicting 1‐, 3‐, and 5‐year(s) survival time was 0.975, 0.913, and 0.915 respectively (Figure 4B–D).

FIGURE 4

(A) KM survival analysis of high‐ and low‐risk groups. Yellow curve represents high‐risk patient group; blue curve represents low‐risk patient group. (B–D). Time‐dependent ROC analysis for the prognostic model to predict 1‐ and 3‐, and 5‐ year(s) survival. Area under the curve (AUC) values are shown

Data from the GSE19750 cohort were used to perform external validation of the predictive value of the model. Consistent with the results in the TCGA cohort, patients in the high‐risk group had significantly poorer survival probability than the low‐risk group (p = 0.011, Figure 5A). The AUCs for 1‐year, 3‐year, and 5‐year OS were 0.765, 0.773, and 0.805, respectively (Figure 5B–D).

FIGURE 5

External Validation. (A). Nomogram predicting survival probability. (B–D). Time‐dependent ROC analysis for the prognostic model to predict 1‐ and 3‐, and 5‐ year(s) survival using the GSE19750 cohort. Area under the curve (AUC) values are shown in the figure

We constructed the nomogram prognostic evaluation model to predict the 1‐, 3‐, or 5‐year OS time in patients by combining riskScores and pathological information (Figure 5A). The predictive accuracy of 1‐, 3‐, or 5‐year OS is shown in Figure 5B–D. The C‐index of the nomogram was 0.92 (se(C)=0.02). We also compared the prediction model with conventional staging systems. The C‐index for the TNM staging system was 0.75 (se(C)=0.05), which was lower than that of our model. Thus, our prognostic prediction model had better predictive ability.

Figure 6 shows the GO (Figure 6A) and KEGG (Figure 6B) analyses of survival‐associated ferroptosis genes. KEGG analysis showed that the genes were mostly enriched in central carbon metabolism in cancer, cellular senescence, and the NOD‐like receptor signaling pathway.

FIGURE 6

Functional analysis. (A) GO analysis. X‐axis represents three types of GO. The node size is representative of gene count level, and the color represents – log 2 (p‐value). MF: molecular function. CC: cellular component. BP: biological process. (B) KEGG analysis. X‐axis represents gene count. Y‐axis represents pathway involved in the analysis. The color represents – log 10 (p‐value)

DISCUSSION

Adrenocortical carcinoma is a highly malignant cancer with limited therapeutic options. Patients usually exhibit lymph node and distant metastases by the time of diagnosis. Surgery is the primary treatment strategy, while adjuvant therapies are frequently needed. Mitotane is currently the only agent approved. For advanced ACC, a combination of mitotane with a cytotoxic regimen of etoposide, doxorubicin, and cisplatin (EDP‐M) is recommended. However, a narrow therapeutic window and endocrine side effects restrict the clinical use of these drugs. ,  Thus, there is an urgent need to identify drug targets and develop new therapeutic strategies to treat ACC.

High‐throughput biotechnology such as genomics provides a good entry point for basic medicine to clinical medicine. Prognostic and predictive biomarkers selected from high‐throughput genomic data are of critical importance in cancer management.  The question of how to mine valuable information efficiently from vast biological sequences is crucial to researchers. Meanwhile, traditional variable‐selecting methods such as multivariate regression analysis are insufficient when facing big data. LASSO, a regularization method, is a promising solution. LASSO is particularly attractive in prognostic studies due to its capabilities of regression coefficients shrinkage and automatic variable selection. LASSO has been successfully applied in prognostic model studies. , In this study, we focused on candidate ferroptosis genes related to prognosis of ACC for the first time. We constructed a prognosis model based on 17 survival‐associated ferroptosis‐related genes using the machine‐learning method. These efforts may contribute to the development of better treatment strategies in the future. We found that the predictive value of our model is better than that of the conventional staging system. Our study provided a handful of candidate targets for revealing the molecular basis of ACC, as well as novel targets for drug development.

Recent studies have demonstrated that ACC is sensitive to ferroptosis, indicating that induction of ferroptosis could be a promising treatment approach. Therefore, we constructed a prognostic model including 17 survival‐associated ferroptosis‐related genes. Belavgeni's study showed direct inhibition of glutathione peroxidase 4, a key factor in the initiation of ferroptosis, in human ACC NCI‐H295R cells leading to high necrotic populations. High STMN1 expression has been observed in aggressive ACC patients. , Ikeya's recent study shows that overexpression of AURKA, a gene identified in our study, can cause atypical mitosis in adrenocortical carcinoma with the p53 somatic variant.  The p53 protein, an important regulator of ferroptosis, is frequently mutated in ACC. ACSL4, which has been reported to dictate ferroptosis sensitivity by shaping cellular lipid composition, is demonstrated to be highly expressed in mouse adrenal glands.

In our study, ferroptosis gene riskScores showed good predictive value. Nomograms have been well developed as a prognostic assessment tool and proven to be more accurate than conventional staging systems in several cancers. , ,  We constructed a nomogram by combining ferroptosis gene riskScores and clinic‐pathological factors. Our model showed better predictive value than the conventional staging system, a finding supported by C‐index (0.92) and calibration curve. In terms of precision medicine, our model has potential clinical applications.

There are some possible weaknesses in this study. We performed internal validation using k‐fold cross‐validation and bootstrap resampling methods. External and multicenter prospective cohorts with large sample sizes are still needed to validate the clinical application of our model, and basic research needs to be done to clarify the underlying mechanism.

In conclusion, our study identified candidate ferroptosis genes, which were related to clinical outcomes of ACC. We constructed a prognosis prediction model of ACC based on ferroptosis‐related genes. Our model showed better predictive value than the conventional staging system. These efforts provided a handful of underlying targets for revealing the molecular basis of ACC, as well as for drug development.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

AUTHOR CONTRIBUTIONS

Lin, Liang, and Hu contributed to the literature search and the design of the study. Lin and Liang analyzed and interpreted the data. Lin and Liang wrote the study Lin and Sun formatted the figures and tables. Sun revised the article Hu helped perform the analysis with constructive discussions. The final study was approved by all the authors.