Department of Bioinformatics, School of Biomedical Engineering and Informatics, Nanjing Medical University, Nanjing, Jiangsu 211166, China
2.
Department of Hematology of the Affiliated Huai'an No. 1 People's Hospital of Nanjing Medical University, Northern Jiangsu Institute of Clinical Medicine, Huai'an, Jiangsu 223300, China
3.
Collaborative Innovation Center for Personalized Cancer Medicine, Jiangsu Key Lab of Cancer Biomarkers, Prevention and Treatment, Nanjing Medical University, Nanjing, Jiangsu 211166, China
4.
Department of Hematology, the First Affiliated Hospital of Nanjing Medical University, Jiangsu Province Hospital, Nanjing, Jiangsu 210029, China
5.
Key Laboratory of Hematology of Nanjing Medical University, Nanjing, Jiangsu 210029, China
6.
Department of Pharmacology, School of Basic Medical Sciences, Nanjing Medical University, Nanjing, Jiangsu 211166, China
7.
The Affiliated Cancer Hospital of Nanjing Medical University, Jiangsu Cancer Hospital, Jiangsu Institute of Cancer Research, Nanjing, Jiangsu 210002, China
Wentao Liu, Department of Pharmacology, School of Basic Medical Sciences, Nanjing Medical University, 101 Longmian Avenue, Nanjing, Jiangsu 211166, China. E-mail: painresearch@njmu.edu.cn
Qianghu Wang and Kening Li, Department of Bioinformatics, School of Biomedical Engineering and Informatics, Nanjing Medical University, 101 Longmian Avenue, Nanjing, Jiangsu 211166, China. E-mails: wangqh@njmu.edu.cn (Wang) and likening@njmu.edu.cn (Li)
Given the extremely high inter-patient heterogeneity of acute myeloid leukemia (AML), the identification of biomarkers for prognostic assessment and therapeutic guidance is critical. Cell surface markers (CSMs) have been shown to play an important role in AML leukemogenesis and progression. In the current study, we evaluated the prognostic potential of all human CSMs in 130 AML patients from The Cancer Genome Atlas (TCGA) based on differential gene expression analysis and univariable Cox proportional hazards regression analysis. By using multi-model analysis, including Adaptive LASSO regression, LASSO regression, and Elastic Net, we constructed a 9-CSMs prognostic model for risk stratification of the AML patients. The predictive value of the 9-CSMs risk score was further validated at the transcriptome and proteome levels. Multivariable Cox regression analysis showed that the risk score was an independent prognostic factor for the AML patients. The AML patients with high 9-CSMs risk scores had a shorter overall and event-free survival time than those with low scores. Notably, single-cell RNA-sequencing analysis indicated that patients with high 9-CSMs risk scores exhibited chemotherapy resistance. Furthermore, PI3K inhibitors were identified as potential treatments for these high-risk patients. In conclusion, we constructed a 9-CSMs prognostic model that served as an independent prognostic factor for the survival of AML patients and held the potential for guiding drug therapy.
Acute myeloid leukemia (AML) is a malignant clonal disorder characterized by myeloid blast proliferation with expansion and a block in differentiation[1–2]. Chemotherapy has been the predominant treatment for AML for decades. Although most patients initially respond to chemotherapy, approximately 75% of the patients relapse and die from the disease within five years of diagnosis[3]. Additionally, the molecular heterogeneity among AML patients may lead to varying treatment outcomes[4–6]. Therefore, it is essential to identify molecular markers at the time of diagnosis to predict the risk of treatment failure or relapse, as well as survival outcomes of AML patients[7].
Several clinical features are shown to be associated with AML prognosis, including age, white blood cell count, and a history of hematological diseases[8]. On the other hand, certain biological and cytogenetic characteristics of leukemia cells, such as internal tandem repeats of the FLT3 gene, as well as t(8;21), inv(16), and 3q abnormalities, may also help predict prognosis[9]. Although these factors have been applied to patient prognostic stratification, some patients still do not have identifiable European LeukemiaNet (ELN) 2022 guideline-defined mutations or cytogenetic events. As such, the prognostic value of these mutations or cytogenetic characteristics is limited[10]. Therefore, the risk stratification and prognosis evaluation of AML patients need to be further improved.
Cell surface markers (CSMs) are proteins present on the surfaces of cells, serving as crucial identifiers that distinguish and categorize various cell types. Because of their extracellular accessibility and involvement in regulating cellular processes, such as metabolite transport and intercellular communication, CSMs have been widely used as diagnostic and prognostic markers in cancers[11]. For example, CSMs have been used to identify cancer stem cells, whose persistence after chemotherapy causes cancer relapse, and for immunophenotyping in cancers[12]. Numerous studies suggest that CSMs play a crucial role in the prognostic assessment and prediction of disease progression. Moreover, it has been reported that over 66% of drugs approved for treating human diseases target CSMs, based on the data from the DrugBank database[13]. The identification of cancer-specific CSMs may not only suggest disease progression and stratify patients into different prognostic categories, but also provide targets for drug therapy.
In AML, CSMs have been reported to be applicable for screening, diagnosis, prognosis, predicting treatment response, and monitoring disease progression. Notably, several studies showed that CSMs, such as CD34, CD38, CD96, and CD47, might have prognostic significance to further guide treatment decisions[14]. However, these studies primarily examined the prognostic significance of individual CSM in AML and frequently relied on a small size of patient samples[15] and low throughput technologies such as flow cytometry[16], immunohistochemistry, and quantitative polymerase chain reaction[17]. Because of the complexity of AML and its extremely high inter-sample heterogeneity, the expression level of a single CSM may have a limited classification power for risk stratification. Therefore, a comprehensive computational study evaluating the prognostic potential of human CSMs in AML patients may provide more precise information for risk stratification and prognostic predictions for AML patients.
To comprehensively evaluate the prognostic efficiency of all CSMs in AML patients, we conducted multi-model analyses on high-throughput transcriptome data using the LASSO, Adaptive LASSO, and Elastic Net algorithms. Consequently, we developed a 9-CSM prognostic model with potential applications in risk stratification and prognostic prediction for AML patients. Furthermore, this model was validated at the protein level, to provide a novel perspective on prognostic prediction and therapeutic decision-making for AML patients.
Materials and methods
Data acquisition and processing
The information on CSMs was downloaded from the Cell Surface Protein Atlas (CSPA, http://wlab.ethz.ch/cspa/), a free online database that describes 2886 proteins constituting the human in silico surfaceome.
The clinical and microarray data (GSE147515 contains data from both 198 normal subjects and 1534 AML patients, while GSE10358, GSE12417, GSE71014, GSE84334, GSE106291, GSE193094, GSE146173, and GSE183817 contain data from 91, 74, 104, 38, 235, 6, 233, and 7 AML patients, respectively) were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). The gene expression profiles and their corresponding clinical information of the BeatAML were downloaded from the Beat acute myeloid leukemia database (BeatAML, http://www.vizome.org/aml/), which includes 239 AML patients. The proteomics data related to BeatAML were downloaded from the Genomic Data Commons Data Portal (GDC, https://pdc.cancer.gov/pdc/). The proteomics data related to The Cancer Genome Atlas (TCGA)-Acute Myeloid Leukemia-like (LAML) were downloaded from the deep-scale proteomic and phosphoproteomic database of AML (https://proteomics.leylab.org/). The clinical information and the gene expression profiles of the TCGA-LAML were downloaded from TCGA (https://portal.gdc.cancer.gov/). We used the average gene expression level for multiple probes corresponding to one single gene. The information on 130 patients for the TCGA training data and the GEO external validation datasets is shown in Supplementary Table 1 (available online).
The expression and function analyses of CSMs in AML
We used the limma package to identify genes that were differentially expressed between AML and normal samples, calculated the mean absolute deviation, and performed principal component analysis (PCA) using the R in-built function "prcomp" to obtain the matrix of principal components. The "enrichr" function implemented in the "enrichR" package was used to identify significantly enriched biological pathways and processes.
Training and construction of the 9-CSMs prognostic model
The Elastic Net, LASSO, and Adaptive LASSO regression methods were implemented using the "glmnet" R package (version 4.1.3) to select candidate features. The LASSO proposed by Tibshirani et al is a technique used to filter out genes with a significant association with patient survival[18]. The LASSO transforms based on the formula below,
ˆβlasso=argminβ(12n‖
(1)
Then, separately, we subjected the CSMs to the Adaptive LASSO regression, which was proposed by Zou et al[19]. The Adaptive LASSO applies different weights for different parameters, as shown below,
\hat{\beta}^{ {adaplasso }}={}_{\beta}^{arg min }\left(\frac{1}{2 n}\left\|y- X \beta\right\|_2^2+\lambda {\sum}_{j=1}^p \omega_j\left|\beta_j\right|\right).
(2)
The Elastic Net regression method proposed by Hastie combines the benefits of \lambda_1 and \lambda_2 regularizations[20], as shown below,
\begin{split}
\hat{\beta }^{ { elasticNet }}= {}_{\beta}^{arg min }\left(\frac{1}{2 n}\left\|y- X \beta\right\|_2^{2 }+\lambda_1 \|\beta\|_1+\lambda_2\|\beta\|_2^2\right).
\end{split}
(3)
The union of the candidate features selected from each method was then subjected to an automatic model averaging and multi-model selection analysis implemented using the "glmulti" function. The best model was selected based on the minor Akaike information criterion (AIC) value. Finally, the prognostic risk score was calculated using the following equation:
score=\sum_{i=1}^N\beta_i\times Gene_i,
(4)
where βi is the CSM model coefficient obtained from the best model (Supplementary Table 2, available online) and Genei is the gene expression level of each CSM in the prognostic model.
We performed a Kaplan-Meier survival analysis of the overall survival (OS) and event-free survival (EFS) of AML patients using the "Survival" and "Survminer" R packages. For each dataset, the patients were classified into high-risk and low-risk groups according to the median CSMs score, and the significance was evaluated by the log-rank test. To evaluate the predictive performance of our prognostic model, we performed receiver operating characteristic (ROC) curve analysis using the "timeROC" R package (version 0.4) for 1-year, 3-year, and 5-year survival. The area under the curve (AUC) was used to estimate the accuracy of the model.
Evaluation of the independence of the 9-CSMs prognostic model
We performed univariable Cox regression analysis using the "coxph" function to evaluate the prognostic score of the 9-CSMs model. Furthermore, we performed multivariable Cox regression analysis on variables, including the 9-CSMs prognostic score, age, sex, race, cytogenetic risk, molecular subtype, bone marrow blast percentage, PB blast percentage, WBC count, and gene mutations, to examine their potential association with the survival outcome of AML patients in the TCGA cohort.
Construction and calibration of the nomogram model
For the nomogram model, we integrated independent prognostic factors obtained from multivariable Cox regression analysis together with survival time data and event data to fit the Cox proportional hazards model using the "survival" package. The nomogram model for predicting the OS rate of patients with AML was constructed by "foreign" and "rms" R packages. Finally, calibration plots were generated for respective survival times, with the x-axis showing predicted survival and the y-axis showing actual OS.
Gene-gene interaction network analysis
To unravel the gene-gene interactions associated with 9-CSMs, we used the GeneMANIA Database (http://genemania.org/), a comprehensive database that contains information about known and predicted protein-protein interactions, to find genes that may share functions.
Pre-processing and analysis of single-cell RNA sequencing (scRNA-seq) data
The scRNA-seq data was downloaded from the OMIX (https://ngdc.cncb.ac.cn/omix, accession No. OMIX002180). Data integration, unsupervised clustering, and visualization were conducted using the Seurat framework[21]. We identified the cellular states of the leukemia-like cells according to the article by Li et al[22]. The risk score was computed through AUCell.
Drug sensitivity analysis
To evaluate the use of our prognostic model in guiding medication for patients with AML, we examined the association between anticancer drug sensitivity and the prognostic model risk score. We obtained the half maximal inhibitory concentration (IC50) data and the cell line expression data from the Genomics of Drug Sensitivity in Cancer database (GDSC, https://www.cancerrxgene.org/). The Spearman's rank correlation test was performed to compute the correlation coefficient and P-values between the IC50 values and the prognostic model risk score. A Spearman's rank adjusted P-value < 0.05 was considered statistically significant.
Prediction of immunotherapy response
We used the online tool Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) to predict the response of AML patients to the existing immune checkpoint inhibitors (ICIs) targeting PD-1 and CTLA-4[23–24]. We inputted the gene expression matrix of AML patients, and downloaded the corresponding analysis results from the website. The Spearman's rank correlation test was used to compute the correlation coefficient and P-values between the TIDE score and the 9-CSMs risk score.
Results
Construction of the 9-CSMs prognostic model based on multi-model analysis
To identify reliable predictive markers for the prognosis of AML patients, we designed a bulk RNA-seq-based workflow to comprehensively dissect the characteristics of AML patients (Fig. 1). We obtained 2886 CSMs from the CSPA database and then identified 209 CSMs that were differentially expressed in AML (n = 1534) and healthy control (n = 198) samples (Fig. 2A, Supplementary Tables 3 and 4 [available online]). The results of PCA showed that the differentially expressed genes mentioned above effectively distinguished AML from normal tissues (Fig. 2B). Further functional enrichment analysis revealed that CSMs up-regulated in AML were enriched in biological processes associated with the inflammatory response, cytokine-cytokine receptor interaction, IL-6/JAK/STAT3 signaling pathway, and IL-2/STAT5 signaling (Fig. 2C), indicating that the up-regulated cell-surface proteins in AML may be associated with the immune response and tumor development.
Figure
1.
Schematic of the workflow for identifying and constructing the CSM prognostic model.
Firstly, we obtained 2886 cell-surface markers from the Cell Surface Protein Atlas (CSPA) database, and 1965 cell-surface markers were expressed in the integrated GSE147515 dataset. After performing differential gene expression analysis and univariable Cox regression analysis, 53 markers were retained. We conducted feature selection using multi-model (Adaptive LASSO regression, LASSO regression, and Elastic Net) to identify 9-CSMs linked to AML prognosis. Furthermore, the predictive value of the 9-CSMs risk score was validated at the transcriptome and proteome levels. Finally, the 9-CSMs risk score was used to establish a nomogram model. We performed functional enrichment analysis and identified potential drugs for patients with a high-risk 9-CSMs score.
Figure
2.
Differential gene expression and functional enrichment analysis between AML and normal samples.
A: A volcano plot of differentially expressed genes between AML and normal samples. Each dot represents a differentially expressed gene, with yellow indicating up-regulation in AML and green indicating down-regulation. B: PCA plot of differentially expressed genes between AML (n = 1534) and normal (n = 198) samples. Each dot represents a differentially expressed gene, with yellow indicating up-regulation in AML and green indicating down-regulation in normal samples. C: Enrichment bar plots for pathways enriched by upregulated cell surface markers in AML for GO_Biological_Process and KEGG pathway databases. Abbreviations: AML, acute myeloid leukemia; PCA, principal component analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.
The univariable Cox regression analysis was performed to evaluate the association between gene expression levels of these differentially expressed CSMs and OS of 130 AML patients in the TCGA-LAML cohort (Supplementary Table 5, available online). We then performed multi-model feature selection analysis using three regression analysis methods: LASSO regression analysis, Adaptive LASSO, and Elastic Net. We implemented this analysis using a 10-fold cross-validation maximum likelihood penalty estimator in the TCGA-LAML training cohort. The optimal λ values selected by each method are shown in Supplementary Fig. 1 (available online). In addition, CSMs selected by the three methods were combined to perform an automated model average and multi-model selection analysis based on the smallest AIC value (Supplementary Tables 6–9, available online). Using the multivariable Cox regression, the best model (9-CSMs) was determined (Supplementary Table 2, available online). The results of the stepwise multivariable Cox regression analysis showed that the identified markers, except for CD160, were independent prognostic factors for AML patient survival (Fig. 3A). The details of the 9-CSMs model are shown in Supplementary Table 10 (available online). Furthermore, CXCR2,LY6E, SUSD3, S1PR5, and IL1RAP were identified as risk factors, while CD160, CLEC5A, NTNG2, and SLC4A1 were identified as protective factors in the prognostic model.
Figure
3.
Construction and evaluation of the 9-CSMs prognostic models in the training datasets.
A: A forest plot of a multivariable Cox regression analysis of the identified 9-CSMs. Events represent the number of deaths. Akaike information criterion (AIC) is a criterion to measure the goodness of fit of a statistical model. The concordance index was used to evaluate the predictive ability of the model. The log-rank test was used for hypothesis testing. B: Distribution of survival data in the TCGA-LAML training cohort. The 130 patients from the TCGA training cohort were stratified into two groups according to the median of the risk scores. C: The mRNA heatmap of the 9-CSMs in the TCGA-LAML training cohort. D and E: Kaplan-Meier survival curves for the 9-CSMs prognostic model for OS and EFS. The high 9-CSMs risk and low risk groups each consisted of 65 patients. P-values were calculated by log-rank test. F: The ROC curves for 1-year, 3-year, and 5-year overall survival predictions for the 9-CSMs prognostic model. *P < 0.05, **P < 0.01, and ***P < 0.001. Abbreviations: CI, confidence interval; TCGA-LAML, The Cancer Genome Atlas-acute myeloid leukemia; AUC, area under the curve.
The 9-CSMs prognostic model predicted outcomes in the training dataset
To assess the prognostic performance of the model, we assigned a risk score calculated from the 9-CSMs prognostic model to each patient in the TCGA-LAML training cohort (Supplementary Table 11, available online) and classified these patients into the low-risk and high-risk groups based on the median 9-CSMs score (Fig. 3B). In patients with high-risk scores, each risk factor was highly expressed, while each protective factor was highly expressed in patients with low-risk scores (Fig. 3C). The 9-CSMs prognostic model significantly predicted survival outcomes for both OS and event-free survival of 130 AML patients in the Kaplan-Meier survival analysis (Fig. 3D and 3E), demonstrating that patients with high-risk scores had shorter survival times than those with low risk scores in the TCGA-LAML training cohort.
Additionally, we assessed the predictive performance of the 9-CSMs prognostic model by calculating the area under the time-dependent ROC in the TCGA-LAML training cohort. The AUC for OS at 1, 3, and 5 years was 0.84, 0.87, and 0.93, respectively, indicating that the model for the TCGA-LAML training cohort performed well (Fig. 3F).
Validation of the 9-CSMs prognostic model at the transcriptome and proteome levels
We validated the predictive performance of the 9-CSMs prognostic model in three independent validation datasets, including GSE10358, GSE12417, and GSE71014 (Supplementary Tables 12–14, available online). Patients were classified into the low-risk and high-risk groups according to the median 9-CSMs scores (Fig. 4A). The Kaplan-Meier survival analysis consistently revealed that patients had significantly poorer survival outcomes in the high-risk group than in the low-risk group (P = 0.0007 for GSE10358, P = 0.036 for GSE12417, and P = 0.078 for GSE71014, Fig. 4B), indicating that patients with high-risk scores had shorter survival than those with low-risk scores in all three datasets. We then performed time-dependent ROC analysis to further assess the consistency of the predictive ability of the 9-CSMs prognostic model in the three independent datasets and found that the 9-CSMs prognostic model accurately predicted the survival outcomes of AML patients (Fig. 4C).
Figure
4.
Evaluation of the 9-CSMs prognostic model in multiple validation cohorts and its prognostic value at the protein level.
A: Distribution of the risk score in the GSE10358, GS12417, and GSE71014 validation cohorts. Patients are stratified into two groups based on the median of the 9-CSMs risk scores. The x-axis represents patients ordered by their risk scores, with each point corresponding to an individual patient. The y-axis shows the risk scores. B: The Kaplan-Meier survival curves for the 9-CSMs prognostic model for overall survival in the GSE10358, GS12417, and GSE71014 validation cohorts. C: The receiver operating characteristic curves for the 1-year, 3-year, and 5-year overall survival predictions for the 9-CSMs prognostic model in the GSE10358, GS12417, and GSE71014 validation cohorts, respectively. The overall survival of patients in the GSE10358 and GSE12417 datasets was less than 5 years. D–H: Scatter plots demonstrate the correlation between mRNA and protein abundances of CSMs in the TCGA-LAML data. I: Kaplan-Meier survival curves for overall survival according to the CSMs score in the TCGA-LAML cohort at the protein level. J: Kaplan-Meier survival curves for disease-free survival according to the CSMs score in the Beat acute myeloid leukemia (BeatAML) cohort at the protein level. P-values were calculated by log-rank test (B, I, and J). Abbreviation: CSMs, cell surface markers.
Only five genes, including CLEC5A, CXCR2, IL1RAP, LY6E, and SLC4A1, had available proteomic data in the TCGA-LAML dataset, and the results showed a significant correlation between their mRNA and protein abundances (Fig. 4D–4H). Consequently, we calculated the risk score using the five available proteins in the TCGA-LAML dataset (33 AML patients) and found that the CSMs risk score effectively stratified patient outcomes at the protein level (Fig. 4I). We also found a trend towards poorer prognosis in patients with higher 9-CSMs risk scores in another independent cohort, BeatAML, which includes 167 AML patients (Fig. 4J).
The 9-CSMs risk score was an independent prognostic factor for AML patients
To evaluate whether the 9-CSMs prognostic model was an independent prognostic factor for patient survival, we performed univariable and multivariable Cox regression analyses on clinicopathological variables, along with the 9-CSMs risk score in the TCGA-LAML training cohort. The results showed that the 9-CSMs risk score, TP53, DNMT3A, RUNX1, and age were independent prognostic factors for AML patients (Table 1). These independent prognostic factors (i.e., age, TP53, DNMT3A, RUNX1, and 9-CSMs risk score) were further used to construct a nomogram to facilitate AML patient prognosis prediction. The corresponding score for each factor was calculated in the nomogram and the total score might be used as a tool for survival prediction for AML patients (Supplementary Fig. 2A, available online). The prediction results of the nomogram calibration curve of OS were consistent with the observation results of AML patients (Supplementary Fig. 2B).
Table
1.
Univariable and multivariable Cox regression analyses for independent prognostic factors in 130 TCGA-LAML patients
In both the BeatAML (239 patients) and TCGA-LAML (130 patients) cohorts, patients in the adverse risk group had a higher 9-CSMs risk score than those in the favorable and intermediate-risk groups (Fig. 5A–5D). Furthermore, we found that the 9-CSMs risk score did not significantly differ between 34 NPM1 mutated and 93 wild-type AML patients (Fig. 5E), and that the 9-CSMs risk score was shown to effectively predict prognosis in AML patients, regardless of their NPM1 mutation status (Fig. 5F and 5G). Additionally, we compared the 9-CSMs risk score between de novo and secondary AML patients using the GSE183817 (4 de novo and 3 secondary AML patients) and GSE146173 (202 de novo and 31 secondary AML patients) datasets and found that there was no significant difference in the 9-CSMs risk score between patients classified as de novo or secondary AML (Supplementary Fig. 3A and 3B, available online). A higher 9-CSMs score in patients with secondary AML had a relatively poor prognosis but was not statistically significant, possibly because only 31 secondary AML samples were present in the GSE146173 (Fig. 5H). These findings indicate that the 9-CSMs prognostic model may be applicable to both de novo and secondary AML patients.
Figure
5.
Evaluation of the 9-CSMs prognostic model in different AML risk stratifications and subtypes.
A and B: Boxplot showing the distribution of the 9-CSMs risk score across ELN 2017 (A) and ELN 2008 (B) risk categories in 239 BeatAML patients. C and D: Boxplot showing the 9-CSMs risk score of AML patients with different cytogenetic (C) and molecular (D) risks in the TCGA dataset. E: Boxplot showing the 9-CSMs risk score of AML patients with NPM1 mutation or wild-type in the TCGA dataset. Wilcoxon rank-sum test was used to measure the differences between groups. F and G: Kaplan-Meier analysis of overall (F) and event-free (G) survival of 130 TCGA-LAML patients. "+" indicates patients carrying NPM1 mutations or a high 9-CSMs risk score, while "–" represents wild-type or a low 9-CSMs risk score. Patients were classified into four groups according to NPM1 mutation and 9-CSMs risk score. H: Kaplan-Meier survival curves for overall survival according to the 9-CSMs score. Patients from the GSE146173 cohort were classified into two groups based on the median risk score. I: Boxplot showing the 9-CSMs risk scores among FAB classification of AML patients in the TCGA-LAML dataset. The Kruskal-Wallis test was used to measure the differences between groups. J and K: Kaplan-Meier survival curves for overall survival in non-M3 patients (J) and M3 (K) patients according to the 9-CSMs prognostic model score. Non-M3 patients and M3 patients from the TCGA-LAML cohort were classified into two groups based on the median risk score. Data are presented as the median and interquartile range (A–E and I). P-values were calculated by Wilcoxon rank-sum test (A–E) or log-rank test (F–H, J, and K). Abbreviations: AML, acute myeloid leukemia; CSMs, cell surface markers; BeatAML, beat acute myeloid leukemia; TCGA-LAML, The Cancer Genome Atlas-acute myeloid leukemia; ELN, European LeukmiaNet; FAB, French-American-British classification system.
Acute promyelocytic leukemia is classified as the M3 subtype by the French-American-British classification system[25]. This subtype is sensitive to differentiation induction therapy containing all-trans retinoic acid and arsenic trioxide[26]. In the current study, we found that 14 AML patients with M3-type had the lowest 9-CSMs risk score (Fig. 5I). To exclude the influence of patients with M3 type, we performed survival analysis in 116 non-M3 AML patients and found that the 9-CSMs risk score might still be used for prognostic stratification of these patients (Fig. 5J). Notably, the 9-CSMs risk score also provided prognostic stratification in M3 patients (Fig. 5K). These results indicate that the prognostic role of 9-CSMs risk score may be independent of AML subtypes.
Patients with high 9-CSMs risk scores were resistant to chemotherapy
To further characterize patients with high 9-CSMs risk scores, we classified the patients in the TCGA-LAML training cohort according to the median 9-CSMs risk score and performed a functional enrichment analysis. We observed eight up-regulated KEGG pathways in patients with high 9-CSMs risk scores (Fig. 6A). Notably, one of these enriched signaling pathways, the PI3K-AKT signaling pathway, has been reported to be frequently activated in AML patient blasts and contributes to drug resistance of these cells[27]. Additionally, we performed a network analysis based on physical interactions and co-expression and identified that 9-CSMs gene modules were associated with leukocyte migration (Fig. 6B). It is noteworthy that leukocyte migration was reported to be associated with extramedullary metastasis in AML[28], and one study has proposed a correlation between extramedullary infiltration in AML and a poorer prognosis[29].
Figure
6.
Functional enrichment and gene co-expression network analysis of the 9-CSMs prognostic model.
A: Enrichment bar plots for pathways enriched by 65 patients with high 9-CSMs risk scores in AML for KEGG_Human_2021 pathway databases. B: The functional modules involved by the 9-CSMs prognostic model based on GeneMANIA.
Based on the prognostic value of 9-CSMs, we hypothesized that the 9-CSMs prognostic model might help predict the treatment outcome of patients. To test this hypothesis, we investigated the association between 9-CSMs risk scores and drug resistance using real-world medication data of AML patients.
Firstly, we derived malignant cells from published AML scRNA-seq datasets and analyzed the association between 9-CSMs risk scores in these cells and their sensitivity to chemotherapy, the prevailing treatment for AML[22]. The results showed that malignant cells that refractory to chemotherapy (n = 7982) had higher 9-CSMs risk scores than those at the initial diagnosis (n = 9557) (Fig. 7A and 7B). We also found that quiescent stem-like cells (QSCs, n = 1785) had a higher 9-CSMs risk score than proliferating stem/progenitor-like cells (PSPs, n = 2344) (Fig. 7C and 7D), which was consistent with a previous research showing that a subpopulation of QSCs was involved in chemoresistance and poor outcome in AML compared with PSPs[22].
Figure
7.
Relationship between 9-CSMs risk score and drug resistance in real-world medication data of AML patients.
A: UMAP visualization of different responses to chemotherapy in PT3 malignant cells. B: Boxplots of 9-CSMs risk scores at diagnosis and refractory for PT3 malignant cells. C: UMAP visualization of PSP and QSC based on the definition of previous studies from PT3. D: Boxplot of 9-CSMs risk scores for PSP and QSC based on previous studies defined in PT3. E: Boxplots of 9-CSMs risk scores for initial treatment-sensitive or treatment-resistant groups in GSE106291. F: Boxplots of 9-CSMs risk scores for untreated and AraC-resistant groups in GSE193094. G: Boxplot showing the 9-CSMs risk score in the decitabine (DAC) treated group in bulk transcriptomic data (GSE84334) of patient-derived AML samples. This dataset includes 18 patients who responded to decitabine and 20 patients who did not respond to decitabine. H: The scatterplot showing the correlation between the 9-CSMs risk score (x-axis) and the TIDE score (y-axis). Data are presented as the median and interquartile range (B and D–G). P-values were calculated using the Student's t-test (B) or Wilcoxon rank-sum test (D–F). *P < 0.05 and ***P < 0.001. Abbreviations: AML, acute myeloid leukemia; Dx, diagnosis; Refr, refractory; QSC, quiescent stem-like cells; PSP, proliferating stem/progenitor-like cells; AraC, cytarabine; PT3, Patient 03.
Furthermore, we found that the 9-CSMs score was higher in the chemotherapy-resistant group (n = 71) than in the chemotherapy-sensitive group in the GSE106291 dataset (n = 164) (Fig. 7E). Additionally, the 9-CSMs score was higher in the cytarabine-resistant group (n = 3) than in the group without treatment (n = 3) in the GSE193094 dataset (Fig. 7F). We subsequently evaluated the response of patients with high 9-CSMs scores to the hypomethylating agent, decitabine. Patients who did not respond to decitabine (n = 20) had higher 9-CSMs risk scores than those who responded to decitabine (n = 18) in the GSE84334 dataset, but this difference was not statistically significant (Fig. 7G). We also predicted the response of AML patients to the existing ICIs targeting PD-1 and CTLA-4 using the online TIDE tool, which integrates both intrinsic tumor cell characteristics and the tumor microenvironment to assess the response to immunotherapy. The result showed a significant positive correlation between the 9-CSMs risk score and the TIDE score (Fig. 7H), indicating that patients with a high 9-CSMs risk score tended to be resistant to ICIs.
In summary, patients with high 9-CSMs scores may respond poorly to conventional chemotherapy, hypomethylating agents, and the existing ICI treatment. These results indicate that the 9-CSMs prognostic model may be a valuable tool for predicting patient response to the treatments.
Potential drug prediction in patients with higher 9-CSMs scores
To further determine the clinical relevance of the 9-CSMs risk score in terms of AML patient treatment, we analyzed the association between the CSM prognostic score and drug sensitivity to small-molecule drugs using the drug IC50 data and cell line expression data from the GDSC database. The results showed that patients with a high prognostic risk score were more sensitive to Torin 2, WYE-125132, and GSK1059615 (Fig. 8A and 8B, Supplementary Fig. 4A–4C [available online]). Torin 2 and WYE-125132 act as inhibitors of mTOR[30–31], a downstream target of the PI3K/AKT pathway, while GSK1059615 is a dual inhibitor of PI3K and mTOR[32]. In gastric cancer, GSK1059615 was reported to block the PI3K-AKT-mTOR cascade activation; and the administration of GSK1059615 in nude mice significantly suppressed the growth of subcutaneous AGS xenografts[32]. Combined with the aforementioned results that TCGA-AML patients with high 9-CSMs risk scores exhibited a high activation of the PI3K-AKT signaling pathway, we believe that the 9-CSMs prognostic model may play a role in guiding the selection of chemotherapy drugs for patients.
Figure
8.
Prediction of potential drugs in patients with higher 9-CSMs scores.
A: The volcano plot of correlation between Spearman correlation coefficient and –lg(P-values). Drugs with a significant correlation to the 9-CSMs risk score are labeled. B: Correlations between GDSC drug sensitivity values (IC50) and the 9-CSMs risk score for the most sensitive drugs. The top 30 drugs were considered. Abbreviations: CSMs, cell surface markers; GDSC, Genomics of Drug Sensitivity in Cancer.
The diagnostic and screening advantages of CSFs have been found in various cancers, including AML. However, CSFs for predicting the prognosis of AML patients still need to be explored. Here, we constructed a prognostic model based on nine CSMs (SLC4A1, CLEC5A, IL1RAP, CD160, S1PR5, LY6E, NTNG2, CXCR2, and SUSD3) that might provide a reliable prognosis value for clinical application. Our results indicate that the risk score may stratify AML patients into two risk groups exhibiting significantly divergent survival outcomes.
Firstly, we examined the difference in the expression of CSMs between AML patients and healthy controls, which was associated with overall survival in AML patients. The functional enrichment analysis showed that the differentially up-regulated CSMs were correlated with the inflammatory response, IL-6/JAK/STAT3 signaling pathway, and IL-2/STAT5 signaling pathway. These pathways were reported to be dysregulated in many cancers, indicating that CSMs play a crucial role in cellular immunity. Secondly, we constructed a prognostic model based on multi-model analysis, and this 9-CSMs prognostic score was found to be an independent prognostic factor for AML patients. We also evaluated the performance of the prognostic model and validated its accuracy at both the transcriptome and proteome levels. Thirdly, paired single-cell data showed that refractory samples in malignant cells had higher risk scores than those at the initial diagnosis. Finally, we examined the correlation between the risk score of our prognostic model and drug sensitivity. We found that AML patients in a high-risk category were more likely to respond to treatments, such as Torin 2, WYE-125132, and GSK1059615, demonstrating the clinical importance of the 9-CSMs prognostic model in the treatment selection of AML patients in the future. Collectively, we constructed an effective prognostic model that might significantly contribute to prognosis prediction and serve as a crucial reference in clinical treatment.
The use of immunotherapies, especially ICIs, has demonstrated effects in enhancing anti-tumor responses in certain solid tumors and hematological cancers[33]. However, several studies have reported that the existing ICIs, such as PD-1 and CTLA-4 inhibitors, have limited efficacy in the treatment of AML patients. Moreover, it was reported that AML patients who received hematopoietic stem cell transplantation developed severe graft-versus-host disease after ICI treatment[34–35]. Thus, it is urgent to identify new immunotherapies for AML patients. To address this need, our model included five risk factors: CXCR2, LY6E, SUSD3, S1PR5, and IL1RAP. High expression of IL1RAP is independently associated with low overall survival in patients with AML, and IL1RAP may serve as a potential immunotherapeutic target for malignancies[36]. Meanwhile, CXCR2, produced by tumors, induces neutrophil extracellular traps that interfere with immune cytotoxicity[37]; thus, inhibiting CXCR2 may significantly reduce neutrophil infiltration and enhance anti-tumor T cell activity via promoting CD8 T cell activation[38]. Furthermore, the increased expression of LY6E was also correlated with the increased expression of immune checkpoint molecules, PD-L1 and CTLA-4, as well as the reduced activation of natural killer cells[39]. These findings suggest that CXCR2, LY6E, and IL1RAP may provide new targets for immunotherapy in AML patients.
Conversely, our model also included four protective factors: CD160, CLEC5A, NTNG2, and SLC4A1. CD160, which is highly expressed in natural killer cells, plays a significant role in killing tumor cells. Studies have shown that CD160 factors regulate functions of natural killer cells and control potential early cancer therapy[40]. CLEC5A is a cell surface receptor associated with the activation and differentiation of myeloid cells. In the current study, we report for the first time that NTNG2 is highly expressed in AML. Although NTNG2 is not well reported in AML, it is associated with the development of other cancers, such as breast and thyroid cancer[41].
In the current study, we used high-throughput transcriptomic data to comprehensively assess the effects of all CSMs on the prognosis of 130 AML patients from TCGA-LAML, and ultimately constructed a 9-CSMs prognostic model associated with AML patient outcomes. In clinical practice, the expression levels of each CSM may be obtained by some protein quantitative methods, such as cytometry by time-of-flight (CyTOF). Then, protein quantification data may be used to calculate the 9-CSMs risk score using the coefficients reflecting the correlation of each marker with disease prognosis derived from extensive clinical research and statistical analysis in the current study. The 9-CSMs score may help us stratify patients by risk, provide personalized treatment plans, aid in the development of targeted new drugs, and accelerate the process of bringing new therapies to market (Fig. 9).
Figure
9.
Clinical applications of the 9-CSMs risk score.
In clinical practice, the expression level of each CSM can be obtained by some protein quantitative methods, such as cytometry by time-of-flight (CyTOF). The protein quantitative data then be used to calculate a 9-CSMs risk score that uses coefficients reflecting the association of markers with widely clinically-derived disease outcomes. In conclusion, integrating 9-CSMs into prognostic models holds an important promise for advancing personalized medicine in oncology and improving patient outcomes.
However, the current study has several limitations. 1) We did not specifically identify the direct molecular relationships between CSMs and the progress in AML patients. The deep mechanisms of how CSMs reflect prognosis, such as the cell-cell interaction, biological processes in cells, and cellular signal transduction, remain unclear. 2) In addition, we trained and validated the prognostic model entirely based on public databases. We would focus on conducting a prospective study to validate our findings experimentally in the future. 3) Moreover, the risk definition of 9-CSMs risk score in individual patients is unclear, and we need to collect more quantitative protein data to evaluate the association between 9-CSMs risk score value and prognosis to determine the appropriate threshold for clinical application.
In summary, we identify a 9-CSMs prognostic model that may compute the risk score and predict AML prognosis accurately. The predictive value of the 9-CSMs prognostic model is validated at both transcriptomic and proteomic levels, establishing it as an independent prognostic factor for AML patients. Further analysis indicates that patients with high 9-CSMs risk scores are resistant to chemotherapy and that inhibitors of PI3K are potential therapeutic drugs for patients with high 9-CSMs scores. Overall, our research will provide valuable information for predicting prognosis and guiding therapeutic decisions for AML patients.
Fundings
The current study was supported by the National Natural Science Foundation of China (Grant Nos. 32200590 to K.L., 81972358 to Q.W., 91959113 to Q.W., and 82372897 to Q.W.) and the Natural Science Foundation of Jiangsu Province (Grant No. BK20210530 to K.L.).
De Kouchkovsky I, Abdul-Hay M. 'Acute myeloid leukemia: a comprehensive review and 2016 update'[J]. Blood Cancer J, 2016, 6(7): e441. doi: 10.1038/bcj.2016.50
[3]
Pollyea DA, Jordan CT. Therapeutic targeting of acute myeloid leukemia stem cells[J]. Blood, 2017, 129(12): 1627–1635. doi: 10.1182/blood-2016-10-696039
[4]
Vadakekolathu J, Minden MD, Hood T, et al. Immune landscapes predict chemotherapy resistance and immunotherapy response in acute myeloid leukemia[J]. Sci Transl Med, 2020, 12(546): eaaz0463. doi: 10.1126/scitranslmed.aaz0463
[5]
The Cancer Genome Atlas Research Network. Genomic and epigenomic landscapes of adult de novo acute myeloid leukemia[J]. N Engl J Med, 2013, 368(22): 2059–2074. doi: 10.1056/NEJMoa1301689
[6]
Zeng AGX, Bansal S, Jin L, et al. A cellular hierarchy framework for understanding heterogeneity and predicting drug response in acute myeloid leukemia[J]. Nat Med, 2022, 28(6): 1212–1223. doi: 10.1038/s41591-022-01819-x
[7]
Kantarjian H, Ravandi F, O'Brien S, et al. Intensive chemotherapy does not benefit most older patients (age 70 years or older) with acute myeloid leukemia[J]. Blood, 2010, 116(22): 4422–4429. doi: 10.1182/blood-2010-03-276485
[8]
Lowenberg B, Downing JR, Burnett A. Acute myeloid leukemia[J]. N Engl J Med, 1999, 341(14): 1051–1062. doi: 10.1056/NEJM199909303411407
[9]
Grob T, Sanders MA, Vonk CM, et al. Prognostic value of FLT3-internal tandem duplication residual disease in acute myeloid leukemia[J]. J Clin Oncol, 2023, 41(4): 756–765. doi: 10.1200/JCO.22.00715
[10]
Döhner H, Wei AH, Appelbaum FR, et al. Diagnosis and management of AML in adults: 2022 recommendations from an international expert panel on behalf of the ELN[J]. Blood, 2022, 140(12): 1345–1377. doi: 10.1182/blood.2022016867
[11]
Mercatelli D, Cabrelle C, Veltri P, et al. Detection of pan-cancer surface protein biomarkers via a network-based approach on transcriptomics data[J]. Brief Bioinform, 2022, 23(5): bbac400. doi: 10.1093/bib/bbac400
[12]
Andrews TE, Wang D, Harki DA. Cell surface markers of cancer stem cells: diagnostic macromolecules and targets for drug delivery[J]. Drug Deliv Transl Res, 2013, 3(2): 121–142. doi: 10.1007/s13346-012-0075-1
[13]
Hu Z, Yuan J, Long M, et al. The Cancer Surfaceome Atlas integrates genomic, functional and drug response data to identify actionable targets[J]. Nat Cancer, 2021, 2(12): 1406–1422. doi: 10.1038/s43018-021-00282-w
[14]
Sadovnik I, Herrmann H, Blatt K, et al. Evaluation of cell surface markers and targets in leukemic stem cells (LSC) reveals distinct expression profiles, unique drug effects, and specific checkpoint regulation in AML LSC and CML LSC[J]. Blood, 2016, 128(22): 4234. doi: 10.1182/blood.V128.22.4234.4234
[15]
Wu Z, Ou J, Liu N, et al. Upregulation of Tim-3 is associated with poor prognosis in acute myeloid leukemia[J]. Cancer Med, 2023, 12(7): 8956–8969. doi: 10.1002/cam4.5549
[16]
Heitmann JS, Hagelstein I, Hinterleitner C, et al. Identification of CD318 (CDCP1) as novel prognostic marker in AML[J]. Ann Hematol, 2020, 99(3): 477–486. doi: 10.1007/s00277-020-03907-9
[17]
Metzeler KH, Heilmeier B, Edmaier KE, et al. High expression of lymphoid enhancer-binding factor-1 (LEF1) is a novel favorable prognostic factor in cytogenetically normal acute myeloid leukemia[J]. Blood, 2012, 120(10): 2118–2126. doi: 10.1182/blood-2012-02-411827
[18]
Tibshirani R. Regression shrinkage and selection via the lasso[J]. J R Stat Soc Series B Stat Methodol, 1996, 58(1): 267–288. doi: 10.1111/j.2517-6161.1996.tb02080.x
[19]
Zou H. The adaptive lasso and its oracle properties[J]. J Am Stat Assoc, 2006, 101(476): 1418–1429. doi: 10.1198/016214506000000735
[20]
Tay JK, Narasimhan B, Hastie T. Elastic net regularization paths for all generalized linear models[J]. J Stat Softw, 2023, 106(1). doi: 10.18637/jss.v106.i01
[21]
Butler A, Hoffman P, Smibert P, et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species[J]. Nat Biotechnol, 2018, 36(5): 411–420. doi: 10.1038/nbt.4096
[22]
Li K, Du Y, Cai Y, et al. Single-cell analysis reveals the chemotherapy-induced cellular reprogramming and novel therapeutic targets in relapsed/refractory acute myeloid leukemia[J]. Leukemia, 2023, 37(2): 308–325. doi: 10.1038/s41375-022-01789-6
[23]
Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response[J]. Nat Med, 2018, 24(10): 1550–1558. doi: 10.1038/s41591-018-0136-1
[24]
Fu J, Li K, Zhang W, et al. Large-scale public data reuse to model immunotherapy response and resistance[J]. Genome Med, 2020, 12(1): 21. doi: 10.1186/s13073-020-0721-z
[25]
Yilmaz M, Kantarjian H, Ravandi F. Acute promyelocytic leukemia current treatment algorithms[J]. Blood Cancer J, 2021, 11(6): 123. doi: 10.1038/s41408-021-00514-3
[26]
Zhang X, Sun J, Yu W, et al. Current views on the genetic landscape and management of variant acute promyelocytic leukemia[J]. Biomark Res, 2021, 9(1): 33. doi: 10.1186/s40364-021-00284-x
[27]
Martelli AM, Nyåkern M, Tabellini G, et al. Phosphoinositide 3-kinase/Akt signaling pathway and its therapeutical implications for human acute myeloid leukemia[J]. Leukemia, 2006, 20(6): 911–928. doi: 10.1038/sj.leu.2404245
[28]
Stefanidakis M, Karjalainen K, Jaalouk DE, et al. Role of leukemia cell invadosome in extramedullary infiltration[J]. Blood, 2009, 114(14): 3008–3017. doi: 10.1182/blood-2008-04-148643
[29]
Fianchi L, Quattrone M, Criscuolo M, et al. Extramedullary involvement in acute myeloid leukemia. A single center ten years' experience[J]. Mediterr J Hematol Infect Dis, 2021, 13(1): e2021030. https://pubmed.ncbi.nlm.nih.gov/34007418/
[30]
Sun Y, Berleth N, Wu W, et al. Fin56-induced ferroptosis is supported by autophagy-mediated GPX4 degradation and functions synergistically with mTOR inhibition to kill bladder cancer cells[J]. Cell Death Dis, 2021, 12(11): 1028. doi: 10.1038/s41419-021-04306-2
[31]
Yu K, Shi C, Toral-Barza L, et al. Beyond rapalog therapy: preclinical pharmacology and antitumor activity of WYE-125132, an ATP-competitive and specific inhibitor of mTORC1 and mTORC2[J]. Cancer Res, 2010, 70(2): 621–631. doi: 10.1158/0008-5472.CAN-09-2340
[32]
Bei S, Li F, Li H, et al. Inhibition of gastric cancer cell growth by a PI3K-mTOR dual inhibitor GSK1059615[J]. Biochem Biophys Res Commun, 2019, 511(1): 13–20. doi: 10.1158/1535-7163.MCT-17-1178
[33]
Ramchandren R, Domingo-Domènech E, Rueda A, et al. Nivolumab for newly diagnosed advanced-stage classic hodgkin lymphoma: safety and efficacy in the phase II CheckMate 205 study[J]. J Clin Oncol, 2019, 37(23): 1997–2007. doi: 10.1200/JCO.19.00315
[34]
Haverkos BM, Abbott D, Hamadani M, et al. PD-1 blockade for relapsed lymphoma post-allogeneic hematopoietic cell transplant: high response rate but frequent GVHD[J]. Blood, 2017, 130(2): 221–228. doi: 10.1182/blood-2017-01-761346
[35]
Meric-Bernstam F, Larkin J, Tabernero J, et al. Enhancing anti-tumour efficacy with immunotherapy combinations[J]. Lancet, 2021, 397(10278): 1010–1022. doi: 10.1016/S0140-6736(20)32598-8
[36]
Barreyro L, Will B, Bartholdy B, et al. Overexpression of IL-1 receptor accessory protein in stem and progenitor cells and outcome correlation in AML and MDS[J]. Blood, 2012, 120(6): 1290–1298. doi: 10.1182/blood-2012-01-404699
[37]
Teijeira Á, Garasa S, Gato M, et al. CXCR1 and CXCR2 chemokine receptor agonists produced by tumors induce neutrophil extracellular traps that interfere with immune cytotoxicity[J]. Immunity, 2020, 52(5): 856–871. e8.
[38]
Cheng Y, Mo F, Li Q, et al. Targeting CXCR2 inhibits the progression of lung cancer and promotes therapeutic effect of cisplatin[J]. Mol Cancer, 2021, 20(1): 62. doi: 10.1186/s12943-021-01355-1
[39]
AlHossiny M, Luo L, Frazier WR, et al. Ly6E/K signaling to TGFβ promotes breast cancer progression, immune escape, and drug resistance[J]. Cancer Res, 2016, 76(11): 3376–3386. doi: 10.1158/0008-5472.CAN-15-2654
[40]
Tu TC, Brown NK, Kim TJ, et al. CD160 is essential for NK-mediated IFN-γ production[J]. J Exp Med, 2015, 212(3): 415–429. doi: 10.1084/jem.20131601
[41]
Hao W, Yu M, Lin J, et al. The pan-cancer landscape of netrin family reveals potential oncogenic biomarkers[J]. Sci Rep, 2020, 10(1): 5224. doi: 10.1038/s41598-020-62117-5
Authors and Reviewers
DownLoad:
