The association between the immune score and the response to imatinib therapy
We obtained gene expression profiles and clinical information of all patients from the GSE130404 dataset who had received an initial pathologic diagnosis. Out of these patients, 58 (60.4%) were female, and 38 (39.6%) were male. Among them, 83 patients responded to the therapy, while 13 patients did not. We included a total of 96 cases with complete gene expression data and clinical information in our analysis. By using the ESTIMATE algorithm, we calculated the immune score for each patient, which ranged from 1974.5 to 3828.2 (Fig. 1A). There was a significant difference in the immune score between the group of patients who responded to the therapy and the group who did not respond. Specifically, the immune score was higher in the responding group compared to the non-responding group (p = 0.0015). We also evaluated the predictive potential of the immune score using the area under the curve of the receiver operating characteristic (ROC), which was 0.753 (Fig. 1B), suggesting a moderate potential as a predictive biomarker.
Identification of immune-related genes (IRGs)
To identify immune-related genes (IRGs), we performed differential expression analysis using the Linear Models for Microarray Data (LIMMA) method, using the median value as the cutoff point. This analysis revealed 428 differentially expressed genes (DEGs), with 134 down-regulated and 294 up-regulated genes, in the comparison between samples with higher immune scores and samples with lower immune scores in CML patients. The volcano plot of these differentially expressed IRGs is shown in Fig. 2A, and their expression patterns are depicted in the heat map in Fig. 2B.
To gain further insights into the potential functions of these IRGs, we conducted pathway enrichment analysis using the clusterProfiler package. The KEGG analysis showed that these IRGs were primarily enriched in pathways such as Hematopoietic cell lineage, Th17 cell differentiation, Th1 and Th2 cell differentiation, Cell cycle, T cell receptor signaling pathway, and Cytokine-cytokine receptor interaction (Fig. 2C). The top gene ontology (GO) terms related to biological processes (BP) included mononuclear cell differentiation, regulation of T cell activation, lymphocyte differentiation, leukocyte-mediated immunity, and leukocyte cell-cell adhesion. The top GO terms related to cellular components (CC) included the external side of the plasma membrane, spindle, and membrane raft. The top GO terms related to molecular functions (MF) included MHC protein complex binding, immune receptor activity, and cytokine binding (Fig. 2D).
Construction of the weighted gene-co-expression network
We constructed a gene co-expression network using the weighted gene co-expression method. To determine the appropriate soft-thresholding power for the network construction, we calculated the power value by predicting the scale-free topology fit index (scale-free R2 = 0.90) for different powers. Among the tested powers, we chose β = 17 as the soft-thresholding power to ensure a scale-free network, as indicated by the analysis of the mean connectivity for different powers (Fig. 3A). By setting the cut line at 0.25, we merged similar modules, resulting in the formation of several modules including grey, green, brown, yellow, blue, and turquoise (Fig. 3B).
Identification of clinically significant modules
To examine the relationship between the phenotype and module eigengenes, we analyzed the connectivity of eigengenes (Fig. 4A). We calculated the module-trait relationship using the Pearson correlation coefficient and visualized it in a heat map. In the module-trait relationships, we observed that the turquoise module showed a stronger positive correlation with response to imatinib (r = 0.27, p = 0.008) and immune score (r = 0.84, p = 3e-26) compared to other modules (Fig. 4B). Similarly, the grey module exhibited a positive correlation with response to imatinib (r = 0.24, p = 0.02) and immune score (r = 0.83, p = 9e-26).
Functional annotation of modules and identification of bub genes
We used the ggraph and igraph R packages to visualize different modules among 389 immune-related genes (IRGs) (Fig. 5A). To compare the functional profiles of different modules, we performed KEGG pathway analysis using the IRGs within each module, utilizing the compareCluster function in the clusterProfiler package. The blue module was significantly enriched in the Cell cycle pathway. The grey module was enriched in the Toll-like receptor signaling pathway, Chemokine signaling pathway, and Cytokine-cytokine receptor interaction. The turquoise module was mainly enriched in the T cell receptor signaling pathway (Fig. 5B). We analyzed the hub genes based on the node degree values in the network and identified the top 10 genes with the highest scores as hub genes (Fig. 5C). Among these genes, PLEKHA1, PRKCH, CD247, BCL11B, CD3D, KLF12, CD7, SBK1, ITPR3, and IL2RB1 showed lower gene expression in non-responders compared to the responding group, suggesting their potential involvement in immune-related molecular mechanisms associated with the response to imatinib (Fig. 5D).
Validating biomarkers of response to imatinib therapy
To develop a practical and cost-effective predictive assay for customizing therapy in newly diagnosed chronic myeloid leukemia (CML) patients, we narrowed down our candidate IRGs by selecting genes with an AUC value greater than 0.7. Among them, three genes were validated in an independent cohort of 59 CML patients who exhibited a major cytogenetic response after 12 months of imatinib therapy. In the GSE130404 dataset, the AUC values for IL10RA, SCN9A, and SLC26A11 were 0.722, 0.728, and 0.778, respectively (Fig. 6C, n = 96, p < 0.05). In the validating datasets, the AUC values were 0.732, 0.710, and 0.741, respectively (Fig. 6D). IL10RA and SLC26A11 were up-regulated in non-responders compared to responders, while SCN9A was down-regulated in non-responders (Fig. 6A and B). Additionally, the immune score exhibited a similar expression difference between non-responders and responders in both datasets (Fig. 6). These findings suggest that the immune score and the three identified genes may serve as biomarkers for predicting the response to imatinib therapy.