Immune regulatory genes impact the hot/cold tumor microenvironment, affecting cancer treatment and patient outcomes

The transcripts per million expression data of 33 cancer types were obtained from the UCSC database and normalized. Only data of cancer samples were used in this investigation, and all normal data were excluded. Clinical metadata, mutation annotation data, and copy number variation data were also obtained from the UCSC database. Single-cell RNA sequencing (scRNA-seq) datasets of PAAD, including CRA001160 (23), GSE111672↗ (24), GSE141017↗ (25), GSE148673↗ (26), GSE154778↗ (27), GSE158356↗ (28), GSE162708↗ (29), and GSE165399↗ (30), were downloaded from the Gene Expression Omnibus (GEO: https://www.ncbi.nlm.nih.gov/geo/↗) database and to identify the gene expression in various cell types. Immunohistochemistry (IHC) data were further downloaded from the Human Protein Atlas (HPA, http://www.proteinatlas.org↗) database to verify protein expression in PAAD.

The CIBERSORT algorithm (31) of the IOBR packages (v.0.99.9) (32) was employed to measure the infiltration of 22 immune cell types. Using the R package “ConsensusClusterPlus” (v 1.64.0) (33), we identified various clusters that differ in terms of immune infiltration by consensus clustering for each cancer type. The number of clusters was determined by the k value and the area under the cumulative distribution function curve. To ensure the accuracy of our classification results, we repeated this step 1,000 times. The Single-sample Gene Set Enrichment Analysis (ssGSEA), as implemented in the R package GSVA (version 0.99.9) (34), was utilized to quantify 13 immune function scores derived from the work of He et al. (35), including T cell proliferation (36) and MDSCs (37). Based on the infiltration levels of CD8+ T cells, activated NK cells, and M2-type macrophages, as well as the scores for cytolytic activity and T cell proliferation, we categorized the clusters into two distinct groups: “hot-immune” and “cold-immune.” Additionally, we utilized several computational algorithms, including the Tumor Immune Estimation Resource (TIMER) (38), EPIC (39), Microenvironment Cell Populations-counter (MCP-counter) (40), xCELL (41), and quanTIseq (42), to quantify immune cell infiltration and to identify key immune cell types associated with hot and cold tumor phenotypes.

Based on the identified cold and hot tumor types, we performed a univariate Cox regression model in the survival (v 3.2-7) package (43) to analyze the prognostic relationship between cold and hot tumors. Tumors with p-values <0.05 were retained for further analysis.

We sourced immune regulatory genes, including chemokines, receptors, major histocompatibility complex (MHC) genes, immunoinhibitors, and immunostimulators, from the Sangerbox database (44). Additionally, we obtained immune checkpoint genes, encompassing both inhibitory and stimulatory genes, based on the study by Thorsson et al. (45). We conducted an analysis of the correlation between immune regulatory and checkpoint genes and the infiltration of 22 immune cell types, as estimated by the CIBERSORT algorithm, across various cancer types.

Gene expression enrichment was evaluated through GSVA analysis under both unsupervised and parameter-free conditions (34). GSVA was employed to examine the various Kyoto Encyclopedia of Genes and Genomes (KEGG) and hallmark pathways between hot and cold tumor types. Downloads for KEGG and hallmark gene sets were obtained from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp↗), and hypergeometric p-values were adjusted using Benjamini–Hochberg multiple testing correction.

The GSEA algorithm was applied to identify expression profiles that could either activate or suppress hallmark pathways between high- and low-survival groups. After 100 permutations, an enriched gene set was obtained based on a p-value <0.05 and a false discovery rate of 0.25.

