An inflammatory response-related gene signature associated with immune status and prognosis of acute myeloid leukemia

Data collection and processing

The transcriptomic dataset and associated clinical information of 151 AML patients were acquired from the TCGA web site (TCGA-AML Cohort, https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), and 264 normal samples were acquired from GTEx web site (GTEx-Normal Cohort, https://xena.ucsc.edu/). The transcriptomic dataset and clinical data of 356 AML specimens were acquired from the ICGC website (LAML-US Cohort, https://dcc.icgc.org/). The patients whose data were included in these sets were primarily American citizens, and all were diagnosed with AML. The information from TCGA, GTEx, and ICGC is publicly available in accordance with their data acquisition and publication protocols.

Establishment and verification of a prognostic inflammation reaction-associated genetic signature

First, as presented in Supplementary Table 1, we acquired 200 inflammation reaction-related genes (IRRGs) from the Molecule Signature data base. Then, per the standards of fold change > 2 and a false discovery rate (FDR) < 0.05, we determined the DEGs between AML specimens and healthy specimens within the TCGA cohort with “limma” R software. Least absolute shrinkage and selection operator (LASSO)-penalized Cox regressive analysis was used to establish a prognosis model and reduce the risk of overfitting to a minimum [15]. The “glmnet” R package was used to select and shrink variates, equaling certain regression coefficients to 0 and yielding a rational model.

Patient risk scores were computed based on the expression level of all evaluated inflammatory reaction-associated genes and their relevant regression coefficients. The calculation was carried out according to the following equation: risk score= e^{sum (every genetic expression × relevant coefficient)}. Then, patients were separated into risk_high and risk_low groups based on the mid-value of the risk score. PCA analysis was further employed to study the distributional status of prognostic genes in various groups with “Rtsne” and “ggplot2” R packages.

To study the overall survival (OS) of the risk_high group and risk_low group, survival analysis was performed using the “survminer” R package. Moreover, the “survival” and “timeROC” R programs were implemented to complete time-reliant receiver operating characteristic (ROC) curve analyses to assess the prediction potential of the prognosis hallmark.

Functional enrichment analysis of inflammatory reaction-related DEGs

To reveal the potential roles of the prognostic signature, we employed functional enrichment analyses, like gene ontology, Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene set enrichment analysis (GSEA) with program 4.1 according to the DEGs between the risk_high group and risk_low group [16]. The Banjamini-Hochberg (BH) adjusted P < 0.05 was considered statistically significant.

Tumor microenvironment (TME) and immune response analysis

Immunity and stroma score was carried out to measure the infiltrative status of immunocytes and stromal cells across the two groups [17]. Spearman’s correlation was employed to examine the association of risk score with immune and stromal score. Moreover, tumor stem cell-like characteristics, including RNA stemness score (RNAss) and DNA stemness score (DNAss), were analyzed using the PanCancer TCGA database of tumor stem cell transcriptome and epigenetic data [18]. The stemness indices of RNA and DNA were calculated according to mRNA expression and DNA methylation, respectively [19]. The correlation between cancer stemness and risk score was identified using Spearman correlative analyses.

Analysis of chemosensitivity

The NCI-60 database covers 60 tumor lineage cell lines from nine different cancer types hosted on the Cell Miner website (https://discover.nci.nih.gov/cellminer). A Pearson correlative assay was performed to explore the relationship between genetic expression prognosis and drug susceptibility. The association assay was performed on 216 medicines, selected from medications with FDA approval or those undergoing clinical trial testing (Supplementary Table 2).

Validation of prognosis gene expression between AML and control samples by qRT-PCR

Ten paired AML tissue and corresponding normal peripheral blood samples were acquired from Qinghai Provincial People’s Hospital. The mRNA expression levels of 11 prognosis genes were measured in all samples using qRT-PCR. Based on the manufacturer’s specifications, the target RNAs were extracted using TRIzol reagent (Servicebio). Extracted RNAs were then converted into cDNA by reverse transcription using the Revert Aid First Strand cDNA Synthetic Tool (Thermo). Fast Start Universal SYBR Green Master Mix (Roche) was used to complete the qRT-PCR assay using StepOne (Applied Biosystems). The relevant primer sequences are listed in Supplementary Table 3. Each RNA sampling procedure was repeated three times. The comparative expression levels of IRRGs were computed using the 2^-ΔΔCt approach to compare expression levels across samples.

Validation of prognostic gene expression between AML and control samples by immunofluorescence (IF)

Paraffin-embedded tumor tissues were deparaffinized for antigen repair. IF staining was performed using anti-ACVR2A, CCL22, EBI3, EDN1, FFAR2, HRH1, ICOSLG, IL-10, INHBA, ITGB3, and LAMP3 antibodies, according to the standard protocol. The nuclei were then stained with DAPI (Invitrogen), and imaging results were observed and analyzed via fluorescence microscopy.

Statistics

The Wilcoxon test was used to compare DEGs between AML and control samples, and the chi-square test was used to compare the percentages. The Mann-Whitney test was employed to contrast the ssGSEA score of immunocytes or immunity paths between the risk_high and risk_low groups, and the BH approach was utilized to modify the P result. K-M analyses were used to compare the OS among the diverse groups. The association between prognosis model risk score, prognosis genetic expression, stemness and stroma score, immunity score, and medicine susceptibility was further determined via Spearman and Pearson correlative analyses. To generate plots, R program 3.6.3, including the various packages of venn, igraph, ggplot2, pheatmap, ggpubr, corrplot, and survminer were employed. Finally, a two-tailed P value < 0.05 was considered statistically significant.

Determination of prognostic inflammation-associated DEGs in the TCGA cohort

A flow diagram of this study is presented in Figure 1. The research subjects comprised 151 AML patients from the TCGA-AML cohort, 264 normal samples from the GTEx cohort, and 356 AML patients from the ICGC (LAML-US) cohort. In AML and neighboring normal samples, 79 DEGs associated with inflammatory response and 49 prognostic genes were identified (Figure 2A). The heatmap shows the expression levels of the 21 overlapping genes in AML and normal samples (Figure 2B). Univariable Cox analysis revealed an association between the expression levels of 49 prognosis genes and OS (Supplementary Figure 1A). The 21 inflammatory response-related genes were correlated with OS and were expected to be potential prognostic markers. The hazard rate of CCL22 gene was 18.178 (95% CI = 4.296-76.910, P < 0.001, Figure 2C). As shown in Figure 2D, we determined the correlative network of these genes.

Establishment of a prognostic model for the TCGA cohort

LASSO-Cox regressive analysis was performed to construct the prognostic model of the aforementioned 21 genes, which were identified based on the optimum value of λ (Supplementary Figure 1B, 1C). The patients’ risk scores were assessed as follows: risk score = 0.120 * expression level of ACVR2A + 1.530 * expression level of CCL22 + 0.168 * expression level of EBI3 + 0.691 * expression level of EDN1 + 0.066 * expression level of FFAR2 + 0.661 * expression level of HRH1 + 0.313 * expression level of ICOSLG + 0.033 * expression level of IL10 + 0.095 * expression level of INHBA + 0.019 * expression level of ITGB3 + 1.103 * expression level of LAMP3.

Subsequently, patients were separated into risk_high and risk_low groups based on the mid-value cut-off (Figure 3A). Early mortality was observed more frequently for Risk_high patients than for risk_low patients (Figure 3B). Furthermore, PCA analyses revealed that patients in the risk_high and risk_low groups were mainly distributed in two orientations (Figure 3C). Similarly, the K-M curve revealed that risk_high patients displayed a remarkably inferior OS compared to risk_low patients (Figure 3D, P < 0.001). ROC curves were used to forecast OS in the constructed prognostic model, and the AUC was 0.799 at 1 year, 0.803 at 2 years, and 0.780 at 3 years (Figure 3E). To reveal the association between the 11 prognosis genes and prognoses, we performed survival analyses according to the optimum threshold of each prognosis gene. The results revealed that for all genes, high expression levels were significantly related to worse OS (Supplementary Figure 2A-K, P < 0.001). As presented in Supplementary Figure 3, the expression levels of each of the 11 prognostic genes were higher in AML samples than in the normal controls, which is consistent with the results shown in Figure 2B.

Validation of gene signatures in the ICGC cohort

AML patients within the ICGC cohort were separated into risk_high and risk_low groups according to the mid-value of the TCGA cohort to determine the validity of the established prognosis model (Figure 3F). However, all patients with survival data in the ICGC cohort died, so the data were unsuitable for comparing mortality in the different risk groups (Figure 3G). As per PCA analysis, the results verified the scattered distribution of AML patients in risk_high and risk_low groups, showing similar results to those obtained from the TCGA cohort (Figure 3H). As shown in Figure 3I, patients in the risk_high group exhibited a lower survival rate than those in the risk_low group. In addition, the AUC of the 11-gene hallmark reached 0.738 at 1 year, 0.745 at 2 years, and 0.789 at 3 years, indicating that the risk score yielded satisfactory prognostic predictability (Figure 3J).

Biological function and pathway analyses

The different OS rates of patients in the risk_high group and risk_low group were studied, and GSEA was further employed to explore the underlying diversity of biological functions and pathways between them. GO term analysis revealed that the biological process of T cell proliferation was remarkably enriched in risk_high patients (Figure 4A; Supplementary Figure 4). In addition, five KEGG pathways were mainly sponged in the risk_high patients, with a FDR < 0.05, which is presented in Figure 4B and Supplementary Figure 5. The pathways of the KEGG database demonstrated that the remarkably expressed DEGs in the risk_high patients were predominantly associated with the pathways relevant to the inflammatory response, including B cell receptor, chemokine, and Fc-g receptor-mediated phagocytosis. These data indicate that these DEGs may be highly influential on the immunoresponse.

Potential relevance of gene signature in TIME

The TME immune status is mainly associated with cancer development, drug tolerance, and clinical results. To evaluate the association between risk score and immune status, ssGSEA was performed to study the enrichment score of diverse immunocyte subset-related roles and pathways. All genes relevant to ssGSEA analysis are presented in Supplementary Table 4. We discovered that the process of antigenic presentation within the TCGA cohort, such as iDCs, pDCs, APC co-suppression, APC co-activation, HLA, and MHC class I, was significantly promoted in risk_high patients with a modified P < 0.05 (Figure 5A, 5B). Moreover, the ratio of Tfh cells, Th1 cells, Treg cells, T cell co-activation, and T cell co-suppression was higher in risk_high patients than in risk_low patients, suggesting that differences in T cell regulation may play a potential role in observed differences in risk_high and risk_low groups. In addition, our team observed that CCR, checkpoint, macrophage, neutrophilic cell, inflammation-facilitating activity, and type II IFN reaction activity scores were also higher in risk_high patients than in risk_low patients with a modified P < 0.05 (Figure 5A, 5B). The comparisons of the above indicators within the ICGC cohort between the two risk groups resembled that within the TCGA cohort with a modified P < 0.05 (Figure 5C, 5D).

To further study whether risk score was associated with cancer stemness and TIME, RNAss and DNAss were performed to measure tumor stemness. Stroma and immunity score were also employed to assess TIME. As presented in Figure 5E, our team discovered that the risk score was not remarkably associated with DNAss but remarkably positively correlated with stromal and immune scores (P < 0.001) and negatively correlated with RNAss (P < 0.001).

In addition, PD-1/PD-L1 and PD-1/PD-L2 paths are critical modulators of tumor immunoescape. Thus, testing the expression levels of immunity checkpoints like PD-L1 and PD-L2 is a vital approach for immunological therapy. In comparison with the risk_low patients, the expression levels of PD-L1 and PD-L2 were higher in the risk_high patients (Figure 6A, 6B), and the relationship between the expression levels of PD-L1 and PD-L2 and the risk score was positive (Figure 6C, 6D).

Prognostic gene expression and oncocyte susceptibility to chemotherapy

The multidrug resistance-associated protein (MRP) family is essential for modulating drug resistance-associated tumor genes. We found that patients in risk_high group exhibited higher expression levels of MRP2 and MRP3 than risk_low patients (Figure 6E, 6F). Additionally, the expression levels of MRP2 and MRP3 were positively related to the risk score (Figure 6G, 6H).

Subsequently, in NCI-60 lineage cells, we measured the expression of prognosis genes and investigated the association between their expression levels and medicine susceptibility (Supplementary Table 5). The results suggested that the prognostic genes were related to certain chemotherapeutic drug sensitivities (P < 0.01). Representative results are as follows: elevated expression levels of INHBA, HRH1, and IL-10 were associated with elevated susceptibility of oncocytes to dasatinib, midostaurin, lomustine, and azacitidine. In contrast, the increased expression levels of CCL22 and ACVR2A were associated with decreased susceptibility of oncocytes to midostaurin and etoposide (Figure 7). Notably, the elevated expression of both INHBA and HRH1 was related to increased drug susceptibility of oncocytes to dasatinib, which was accepted by the FDA as a first-line therapy for AML in 2006.

Verification of prognosis gene expression between AML and control samples by qRT-PCR and IF

To validate the diverse expression of the 11 prognostic genes (ACVR2A, CCL22, EBI3, EDN1, FFAR2, HRH1, ICOSLG, IL-10, INHBA, ITGB3, and LAMP3) between AML and corresponding healthy samples, the expression level of mRNA and protein was analyzed by qRT-PCR and IF. As expected, the prognosis gene expression levels were higher in AML specimens than in healthy specimens based on the qRT-PCR assay (Figure 8, P < 0.001). As depicted in Figure 9, the results of IF were similar to those of qRT-PCR (P < 0.01). These experimental results further confirmed the mRNA levels of 11 prognostic genes in the TCGA database (Supplementary Figure 3).