3.1 The differential mRNA expression and prognostic value of chemokine genes
A total of 41 chemokine genes were identified in this study, the differential expression of which we initially examined in 33 tumor types based on information obtained from the TCGA and GTEx databases. As shown in Fig. 1A, with the exception of pheochromocytoma (PCPG) and sarcoma (SARC), chemokine genes were commonly differentially expressed in most tumor types. CCL18 was highly expressed in 28 of the 33 tumors, whereas CCL16 was expressed at low levels in 23 of the 33 tumors. Using TCGA pan-cancer data, we also performed mRNA correlation analysis (Fig. 1B), and also conducted pan-cancer univariate analysis to determine the prognostic value of each gene (Fig. 2A). The results thus obtained indicated that CCL20 was a risk factor in nine of the 33 tumors, whereas CCL17 was a protective factor in seven
(Fig. 2B).
3.2 The alteration of chemokine genes
For each of the assessed chemokine genes, we obtained the following information: variant classification, variant type, SNV class, variants per sample, variant classification summary, and the top 10 mutated genes (Fig. 3A). With regards to the SNV percentage, the total deleterious mutation percentage (the number of samples with at least one deleterious mutation site/the number of samples with SNV mutation data) revealed the highest deleterious mutation frequency in CXCL1 and CXCL12 in uterine corpus endometrial carcinoma (UCEC) (Fig. 3B). Moreover, pan-cancer analysis revealed that CX3CL1 has the highest mutation frequency (10%) among the assessed chemokine genes (Fig. 3C).
We further analyzed the copy number variants (CNVs) of chemokine genes, the proportions of different types of which, including heterozygous amplification, heterozygous deletion, homozygous amplification, and homozygous deletion, in pan-cancer are shown in Fig. 4A. We accordingly found that the CNV of CXCL16 was positively correlated with the mRNA expression of this gene in 22 of the 33 tumor types (Fig. 4B). Furthermore, the CNV of homozygous or heterozygous amplification was established to be positively correlated with mRNA expression, whereas the CNV of homozygous or heterozygous deletions was negatively correlated with mRNA expression (Fig. 4C and D).
We also determined the methylation status of chemokine genes in the 33 tumor types, with the results indicating that the methylation levels were negatively correlated with their mRNA expression (Fig. 5A). The differential methylation levels of genes between the tumor and normal sample groups are shown in Fig. 5B.
3.3 The calculation, differential distribution, and survival analysis of chemokine scores
We performed ssGSEA to calculate chemokine scores for the 33 tumor types identified in the TCGA cohort. These chemokine scores were found to be highest and lowest for lung adenocarcinoma (LUAD) and acute myeloid leukemia (LAML), respectively (Fig. 6A). Moreover, we found that chemokine scores were lower for tumor tissues than in the adjacent normal tissues with respect to bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), kidney chromophobe (KICH), liver hepatocellular carcinoma (LIHC), lung squamous cell carcinoma (LUSC), and UCEC (Fig. 6B-G), although were higher in esophageal carcinoma (ESCA), kidney renal clear cell carcinoma (KIRC), and stomach adenocarcinoma (STAD) (Fig. 6H-J).
The results of uniCox analysis revealed the following: (1) for OS, chemokine score was a risk factor in KIRC, UVM, THYM, GBM, LAML, and PAAD, and a protective factor in SKCM, SARC, BRCA, and OV (Fig. 7A); (2) for DSS, the chemokine score was a risk factor for KIRC, UVM, THYM, and GBM, and a protective factor in SKCM, SARC, and OV (Fig. 7B); (3) for DFI, the chemokine score was a protective factor in LIHC, COAD, and BLCA (Fig. 7C); and (4) for the PFI, the chemokine score was a risk factor for KIRC, GBM, and THYM, and a protective factor in LIHC, SKCM, ACC, BRCA, and HNSC (Fig. 7D).
3.4 Gene set variation analysis of chemokine score
To analyze the pathways potentially associated with chemokines, we performed GSVA based on 50 HALLMARK pathways. The correlations between chemokine expression and GSVA scores in pan-cancer are shown in Fig. 8. We observed that chemokine scores were positively associated with numerous immune-related pathways in pan-cancer, including IL6-JAK-STAT3 signaling, inflammatory response, IL2 STAT5 signaling, and interferon gamma response.
3.5 Relationships between chemokines and the tumor microenvironment
We also examined the association between chemokine scores and the stromal and immune scores in pan-cancer (Fig. 9A), the results of which revealed that the chemokine score was positively associated with the immune, stromal, and ESTIMATE scores in pan-cancer. We further obtained and analyzed TME-related pathways based on reference to a previously published paper (13), including immune-related, stromal-related, and DNA repair-related pathways (Fig. 9B). For example, in BLCA, immune-related signaling (immune checkpoint, CD8 T effector, and antigen processing machinery) was higher in the high-chemokine score group (Fig. 9C).
These findings indicate that chemokines are highly correlated with immune-related pathways, and we accordingly further examined the correlations between chemokine scores and immune cells in the tumor immune microenvironment. On the basis of the immune cell infiltration results obtained from TIMER2, we found chemokine scores to be positively correlated with the infiltration level of most immune cells (Fig. 10A), and reference to the ImmuCellAI database indicated that, with the exception of CD8 naive cells and neutrophils, chemokine scores were highly correlated with most immune cells at the pan-cancer level (Fig. 10B). We also analyzed the correlations between chemokine scores and immune regulation-related genes and observed a positive correlation in these gene sets, including immunosuppressive (Fig. 11A), immune-activating (Fig. 11B), and MHC (Fig. 11C) genes. These results indicated that tumor samples with high chemokine scores were rich in immune cells.
3.6. The association between chemokines and immunotherapeutic response
On the basis of our observations indicating that chemokine scores were positively correlated with most immune cell types, we hypothesized that chemokines play a significant role in determining the efficacy of immunotherapy and that patients with a high chemokine score may be sensitive to ICI treatment. To verify these suppositions, we obtained immunotherapy data and calculated the corresponding chemokine scores, which revealed that GSE176307 cohort patients with high chemokine scores had better survival status (Fig. 12A). Furthermore, Kaplan–Meier analysis of the different groups revealed that the proportion of SD/PD patients in the high-chemokine group (71%) was lower than that in the low-chemokine group (85%) (Fig. 12B), and similar results were obtained for the ICB.Riaz2017 (Fig. 12C and D) and IMvigor210CoreBiologies (Fig. 12E and F) cohorts.