Data on drug sensitivity were obtained from The Genomics of Drug Sensitivity in Cancer database (https://www.cancerrxgene.org/↗). The half maximal inhibitory concentration (IC₅₀) values of each drug were downloaded using the R package “oncoPredict” (v 0.2) (46). Subsequently, we conducted a correlation analysis between drug sensitivity and the genes that control immune responses. Furthermore, we computed the dissimilarities in drug sensitivity among patients with hot and cold tumors. This investigation of the correlation between drug sensitivity and immune regulatory gene expression in various patient groups and subtypes was conducted with the aim of identifying potential treatments to provide new insights into personalized therapeutic strategies for patients with PAAD. The tumor immune dysfunction and exclusion (TIDE) algorithm was used to model the tumor immune evasion of hot/cold-immune tumors (47). The processed RNA expression levels of five patients with cancer were uploaded to the online TIDE database website (http://tide.dfci.harvard.edu/↗) to derive the TIDE score of each patient for predicting immunotherapy response.

The “ESTIMATE” R package (v 1.0.13) (48) was further applied to calculate the three immune-related scores, namely StromaScore, ImmuneScore, and ESTIMATEScore, for each patient. Additionally, we conducted an analysis of the anticancer immune response using the Tracking Tumor Immunophenotype (TIP) database (http://biocc.hrbmu.edu.cn/TIP/↗) (49).

The protein structures corresponding to several identified genes were downloaded from the Protein Data Bank database (https://www.rcsb.org/↗) (50) and pretreated with the UCSF Chimera (v 1.15). This included adding hydrogen, assigning partial charges and protonation states, and energy minimization (51). The chemical structure of active drug compounds was downloaded from the ZINC15 database (https://zinc15.docking.org/↗) (52). All compounds were subsequently docked into the binding sites of target proteins using the software DOCK (v 6.10), and visualized using the UCSF Chimera (v 1.14) and LigPlus (v 2019).

Paraffin-embedded sections of PAAD were obtained from the Affiliated Hospital of Nantong University (Nantong, China). This study was approved by the Ethical Committee of the Affiliated Hospital of Nantong University. Immunofluorescent analysis was performed to identify and assess the protein expression of pancytokeratin (PANCK; tumor epithelium marker), CD8, 5’-nucleotidase ecto (NT5E; CD73), and CD163 (M2 macrophage marker) in tumor tissues. All antibodies were purchased from AiFang Biological, Changsha City, China (product numbers AF20164, AF20211, AF301239↗, and AF20010 for PANCK, CD8, NT5E, and CD163, respectively). Paraffin-embedded sections were treated with xylene, briefly washed in a graded series of ethanol (100%, 95%, 85%, 80%, 75% ethanol), and subsequently washed with distilled water. The tissue slices were treated with antigen retrieval buffer with ethylenediaminetetraacetic acid (pH 9.0), and heated in a microwave oven. Next, the slices were immersed in a 3% hydrogen peroxide solution and incubated at room temperature for 15 min to eliminate the endogenous peroxidase. Subsequently, they were washed with phosphate-buffered saline (pH 7.4) in a decolorization shaker for 5 min, and treated with goat serum for blocking at room temperature for 30 min. Thereafter, staining was performed according to the instructions provided by the manufacturer (3-Color Multiple fluorescence Kit; AiFang Biological). Images were captured using the ECLIPSE Ci series microscope (Nikon, Tokyo, Japan), and were analyzed with the HALO image analysis platform (Indica labs, Albuquerque, NM).

R 4.0.5 software (R Project for Statistical Computing, Vienna, Austria) was used for data processing, statistical analysis, and plotting. The correlation between two continuous variables was evaluated using Pearson’s correlation coefficients, while the chi-squared test was used to compare categorical variables, and the Wilcoxon rank-sum test or t-test was used to compare continuous variables.

The CIBERSORT algorithm was used to assess the infiltration levels of 22 immune cell types and subsequently perform a consensus cluster analysis. The 33 types of cancer samples were grouped into 2–9 clusters based on the type of cancer. Cumulative distribution function curves of the consensus score matrix and proportion of ambiguous clustering statistics were used to determine the optimal number for the 33 cancer types (Figure 1A). Using the Cox model, we estimated the survival probability and determined the significance of clustering in each cancer type. Kaplan–Meier (KM) survival analysis revealed statistically significant findings for only eight of the 33 cancer types: skin cutaneous melanoma (SKCM), bladder urothelial carcinoma (BLCA), brain lower grade glioma (LGG), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), kidney renal papillary cell carcinoma (KIRP), PAAD, thymoma (THYM), and sarcoma (SARC) (Figure 1B).

We were able to differentiate the aforementioned clusters into hot and cold tumor types by considering the infiltration of CD8+ T cells, activated NK cells, and M2 macrophages, along with the scores of cytolytic activity and T cell proliferation. Increased numbers of CD8+ T and activated NK cells, higher scores of cytolytic activity and T cell proliferation, and fewer M2 macrophages were observed in hot tumors compared with cold tumors in BLCA (Figure 1C; Supplementary Figure 1A), CESC (Figure 1C; Supplementary Figure 1B), PAAD (Figure 1C, Supplementary Figure 1C), SARC (Figure 1C, Supplementary Figure 1D), and SKCM (Figure 1C, Supplementary Figure 1E). These five types of cancer all exhibited a strong consistency in the cold versus hot immune states, whereas the remaining three types had only a few distinct features; for example, hot LGG tumors exhibited lower cytolytic activity and T cell proliferation than cold LGG tumors (Figure 1C; Supplementary Figures 1F–H). In addition, there was no significant disparity in the abundance of CD8+ T cells between hot and cold THYM tumors, while several immune-related scores were inversely correlated with the presence of hot THYM tumors. KM analysis indicated that hot tumors of seven of the eight cancer types were associated with better survival than cold tumors, whereas the opposite was true for KIRP (Figure 1D). These observations led to the selection of BLCA, CESC, PAAD, SARC, and SKCM for further investigation, due to their manifestation of a comparable hot/cold immune state, with a more favorable prognosis observed for hot tumors. To verify the stability and robustness of the results obtained using the CIBERSORT algorithm, we used five additional algorithms, including TIMER, EPIC, MCP-counter, xCell, and quanTIseq, to ensure that the two consensus clusters were not biased by the analytical algorithm (Figure 1E).

To further verify the distinct immune compositions in hot and cold tumors and their implications on patient survival, we conducted an mIHC analysis using a tissue array of 71 patients with PAAD to identify the expression levels of CD8 (a marker of CD8 T cells), CD163 (a marker of M2 macrophages), and PANCK (a marker of cancer cells) (Figures 2A–D). This study revealed that CD8+ cells were sparsely distributed in the stromal regions, and were scarcely detected in some PAAD samples. Additionally, the total number of CD8+ T cells was found to be inversely correlated with overall survival (Figure 2E), indicating a compromised CTL response in PAAD. CD163+ cells were identified in both the tumor and stroma regions, as depicted in Figure 2F. Remarkably, an increased presence of CD163+ cells in tumor regions was correlated with reduced survival rates, whereas a greater abundance of stromal CD163+ cells correlated with improved patient survival, as shown in Figure 2G. These findings suggest that the CD163+ M2 macrophages within tumor regions are primarily responsible for protumoral functions. As expected, PANCK+ cells were exclusively identified in the tumoral region (Figure 2H). Interestingly, a decrease in PANCK+ cell numbers was found to be correlated with enhanced survival rates (Figure 2H), indicating that a lower tumor purity or a higher stroma component within the TME may be linked to the presence of hot tumors in PAAD. According to these findings, we suggest that the levels of M2 macrophage infiltration and tumor purity, which can be simply assessed through routine immunohistochemistry (IHC) in clinical settings, may serve as biomarkers for differentiating hot and cold tumors and predicting the prognosis of patients with PAAD.

The distribution of somatic mutations in the five aforementioned cancer types was further examined, allowing comparisons of the mutation frequencies in hot and cold tumors (Figure 3). Only the top 15 genes in each cancer type were shown due to their high mutation frequency. Genes with a mutation frequency exceeding 20% were as follows: TTN (58.3%) and RYR2 (21.7%) in BLCA, OBSCN (21.5%) in CESC, KRAS (88.4%), TP53 (74.0%), SMAD4 (27.4%) in PAAD, TP53 (53.0%) and ATRX (21.2%) in SARC, HYDIN (45.7%), MXRA5 (39.5%), ADAM18 (27.6%), TACC2 (27.0%), EPHA6 (23.4%), FREM1 (22.0%), F8 (21.4%), XDH (21.4%), and DOCK3 (21.1%) in SKCM. We identified a significant difference between hot and cold tumors in the mutation frequency of all top 15 genes in BLCA and CESC, with hot tumors containing more mutations (Figures 3A, B). In PAAD, KRAS was the only gene more frequently mutated in cold tumors than in hot tumors (Figure 3C). In SARC, hot tumors had more mutations in FCGBP, but fewer mutations in MUC16 and ADGRV1, than cold tumors (Figure 3D). In SKCM, TACC2 and HEPHL1 were more frequently mutated in hot tumors, whereas most other genes were more frequently mutated in cold tumors (Figure 3E).

We further determined the correlations between diverse immune cells in the five abovementioned cancer types, identifying a different pattern between hot and cold tumors for each cancer type (Figures 4A–E; Supplementary Figures 2A–E). For example, CD8+ T cells in hot BLCA tumors were positively correlated with activated NK cells, whereas this correlation was insignificant in cold BLCA tumors (Figure 4A). Follicular helper T cells universally exhibited a positive correlation with CD8+ T cells, particularly in cold tumors (Figures 4A–E). The TIDE score has been widely used to predict resistance to immunotherapy, with higher TIDE scores indicating a higher potential for immune escape and lower immunotherapy response rates (47). We further calculated the TIDE score to predict responsiveness to immunotherapy, finding that hot tumors were more likely to respond to immunotherapy than cold tumors in all five of the examined cancer types (Figures 4F–J). We further examined the crosstalk of certain critical stimulatory and inhibitory components of cancer immunity in hot and cold tumors. These findings suggested that CD8+ CTLs were positively correlated with most of the immune suppressive signatures, such as MDSC presence, T cell co-inhibition, antigen-presenting cell co-inhibition, and inflammation-promotion (Supplementary Figures 2F–J), suggesting that these factors may collectively protect tumor cells from attack by CTLs.

Furthermore, we evaluated the activities of multiple steps in anticancer immune response using the TIP database. These findings suggested that hot tumors generally displayed higher activity scores than cold tumors, particularly in SARC and SKCM (Figures 4K–O). Further, the priming stages of anticancer immune response were significantly inhibited in cold tumors. For example, hot tumors from multiple cancer types exhibited higher activity scores than cold tumors for step 1 (release of cancer cell antigens), step 2 (cancer antigen presentation), and step 3 (priming and activation). Further, in step 4 (trafficking of immune cells to tumors), we consistently found that hot tumors exhibited higher scores for various immune cell types, including CD8+ T cells and NK cells. In step 5 (infiltration of immune cells into tumors), only cold SKCM tumors exhibited lower activity scores than hot tumors, whereas in the other four types of cancers, cold and hot tumors displayed similar scores, indicating that most cold tumors were not inherently immune-excluded relative to hot tumors. The scores for step 7 (activity of killing cancer cells) were higher in PAAD and SKCM hot tumors, and were comparable between BLCA, CESC, and SARC hot and cold tumors. Interestingly, the scores for step 6 (recognition of cancer cells by T cells) were consistently lower in hot tumors than cold tumors, suggesting that insufficient neoantigen recognition by T cell receptors may be the critical immunosuppression factor in hot tumors despite the increased T cell recruitment (Figures 4K–O).

To better describe the immune landscape of each tumor type, we developed an immunogram by incorporating seven antitumoral immune parameters (CD8+ T cells, activated NK cells, follicular helper T cells, T cells proliferation, cytolytic activity, immunogenic cell death, MHC class I) and five protumoral immune parameters (M2 macrophages, MDSCs, Tregs, T cell co-inhibition, inflammation-promotion). Consequently, a 12-axis radar plot was generated to visualize the immune state of hot and cold tumors (Figures 5A–E). We further found that only the values of axis 7 (M2 macrophages) were higher in all cold tumors (Figures 5A–E), while those of axis 8 (MDSCs) were higher in cold tumors of the BLCA and SARC groups (Figures 5A, D). Surprisingly, the values of axis 11 (MHC class I) were higher in cold tumors of PAAD (Figure 5C), although the difference was not statistically significant (Figures 5F–J). Concerning other parameters, higher values were consistently observed in hot tumors, regardless of cancer type. Collectively, these results indicate a similar pattern in all five cancer types, i.e., greater numbers of CD8+ T cells and follicular helper T cells and fewer M2 macrophages in hot tumors compared with cold tumors (p < 0.0001) (Figures 5F–J). Of note, the pattern was not shared by the other three types of cancer, namely KIRP, LGG, and THYM (Supplementary Figures 3A–C).

We selected 150 immune regulator genes and 60 immune checkpoint genes based on the results of prior studies conducted by Shen et al. (44) and Thorsson et al. (45), and compared their expression between our hot and cold tumor datasets (Figure 6A). The findings revealed that 22 genes, including key immune checkpoint molecules such as PDCD1, TIGIT, and LAG-3, were significantly upregulated in hot tumors across all five cancer types (Figure 6B). Additionally, the expression of CD276 was found to be elevated in cold tumors across all five cancer types, while NT5E exhibited higher expression levels in cold tumors of the CESC, PAAD, SARC, and SKCM groups (Figure 6C). Other DEGs between hot and cold tumors are presented in Figure 6C.

Next, GSVA was employed to investigate the differential enrichment of hallmark pathways in hot and cold tumors (Supplementary Figures 4A–E). Our analysis revealed that the apoptosis and bile acid metabolism pathways were enriched in cold tumors of BLCA, CESC, SARC, and SKCM (Figure 7A); while the angiogenesis, apical junction, and apical surface pathways were enriched in cold tumors of BLCA, CESC, PAAD, and SARC (Figure 7A). Correlation analysis was further performed using the 23 shared DEGs (22 upregulated plus one downregulated gene in hot tumors compared with cold tumors from all five cancer types), in addition to the GSVA scores of the five enriched pathways. Interestingly, almost all DEGs were positively correlated with the apical surface and apoptosis pathways in all cancer types, indicating that these genes were actively involved in the regulation of these two pathways. In BLCA, numerous DEGs, except for CD160 and CD96, were positively correlated with the angiogenesis and apical junction pathways; however, this correlation was not observed in the other four types of cancer (Figure 7B). We further performed GSVA on the KEGG pathways in hot and cold tumors separately, the results of which are shown in Supplementary Figures 5A–H.

Correlations between immune cell infiltration and the levels of previously indicated immune regulatory genes (44, 45) were further examined (Supplementary Figure 6). The genes with strong correlations to the indicated immune cell types are illustrated in Figure 8A. Most selected genes were positively correlated with CD8+ T cells in all cancer types, except for CD276, which was negatively correlated with CD8+ T cells. CD276 was also negatively correlated with activated NK cells, but positively correlated with M2 macrophages in most cancer types. Positive correlations between most genes and activated NK cells were observed in BLCA, CESC, SARC, and SKCM. Nevertheless, an opposite trend was observed in PAAD, indicating an impaired recruitment of activated NK cells under these conditions. M2 macrophages were also found to be negatively associated with most genes in SKCM, PAAD, SARC, and CSEC (Figure 8A). The Venn plots in Figure 8B demonstrate the distribution of DEGs positively or negatively associated immune cells, specifically CD8+ T cells, M2 macrophages, activated NK cells, and follicular helper T cells. In PAAD, we further found that NT5E and VEGFA may negatively regulate CD8+ T cells, while CCL13, CCL18, NT5E, and TNFSF4 might positively regulate M2 macrophages (R > 0.2). Other genes that may regulate other immune cells are listed in Supplementary Table 1.

Overall, in the aforementioned experiments, we demonstrated that cold tumors in PAAD express higher levels of NT5E and CD276. Interestingly, a robust linear relationship between the expression of NT5E and CD276 was noted in PAAD (R = 0.54, p = 2.6e-14) (Figure 9A). In addition, levels of NT5E and CD276 were positively associated with tumor grades, but not other clinicopathological characteristics; notably, tumors from elderly patients expressed higher levels of CD276 (Figure 9B). We further examined the links between tumor heterogeneity and the expression of NT5E and CD276, finding positive correlations between NT5E or CD276 expression and the degree of homologous recombination deficiency and loss of heterozygosity (Supplementary Figures 7A–C).

KM survival analysis of the TCGA dataset and the other two independent PAAD cohorts revealed that NT5E and CD276 negatively influenced patient survival (Figures 9C–H). Subsequently, patients with PAAD were further classified into four groups based on their immune subtypes and NT5E expression: NT5Ehigh-hot, NT5Ehigh-cold, NT5Elow-hot, and NT5Elow-cold. Subsequent KM analysis revealed that the NT5Elow-hot group had the best survival rate among these groups (Figures 10A, B). Similarly, the CD276low-hot group was associated with a better survival rate compared with that of the CD276high-hot, CD276high-cold, and CD276low-cold groups (Figures 10C, D). GSEA between the NT5Ehigh-cold and NT5Elow-hot groups (Figure 10E), and between the CD276high-cold and CD276low-hot groups (Figure 10F) revealed that the hypoxia pathway was enriched in the NT5Ehigh-cold and CD276high-cold groups (Figures 10E, F). Positive correlations between the hypoxia pathway and NT5E or CD276 were prominent, as shown in Figures 10G, H. As expected, KM survival analysis showed that hypoxia negatively influenced the prognosis of PAAD (Figure 10I). Moreover, we identified significant differences in overall survival between the NT5Ehigh-Hypoxiahigh and NT5Elow-Hypoxialow groups (Supplementary Figure 8A), as well as between the CD276high-Hypoxiahigh and CD276low-Hypoxialow groups (Supplementary Figure 8B).

Subsequently, multiple single-cell RNA sequencing datasets of PAAD were applied to examine the expression patterns of NT5E and CD276 across different cell populations within the TME. The findings revealed predominant expression of CD276 and NT5E in monocytes/macrophages, fibroblasts, and malignant cells, as illustrated in Figures 11A, B. Additionally, protein expression and subcellular localization of CD276 and NT5E were examined using the results of immunohistochemistry analysis of PAAD tissue sections extracted from the HPA database (Figure 11C). Next, an mIHC assay was conducted to examine NT5E protein expression and its correlation with hot and cold TME in patients with PAAD, as well as its impact on patient survival. This assay detected NT5E, CD163, and PANCK in the PAAD tissue array, revealing the presence of NT5E+ cells in both tumoral and stromal regions (Figures 2A–E). Moreover, a positive correlation was observed between the numbers of NT5E+ cells and CD163+ cells in tumor regions, indicating a potential association between increased NT5E expression and the cold TME in PAAD (Supplementary Figure 9). Ultimately, we demonstrated that a greater abundance of NT5E+ cells, irrespective of their spatial distribution, were correlated with diminished patient survival in PAAD, thereby confirming the findings of transcriptomic analysis (Figure 2H). The overview image of mIHC staining of the PAAD tissue array is shown in Supplementary Figure 10.

Herein, we sought to investigate the potential of the DEGs between hot and cold tumors as candidate drug targets in PAAD. Therefore, we employed the oncoPredict to forecast gene correlations with drug responses. Our findings revealed that two drugs (dasatinib and tozasertib) exerted a positive influence on the majority of genes found to be upregulated in hot tumors. Importantly, these agents also negatively influenced most genes that were upregulated in cold tumors, which is a desirable effect (Figure 12A). The regulatory mechanisms underlying the effects of these two drugs may differ between hot and cold tumors, as specific genes, such as endothelin receptor type B (EDNRB), CCL14, and C-X-C motif chemokine ligand 12 (CXCL12), were regulated more strongly by both drugs in hot tumors (Figure 12B). In addition, both drugs negatively affected the expression of NT5E in cold tumors, but minimally affected its levels in hot tumors (Figure 12B). The effects of these two drugs were more marked in cold tumors versus hot tumors (Figure 12C), as indicated by the lower IC₅₀ values. This evidence suggested that these agents could exert more beneficial effects on cold tumors. The beneficial effects of these two drugs were further indicated by the positive correlations between the differentially expressed immune genes and the ImmuneScore, StromaScore, and ESTIMATEScore (three scores used to estimate the extent of tumor immune infiltration) (Figure 12D). The effects on various immune cells were further evaluated, and the results indicated that both drugs positively regulated T cell subtypes, particularly CD4+ memory-activated cells in hot tumors (Figures 12E, F).

We further downloaded the chemical structures of dasatinib and tozasertib from the ZINC15 database. The complete protein structures of the products of the 44 DEGs between hot and cold PAAD were successfully retrieved and preprocessed (Supplementary Table 2). Next, we explored the binding potential between proteins and drugs. According to the docking scores (Supplementary Table 2), the top three molecules that exhibited the most potent binding affinity with dasatinib were PRF1, CXCR6, and ADORA2A. Further, 2D and 3D molecular visualization of drug-protein interactions demonstrated that dasatinib displayed strong interactions with the PRF1 amino acid residue Ala121 (Figure 13A) and the CXCR6 residue Tyr278 (Figure 13B). Furthermore, pockets were identified on the surface of the PRF1, CXCR6, and ADORA2A which would allow the formation of a stable complex with dasatinib (Figures 13A–C). The top three molecules exhibiting the most robust binding affinity with tozasertib were ADORA2A, PRF1, and CCR6 (Figures 13D–F). Tozasertib exhibited robust interactions with the PRF1 residues Ser234, Arg232, and Asp120 amino acids, as well as the Thr220 residue of CCR6 through hydrogen bonds (Figures 13E, F). In addition, the pockets on the surface of ADORA2A and CCR6 interacted with tozasertib to form a complex (Figures 13D, F).