Gene expression and relevant clinicopathological information on AML were from databases named The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). One GEO AML cohort named GSE10358 and the TCGA-AML cohort were obtained for further relative analyses. The fragments per kilobase million (FPKM) values from the TCGA-AML set were converted to transcripts per kilobase million (TPM) and were then assumed to be the same as those from microarrays (Conesa et al., 2016; Zhao et al., 2021). The two chosen AML cohorts were combined. We removed data from individuals for whom we had no information on their overall survival; hence, 242 AML patients were included in our analysis. The clinical factors shared by the two AML groups were gender, age, duration of follow-up, and survival status. To confirm our proposed prognostic risk model, we also gathered expression data of pediatric AML samples from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database (https://ocg.cancer.gov/programs/target) as an external validation cohort. Whole blood cohorts from GTEx downloaded from the UCSC Xena database (https://xenabrowser.net/datapages/) served as control samples for the analyses on the TCGA-AML and TARGET-AML samples.

A total of 60 CMRGs were retrieved from previous publications (Zhang et al., 2020), as shown in Supplementary Table S1. We used the package named “ConsensusClusterPlus” in R software to perform consensus clustering analysis. The criteria were as follows: Firstly, the cumulative distribution function (CDF) curve should increase gradually and smoothly. Then, the size of all groups was large enough. Finally, the intra-group correlation should increase after clustering, while the inter-group correlation decreased. Also, gene set variation analysis (GSVA) was carried out using the hallmark gene set (c2. cp.kegg.v7.2) to investigate differences in chosen CMRGs in the biological processes.

We evaluated the relationships between clustered costimulatory molecule subtypes, clinical features, and outcomes. The characteristics mainly included gender and age. In addition, Kaplan-Meier curves, which were created using the “survival” and “survminer” packages, were used to evaluate the differences in overall survival across various costimulatory molecule subtypes.

We assessed the immune, estimate and stromal scores of each AML sample in this study by using the ESTIMATE algorithm. Using the CIBERSORT algorithm, the fractions of 22 immune cell subgroups of each AML sample were estimated (Newman et al., 2015). The immune cell infiltration in the TME of AML was also identified by using an algorithm named single-sample gene set enrichment analysis (ssGSEA) (Rooney et al., 2015).

Differentially expressed genes between the two different costimulatory molecule subtypes were identified in R software. To explore the probable activities of costimulatory molecule profile-related differentially expressed genes and determine associated functions, we used “clusterprofiler” package to carry out functional enrichment analyses on the chosen gene.

The CMRG-related risk model was generated, and the score was estimated for each AML sample to quantify the costimulatory molecule patterns. Differentially expressed genes were subjected to univariate Cox analysis. Using an unsupervised clustering technique, AML patients were divided into three distinct CMRG-related gene subtypes (gene subtypes A, B, and C) based on the findings of the univariate Cox analysis. The chosen 242 patients with AML were then randomly categorized into two sets at a ratio of 1:1, a training AML set (n = 121) and a testing AML set (n = 121). Based on CMRG-related genes with prognostic value, we carried out the Lasso Cox algorithm to minimize the over-fitting risk in the training AML set. We next analyzed each independent variable’s change and established a risk model by using 10-fold cross-validation. Candidate genes were further analyzed and chosen in the training AML set based on the using multivariate Cox regression results. The costimulatory molecule signature was calculated using the data of each gene’s coefficient and expression. Then, based on the median score of whole samples in the AML training set, samples of two sets were respectively divided into low- or high-risk groups. Survival analysis and the creation of receiver operating characteristic (ROC) curves were performed on AML patients belonging to two risk categories.

To explore differences in the chemotherapeutic drug curative effect in AML patients between two risk groups from the whole set, we calculated the values of semi-inhibitory concentration (IC50) of drugs using the “pRRophetic” package in R software.

All statistics were analyzed in R software (version 4.1.0). Statistical significance was all set at p < 0.05.

This research contained 60 CMRGs in total. We firstly investigated the somatic copy number variations in the 60 CMRGs. Figure 1A showed the locations of the copy number variation (CNV) alterations in the CMRGs on their respective chromosomes. Among all the CMRGs, ICOSLG, TNFRSF14, TNFRSF4, TNFRSF18, TNFRSF8, TNFRSF1B, RELT, TMIGD2, CD70, TNFSF14, and TNFRSF6B had widespread CNV increases, while TNFRSF25, TNFRSF9, TNFRSF11B, VTCN1, TNFRSF10D, TNFRSF13B, CD40, TNFRSF13C, and EDA showed CNV decreases (Figure 1B). A comparison on the mRNA levels of CMRGs between AML and normal samples revealed that the majority of CMRGs were positively linked with CNV alteration. Some CMRGs with CNV gain, including TNFRSF18, TNFRSF1B, RELT, TMIGD2, and CD70, were significantly elevated in samples from the AML cohort. Meanwhile, CMRGs with CNV loss, including TNFRSF25, TNFRSF11B, VTCN1, CD40, and EDA, were expressed at lower levels in patients with AML when compared to normal samples (Figure 1C), suggesting CNV might participate in regulating the CMRGs’ expression. However, several CMRGs with a high frequency of CNV gain or loss did not vary between AML and normal samples. Thus, CNV is not the sole factor involved in CMRGs’ expression regulation. Both the genetic landscape and expression levels of CMRGs were significantly different between AML samples and controls, showing that CMRGs play a latent role in the oncogenesis of AML.

A total of 455 patients from two AML cohorts (TCGA-AML and GSE10358) were integrated into our study at the beginning. Since only 242 AML samples had information on clinical results, we extracted their data of them for further analysis. The univariate Cox regression and Kaplan-Meier analyses revealed that 44 CMRGs were found might be related to the prognosis of AML (Supplementary Figure S1). A costimulatory molecule network illustrated the full picture of CMRG connections, regulator linkages, and their predictive relevance in AML patients (Figure 2A). To further investigate the expression features of CMRGs in AML, we categorized the AML samples based on the expression profiles of 60 CMRGs using a consensus clustering technique. We found that k = 2 seemed to be the ideal choice for classifying the complete cohort into two categories (Figure 2B). Next, the results of PCA analysis revealed the obvious costimulatory molecule differences between the two costimulatory molecule subtypes (Figure 2C). Kaplan-Meier curves then showed that subtype A patients had longer overall survival than subtype B patients (Figure 2D). Furthermore, we made a comparison of the clinical features between the two various subtypes of AML. As demonstrated in Figure 2E, cluster A was more strongly associated with a younger age than cluster B.

We found that subtype B was enriched in activated immune pathways, such as cytokine receptor interaction, cell adhesion molecules, graft versus host disease, antigen processing and presentation, B cell receptor signaling pathway, allograft rejection, Toll-like receptor signaling pathway, regulation of actin cytoskeleton, and chemokine signaling pathway (Figure 3A). To investigate the functional role of CMRGs in the TME of AML, we used the CIBERSORT algorithm to analyze the correlations between the costimulatory molecule subtype and the 22 immune cell subsets of each AML sample. As shown in Figure 3B, there had been significant differences in most immune cell infiltration between subtype A and subtype B. The infiltration levels of activated CD8+ T cells, activated B cells, activated CD4+ T cells, natural killer cells, immature B cells, activated dendritic cells, monocytes, macrophages, mast cells, neutrophilia, regulatory T cells, follicular helper T cells, type 1 helper T cells, type 2 helper T cells and type 17 helper T cells were higher in patients in subtype B than in those in subtype A.

The biological role of costimulatory molecule patterns in AML remains unknown. We used the “limma” package in R software to identify costimulatory molecule subtype-related differential expressed genes and then carried out the functional enrichment analysis based on them (Figures 4A,B). KEGG analysis indicated these identified costimulatory molecule subtype-related genes enriched in pathways related to cancers and immunology, suggesting that costimulatory molecules play vital roles in TME regulation (Figure 4A). The GO chord in Figure 4B also showed the GO mainly enriched terms and the significantly involved differential expressed genes. We then carried out the univariate Cox regression analysis on the subtype-related genes to screen out 501 genes that were closely related to overall survival time (Supplementary Table S2). To further explore potential regulation mechanisms, we performed a consensus clustering algorithm, which divided all AML patients into three subtypes named gene subtypes A-C. Kaplan-Meier curve in Figure 4C showed that patients in the subtype C group had the best favorable clinical outcome, whereas patients in subtype B showed the worst overall survival. The costimulatory molecule subtypes showed significant differences in the expression levels of CMRGs, which was consistent with the results of the costimulatory molecule patterns (Figure 4D). Moreover, costimulatory molecule gene subtype B patterns seemed more associated with elder age (Figure 4E).

The CMRG-related prognostic signature was then established. First, we randomly divided the AML patients into two groups, including a training group (n = 121) and a testing group (n = 121) at a ratio of 1:1. Following, LASSO and multivariate Cox analyses for costimulatory molecule subtype-related differentially expressed genes were carried out to identify CMRG-related genes with prognostic value in AML patients. According to the minimum partial probability of deviance (Snyder et al., 2014), hub overall survival-related genes remained after LASSO regression analysis (Figure 5A). We then carried out a multivariate Cox analysis to finally obtain 15 ones (GPR18, LGALS1, AOAH, DNMT3B, CBR1, ANKRD55, SIRPB2, DPY19L2, IL1R2, ST8SIA4, DOC2A, SERPINI2, GZMB, TNNT1, and SORCS2), which included nine high-risk genes and six low-risk genes (Supplementary Table S3). Hence, the CMRG-related risk model was constructed as follows: Risk score = -0.4339 * GPR18 + 0.3393 * LGALS1 - 0.3897 * AOAH + 0.3038 * DNMT3B + 0.2168 * CBR1 + 0.3566 * ANKRD55–0.3112 * SIRPB2 - 0.2369 * DPY19L2 + 0.1176 * IL1R2 -0.3160 * ST8SIA4 -0.2405 * DOC2A + 0.4186 * SERPINI2 + 0.3901 * GZMB + 0.1895 * TNNT1 + 0.3065 * SORCS2. The score of each AML patient was measured according to the signature, the median of which was chosen as the cutoff to divide AML patients into two risk score groups. The distribution of AML patients in the two costimulatory molecule subtypes, three gene subtypes, two risk score groups, and clinical outcomes were shown in Figure 5B. We then analyzed the risk scores and discovered a statistically significant difference among the three CMRG-related gene subtypes. As shown in Figure 5C, risk scores in the gene subtype B group were the highest, while the risk scores in the gene subtype C group were the lowest, which indicated that lower costimulatory molecule score might be related to immune activation. It was noteworthy that gene subtype B had a higher risk score than gene subtype A. Furthermore, Figure 5D showed the distribution of scores in two costimulatory molecule subtypes, where AML patients in the costimulatory molecule subtype B group had significantly higher scores than others in the costimulatory molecule subtype A group.

To validate the predictive value of the CMRG-related risk model in patients with AML, we separately divided AML patients in each set into two different risk groups according to the cutoff, which was the median value of those in the AML training set. In detail, AML patients whose scores were larger than the chosen cutoff were grouped as low risk, while others whose scores were smaller than the chosen cutoff were grouped as high risk. The survival curves revealed that AML patients in the training set with lower scores significantly had a favorable overall survival time when they were compared to those with higher scores (Figure 6A). In the training set, the 1-, 3-, and 5-year survival rates for AML patients were respectively 0.913, 0.94, and 0.978, which were represented by AUC values (Figure 6B). The heatmap showed the 15 CMRG-related genes’ expression in two risk groups in the AML training set (Figure 6C). The distribution of the costimulatory molecule risk score revealed that the overall survival time of AML patients in the training set decreased with the scores increase (Figure 6D). Patients in the AML testing set were similarly divided into two risk categories using the same algorithm and cutoff as the training set. Survival analysis revealed that the group with lower scores had considerably better clinical results (Figure 6E). While studies of 1-, 3-, and 5-year prognosis prediction efficiencies revealed that the CMRG-related risk score maintained high AUC values, this was not the case for the 5-year prognostic efficiency (Figure 6F). The heatmap and survival status of AML patients showing the variation tendencies of two risk groups from the testing set were shown in Figures 6G,H, respectively, further confirming the costimulatory molecule signature had an excellent ability to predict the clinical outcome of patients with AML.

To investigate the difference in TME and expression of checkpoints between two risk groups divided by constructed signature, we divided all the 242 AML patients into two groups according to cutoff and perform the CIBERSORT algorithm, by which we could assess the correlations between risk scores and immune cell abundance. As the scatter diagrams in Figure 7A showed the constructed costimulatory molecule risk score was positively correlated with monocytes and was negatively correlated with plasma cells, resting memory CD4+ T cells, resting mast cells, memory B cells, and gamma delta T cells (Figures 7A–F). Higher risk scores were also closely linked with higher immune scores as well as higher estimate scores (Figure 7G). Figure 7H shows that a total of 25 immune checkpoints, including LAG3, IDO1, and PDCD1, were different between the two groups from the total AML set. Moreover, the majority of immune cells were linked with the genes in the proposed model (Figure 7I).

We then selected chemotherapy drugs to evaluate the treatment sensitivities of patients in different risk AML groups. Interestingly, we found that the patients with higher costimulatory molecule scores had lower IC50 values for WZ-1-84, WO2009093972, SL 0101–1, S-trityl-L-cysteine, Roscovitine, Rapamycin, Parthenolide, NVP-TAE684, Kin001-135, GNF-2, Dasatinib, CGP-082996, CGP-60474, Bortezomib, and AZ628, while IC50 values of chemotherapeutics such as VX-702, Vorinostat, Thapsigargin, SB 216763, Plx4720, PF-562271, OSI-906, MK-2206, Mitomycin-C, Midostaurin, Gemcitabine, Embelin, Cytarabine, BX-795, Bosutinib, BI-D1870, AZD7762, Axitinib, AP-24534, AKT inhibitor III, AG014699 and ABT263 were significantly lower in the patients with lower costimulatory molecule scores (Supplementary Figure S2). Moreover, to validate the prognostic performance of the costimulatory molecule signature in pediatric AML, we extracted the data of samples from the TARGET database. It is obvious that most CMRG-related genes in the proposed signature were differentially expressed between pediatric AML samples and control samples (Supplementary Figure S3), suggesting these CMRG-related genes also played important roles in pediatric AML.