Establishing a novel Fanconi anemia signaling pathway-associated prognostic model and tumor clustering for pediatric acute myeloid leukemia patients

2.1 Gene expression and clinical information

Gene expression matrix of the TPM (transcripts per million) type and a set of clinical traits were obtained from The Cancer Genome Atlas (TCGA)-LAML cohort using the “TCGAbiolinks” R package. We analyzed the expression pattern of each FA gene by combining datasets within a genotype-tissue expression (GTEx) and TCGA-LAML cohort. The data of TARGET-AML were obtained from the official website (https://ocg.cancer.gov/). We performed a statistical analysis using the wilcox.test or kruskal.test. The “GEOquery” R package was utilized to obtain the datasets of the GSE71014, GSE37642, and GSE12417 cohorts within the gene expression omnibus (GEO) database. Then, the expression matrix of TARGET-AML (n = 257), TCGA-LAML (n = 137), GSE71014 (n = 104), GSE37642 (n = 136), and GSE12417 (n = 78) was processed using the R language, and the batch correction was performed using the “sva” package. We then used the prcomp() function for the principal component analysis (PCA). The BEAT-AML cohort (n = 399) [22,23,24,25,26] was utilized as an external validation set.

2.2 Tumor clustering

We used the “ConsensusClusterPlus” R package [27] to cluster TARGET-AML samples according to the expression characteristics of FA signatures. Then, PCA was conducted. Additionally, we utilized the “survival” R package to perform the cluster-related prognostic analysis.

2.3 Lasso regression modeling construction

We performed Lasso regression modeling analysis using the datasets of TARGET-AML, TCGA-LAML GSE71014, GSE12417, and GSE37642 cohorts. The cases were divided into a training set (75%) and an internal validation set (25%). We used the “glmnet” R package to conduct a Lasso regression modeling analysis of the training set. The resulting model provided gene coefficients and risk score values. Risk scores were grouped according to their median values to determine high- and low-risk groups. The OS risk values for the internal validation set were also obtained. We plotted the risk profile, survival curves, and the survival status map. The hub gene expression profile was visualized as a heat map. Furthermore, we used the BEAT-AML cohort to externally validate the above FA-associated prognostic model.

2.4 Gene enrichment

Differentially expressed genes in two tumor clusters or risk groups were identified using a “limma” R package with a false discovery rate of 0.05 and a log2 fold change (log2FC) of 1. A “ggplot2” R package was used to generate an MA plot. Also, a heat map of FA genes and the clinical traits of age, gender, French–American–British (FAB), and clustering status were created. Gene enrichment analysis was performed, and the results were visualized using the “GOplot” R package [28]. Furthermore, gene set enrichment analysis (GSEA) was conducted by a “clusterProfiler” R package [29,30], and the results were visualized using the “gseaplot()” function of the enrichplot R package. Common members of differential genes between the risk and clustering groups were obtained, and the related gene ontology (GO)-Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed.

2.5 Cox regression

To assess the OS prognosis, the univariate/multivariate Cox regression analyses were conducted using a “survival” R package, considering the factors of prognostic risk score, FAB, gender, and age. Additionally, we performed external validation of the BEAT-AML cohort based on the multivariate Cox regression model associated with FA, using the predict() function. The forest plots were generated using the plot() R function. To evaluate the receiver operating characteristic (ROC) results for different survival times, we utilized the “survivalROC” R package and provided the value of area under the curve (AUC).

2.6 Nomogram prediction model

A cph () modeling analysis was performed using the “regplot” R package, and a Nomogram was obtained. The calibration curves of AML cases were plotted using the calibrate() function of the “ggstatsplot” R package. The C index value was obtained using the “rms” R package. In addition, the external validation dataset (BEAT-AML) was used to perform calibration curve analysis based on the FA Nomogram. The net reclassification improvement (NRI) and integrated discrimination improvement (IDI) analyses were conducted using the “survIDINRI” R package. Finally, the decision curve analysis (DCA) was performed using the “stdca.R” package from www.decisioncurveanalysis.org.

3.1 Analytic strategy

The schematic diagram of our analytic strategy is presented in Figure 1. We initially analyzed the expression profile of 23 FA genes and the clinical traits. Next, using these FA signatures, we identified two tumor subtypes (cluster 1/2) in AML samples by utilizing the clustering function of the “ConsensusClusterPlus” R package. Additionally, we extracted and performed batch correction on the expression matrix of TARGET-AML, TCGA-LAML GSE71014, GSE12417, and GSE37642 cohorts. This was done with FA-related genes, the FA expression matrix, and corresponding clinical traits. An internal validation and training set were randomly selected from the samples. The training set-based Lasso regression modeling analysis was then conducted. We conducted multivariate Cox regression analyses. Nomograms and related calibrate curves were plotted as well. We assessed the clinical effectiveness of FA models by DCA. Additionally, survival, differential, and enrichment analyses of GO, KEGG, and GSEA were performed after risk or cluster grouping. FA-related models were validated in both the internal and external validation cohorts.

Figure 1

A schematic diagram of the study design.

3.2 Expression profile

We analyzed the expression features of FA genes based on the combined datasets from the GTEx and TCGA-LAML cohorts.

As shown in Figure S1, we observed increased expression levels in the tumor group for FANCA, BRCA2, FANCG, FANCM, SLX4, ERCC4, and FAAP100 (all p < 0.05), while a decreased level was observed for FANCB, FANCC, FANCD2, FANCE, FANCI, FANCL, RAD51C, RAD51, BRCA1, UBE2T, XRCC2, MAD2L2, and RFWD3 (all p < 0.05) compared with the normal group. Additionally, there were expression differences for all FA genes, except FANCG, among different FAB types (Figure 2, all p < 0.05). We found a slightly higher expression of RAD51 (Figure S2, p < 0.05) and UBE2T (p < 0.05) in the male group. Furthermore, as shown in Figure S3, we observed positive correlations between the age factor and expression levels of FANCC and FANCE (p < 0.001).

Figure 2

Correlation between FA gene expression and FAB. Box plots display each FA gene’s expression difference among different FAB types for the TARGET-AML cohort. A kruskal.test() was performed. *p < 0.05, **p < 0.01, ***p < 0.001.

3.3 Tumor clustering analysis

Next, we conducted a tumor clustering analysis based on the expression profile of FA genes within TCGA-LAML (Figure 3a) and TARGET-AML (Figure 3b) datasets. We separated these AML clusters using a PCA approach for both cohorts. However, we found a statistical difference in OS survival between the two clusters only for the TARGET-AML (Figure 3b, p = 1.153 × 10⁻³) but not for the TCGA-LAML (Figure 3a). Specifically, the pediatric AML patients of TARGET in cluster 2 showed a worse OS prognosis compared to those in cluster 1 (Figure 3b). We provided heat map data in Figure 3c, showing the association between FA gene expression and clinical information. Furthermore, a significant difference was observed in all FA gene expressions between clusters 1 and 2 (Figure S4, all p < 0.05).

Figure 3

Tumor clustering analysis of FA genes. (a) AML patients from the TCGA-LAML were clustered according to the expression characteristics of FA-related genes. A PCA was performed to evaluate the classification effect on tumor clusters. The FA-specific prognostic curve of OS is also provided. (b) Similar analyses were conducted for the TARGET-AML cohort. (c) The correlation between FA gene expression and various clinical traits, including age, gender, FAB, and clusters, was analyzed. A heat map is provided. (d and e) The differential gene identification of the two clusters was completed, and a volcano plot is provided (d). The GO/KEGG (e) and GSEA (f) gene enrichment analyses were then conducted.

The differential genes of the two clusters were obtained, and a volcano plot is shown in Figure 3d. The GO/KEGG enrichment analysis data (Figure 3e) revealed that these genes were associated with the events related to the “cell cycle,” “DNA replication,” and peroxide-related events such as “hydrogen peroxide metabolic process,” and “peroxidase activity.” The GSEA result further showed a series of genes linked to DNA repair events, such as “nucleotide excision repair” and “homology-directed repair,” as well as cell cycle events, such as “M_phase” and “cell cycle checkpoints” (Figure 3f). These findings suggest a connection between the pediatric AML FA subtypes and DNA repair, cell cycle, and peroxide-related metabolic events.

3.4 Lasso regression model

We initially conducted batch matching on the TARGET-AML, TCGA-LAML GSE71014, GSE12417, and GSE37642 cohorts (Figure 4a). Subsequently, we used 75% of the samples as a training set. Using the expression feature of FA genes, we performed a Lasso regression modeling analysis (Figure 4b). Our model identified five hub genes with correlation coefficients: BRIP1 (0.284), FANCC (−0.031), FANCL (0.022), MAD2L2 (0.270), and RFWD3 (−0.068). The results, including the expression/risk profile, survival curve, and survival status of the training and internal validation sets, are presented in Figure 4c–d. We observed an increased mortality status in AML cases as the risk value increased. Furthermore, the high-risk group exhibited a worse OS prognosis than the low-risk group (Figure 4c, p = 7.358 × 10⁻⁵; Figure 4d, p = 7.249 × 10⁻³).

Figure 4

Lasso regression modeling analysis of FA. (a) We performed the Lasso regression modeling analysis based on the expression matrix of FA-related genes in TARGET-AML, TCGA-LAML GSE71014, GSE12417, and GSE37642 cohorts. A PCA plot shows the batch correction of the expression matrix. (b) Based on the training set, the Lasso regression modeling analysis was performed. (c and d) The data of the hub gene expression profile, risk profile, survival status map, and survival curve are provided in the training set (c) and internal validation set (d).

3.5 Prognostic analysis

In the Cox regression analyses of TCGA-LAML, the age factor was found to be correlated with worse clinical OS prognosis in adult AML patients (Figure 5a, HR > 1, p < 0.001). For the pediatric AML cases in the TARGET-AML cohort, we observed a relationship between higher risk scores and worse clinical OS prognosis (Figure 5b, all HR > 10, p < 0.001). The relationship between hub gene expression and clinical traits was visualized as a heatmap in Figure 5c. Furthermore, we performed ROC analysis to assess the OS predictive value of our model. Figure 5d shows that the AUC values were greater than 0.62 for 1-, 3-, 5-year survival time. These suggested that our FA model exhibits good predictive value for OS in pediatric AML patients.

Figure 5

Prognostic analyses of TCGA-LAML and TARGET-AML. Targeting the factors (prognostic risk score, FAB, gender, and age), we performed the univariate/multivariate Cox regression analyses to evaluate the prognosis of OS for the TCGA-LAML (a) and TARGET-AML (b). We created a heatmap that combines the expression of hub genes with related clinical traits (c). Based on the risk score, we conducted the ROC analyses to assess the OS prognosis of TARGET-AML for different survival time points (d).

3.6 Nomogram

We conducted the Nomogram and related assessments. Among the hub genes, there were higher expression levels of BRIP1 (Figure 6a, p < 0.001), FANCL (p < 0.05), and MAD2L2 (p < 0.001) in the high-risk group, compared with the low-risk group. By combining the risk score, gender, age, and FAB, we established a Nomogram to predict the OS rates of a pediatric AML case. The Nomogram, shown in Figure 6b, predicted the 1-, 3-, and 5-year survival rates at 0.34, 0.292, and 0.0881, respectively. The calibration plot curve in Figure 6c showed a high overlap between the predicted and observed value lines. The C index value of 0.66 was obtained (Figure 6c). Furthermore, NRI and IDI results suggested that the predictive function of the model improved after incorporating the risk score factor (Figure 6d, all IDI > 0, NRI > 0, p < 0.05). The DCA results (Figure 6e) further showed a better clinical validity of the “age + gender + FAB + risk score” group than other groups.

Figure 6

Nomogram and related assessment analyses. (a) The expression difference of FA hub genes between high- and low-risk groups is presented. (b) A Nomogram is constructed to forecast the 1-, 3-, and 5-year survival rates of a given AML case. Furthermore, the outcomes of the calibration plot curve (c), NRI/IDI assessment (d), and DCA (e) are provided.

3.7 External validation analysis

Based on the FA-associated prognostic models mentioned above, we conducted external validation using the BEAT-AML cohort. As shown in Figure 7a, it can be observed that AML cases with a high FA-associated risk score showed a poorer prognosis (p = 3.612 × 10⁻²). Furthermore, an external validation of the multivariate Cox regression model associated with FA was conducted, revealing a favorable predictive value for AML patients, particularly for the 1-year survival time (AUC = 0.71), as depicted in Figure 7b. A calibration plot curve was also generated using the external validation set, illustrating a substantial concurrence between the predicted and observed value lines (Figure 7c). These results serve to substantiate the good predictive capability of our FA-associated prognostic models.

Figure 7

External validation of FA-associated prognostic models. (a) Based on the Lasso regression modeling of FA genes, we performed the survival curve analysis using the external validation cohorts (BEAT-AML). (b) The ROC results of the external validation set for 1-, 2-, and 3-year survival times were obtained based on the FA-associated multivariate Cox regression model. (c) Based on the FA Nomogram, we also performed the calibration curve analysis of the external validation set.

3.8 Prognostic risk-related differential gene analysis

We identified a series of differential genes (n = 65) based on the high- and low-risk grouping. The results were visualized via an MA plot (Figure 8a) and a heat map (Figure 8b). Further analysis of these differential genes was conducted using GO and KEGG enrichment approaches (Figure 8c). Our findings revealed the biological metabolic processes, such as the “antibiotic metabolic process” and “cellular oxidant detoxification.” Additionally, GSEA results showed the enrichment of events, including “DNA repair,” “G2M checkpoint,” and “oxidative phosphorylation” (Figure 8d). We also obtained 24 common members of differential genes between the risk and clustering groups (Figure 8e) and performed the gene enrichment analyses of GO-KEGG (Figure 8f). This analysis highlighted peroxide-related events, including the “hydrogen peroxide catabolic process” and “peroxidase activity.”

Figure 8

Prognostic risk-related differential gene analysis. (a and b) The differential gene identification between the high- and low-risk groups was performed, and the MA plot (a) and heat map (b) are provided to visualize the data. Gene enrichment analyses of GO-KEGG (c) and GSEA (d) were conducted to gain further insights. (e and f) Common members of differential genes between the risk and clustering groups were obtained, and gene enrichment analyses of GO-KEGG were performed to explore the functional significance.