Immune Cluster and PPI Network Analyses Identified CXCR3 as a Key Node of Immunoregulation in Head and Neck Cancer

Abstract

The tumor microenvironment (TME) is significantly associated with clinical outcomes and therapeutic efficacy. However, the landscape of the head and neck cancer (HNC) microenvironment is not fully understood. Therefore, we divided HNCs into three classes according to differences in the TME to determine effective personalized treatments. We explored the immune landscape of head and neck cancer by analysing the gene expression profile of 501 cases from the Cancer Genome Atlas (TCGA) data portal and validated our findings in 270 cases from the Gene Expression Omnibus (GEO) database. The levels of immune components in the tumor microenvironment were evaluated via single-sample gene set enrichment (ssGSEA) analysis. The HNCs were clustered into an Immunity-H group, Immunity-M group and Immunity-L group according to 40 immune components in the tumor microenvironment. DNA damage and HLA genes play an important role in immune regulation. The patients in the Immunity-H group had a favourable survival compared with patients in the Immunity-M group and the Immunity-L group. The patients in the Immunity-H group and Immunity-M group could benefit from radiotherapy. In addition, the Immunity-L group showed the lowest immunophenoscore and had poor response to anti-PD-1 treatment. CXCR3 was demonstrated to be downregulated in the Immunity-L group, which was related to shorter OS in the TCGA and GEO databases, suggesting CXCR3 as a potential therapeutic target. Taken together, our findings proposed three new microenvironment-related phenotypes of HNCs and suggested that CXCR3 played a major role in immune regulation and could be a novel therapeutic target, providing a reference for clinical decisions and research directions in the future.

Introduction

Head and neck squamous cancer (HNSCC) ranks seventh in global cancer incidence, and 890,000 people were diagnosed with HNSCC in 2018 (1). On the basis of traditional treatments, including surgical resection and chemoradiation, the 5-year survival rate of HNC remains 40%–50% (2). Other than cetuximab, a monoclonal antibody (mAb) targeting epidermal growth factor receptor (EGFR), no other new targeted therapies have been approved for HNSCC for decades. Cetuximab monotherapy efficiency is only 10%~15%, and there are no known biomarkers for predicting response (3, 4). At present, immunotherapy has been approved by the FDA and EMEA to treat recurrent and metastatic patients with HNSCC (5). However, regardless of the type of treatment, some people do not benefit or experience associated side effects. Therefore, it is crucial to explore new therapeutic targets and identify therapeutic methods that are suitable for specific groups of people.

The tumor microenvironment (TME) plays an important role in cancer growth, metastasis and response to therapy. The tumor microenvironment is composed of cancer cells, immune cells, stromal tissues, the extracellular matrix and other components (6). Chemokine and cytokine signalling in the TME regulates tumor behaviour and response to therapy, affecting the interactions among immune cells (7). Immune cells can change their status according to different tumor microenvironments. Sometimes, immune cells fail to clear tumor cells due to the immunosuppressive status of the tumor microenvironment (8).

In our study, we aimed to apply the ssGSEA method to evaluate the levels of various immune components, including immune cells, factors and pathways, in the tumor microenvironment and identify people with different immune statuses through k-means clustering to explore the optimal treatment plans in HNSCC.

Materials and Methods

Data Collection and Collation

We collected the gene expression data and clinical data of the TCGA-HNSC cohort from https://genome-cancer.ucsc.edu/. The CNV data and mutation data of these patients with head and neck squamous cancer were downloaded from the publicly available TCGA database via the GDC Data Portal. We ultimately analyzed 501 TCGA-HNSC patients after excluding patients who did not have survival data.

Single-Sample Gene Set Enrichment Analysis (ssGSEA) and Clustering

We obtained marker gene sets for 40 immune components, including immune cells, immune-related factors and immune pathways, in the tumor environment from a previously published study (9) and the ImmPort database (10). We calculated the enrichment scores of 40 immune-related signatures in 501 patients by using the R package “gsva” (11). According to different immune scores, we used k-means clustering to divide these patients into three groups: the Immunity-high group, the Immunity-moderate group and the Immunity-low group. The immune scores, stromal scores and tumor purity were calculated based on the ESTIMATE algorithm (12). The TIMER database (https://cistrome.shinyapps.io/timer/) was used to calculate the infiltration levels of six immune cells, including CD4⁺ T cells, CD8⁺ T cells, B cells, macrophages, neutrophils, and dendritic cells(DCs) according to different gene sets (13, 14).

Calculation of DNA Damage and Immunogenomic Indicators

Aneuploidy, HRD (homologous recombination deficiency), CNA burden (copy number variation burden), ITH (intratumor heterogeneity), and SNV-related neoantigen data were obtained from a published study (15). The tumor mutational burden (TMB) was calculated as the number of nonsynonymous protein-coding variants divided by the total sequenced genome length according to the mutation data (16).

Response to Therapy

The patients were divided into two groups via k-means clustering based on the expression level of 31 genes associated with radiosensitivity (17). A radioresistant group (RR group) and a radiosensitive group (RS group) were identified based on prognosis under radiotherapy. The immunophenoscore (IPS) (18) and chemosensitivity (19) were determined from previously published studies.

DEG Screening, PPI Network Construction, and Hub Gene Identification

DEGs between the Immunity-H group and Immunity-L group were identified with the “edgeR” package (20) and the criterion |logFC| > 1 and P < 0.05. A web tool, the STRING database (https://string-db.org/), was used to build a protein-protein interaction network for the DEGs (21), and the minimum required interaction score was set as 0.400. The top 10 hub genes were identified by the “cytohubba” tool in Cytoscape software 3.6.1.5 (22, 23).

Gene Ontology and KEGG Pathway Enrichment Analyses

The associated molecular functions (MFs), biological processes (BPs), cellular components (CCs) and pathways for the top 10 hub genes were annotated by the R packages “clusterProfiler” and “org.Hs.eg.db” (24). A false discovery rate (FDR) value <0.05 was set as the enrichment cut-offs to screen for meaningful enrichment results. The enrichment results were visualized via the R package “ggplot2R”. Survival analysis of the top 10 hub genes was performed using the R software package “survival”.

Immunohistochemistry

A total of 53 biopsy samples from patients with primary head and neck cancer were collected at the First Hospital of China Medical University between 2005 and 2011. Using 4 μm-thick sections, immunohistochemistry was performed. The antibody used was rabbit anti-human CXCR3 polyclonal antibody (BA0759, 1:300 dilution; WB-BIO, Wuhan, China). The staining intensity and percentage of positive cells was estimated. The proportion of positive cells greater than 5% is considered positive based on previous studies (25).

Correlation Analysis

We evaluated the correlation between the expression of CXCR3 and other hub genes in HNSCC through GEPIA (http://gepia.cancer-pku.cn/index.html) based on the TCGA database (26). Pearson’s correlation coefficient analysis was used to define correlations.

CNV Analysis

The copy number variation score processed by GISTIC2.0 was downloaded from the TCGA database via the GDC Data Portal. Segment values larger than 0.3 were defined as “1”, and those less than -0.3 were defined as “-1”; segment values between -0.3 and 0.3 were defined as “0”. The relationship between gene copy number variation and immune infiltration was analyzed with TIMER (27).

GEO Validation Dataset

Validation data (GSE65858 and GSE41613) were downloaded from the Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo). The GSE65858 data were analyzed to validate the immune classes. The GSE41613 data were analyzed to validate the prognostic value of CXCR3.

Statistical Analysis

The chi-square test was adopted to identify the significantly different copy number variants (CNVs) between the Immunity-H group and the Immunity-L group. An independent samples t-test was employed to compare two groups, while one-way ANOVA test was employed to compare multiple groups. All analyses were performed with R software (Version 3.6.1) and IBM SPSS Statistics software (Version 26). P-values <0.05 were considered statistically significant.

Results

Identification and Validation of the Immune Class Clustering by 40 Immune Components in the HNSCC Environment

Immune cells, immune factors and immune pathways participate in the immune response together to prevent the occurrence and progression of tumors in the immune microenvironment of HNSCC. We evaluated the levels of 40 immune components, including immune cells, immune factors, and immune pathways, by the ssGSEA method according to the transcriptomes of 501 TCGA-HNSC patients. These patients were clustered into three clusters (Immunity-H: 264; Immunity-M: 185; Immunity-L: 52) based on the 40 immune components by k-means clustering (Figure 1A). To validate the reliability of the immune classes, we first compared the immune scores, stromal scores and tumor purity among the three immunity classes. The Immunity-H group showed the highest immune scores and stromal scores and the lowest tumor purity; in addition, the Immunity-L group showed the lowest immune scores and stromal scores and the highest tumor purity (Figure 1B). The distributions of six immune cells calculated in the TIMER database were in agreement with our immune classes (Figure 1C). We further validated our expression profile-based clustering by the same approach in 270 samples from GEO datasets (Figure 1D). The immune scores, stromal scores, and tumor purity and the distributions of six immune cells were consistent with the immune clustering (Figures 1E, F). Overall, we demonstrated that the HNSCC samples had three different immune statuses when assessed by 40 immune components.

Figure 1

HLA Expression and DNA Damage Were Associated With Immune Class

Considering previous publications (15, 28), there are at least three aspects associated with immune status: HLA gene expression, DNA damage and immunogenicity. HLA (human leukocyte antigen) on the surface of many immune cells plays an important role in activating cellular and humoral immunity (29). We compared HLA gene expression among the three classes and found that all HLA genes were expressed at a high level in the Immunity-H group and at a low level in the Immunity-L group (P<0.001) (Figure 2A), validating that the absence of the HLA gene contributed to the low immune status. In addition to HLA gene expression, DNA damage and immunogenicity are also associated with the immune response (15). To estimate whether DNA damage and immunogenicity affected the immune status of the HNSCC samples, we investigated the correlation between the immune class and these factors. DNA damage can be considered in 4 major categories: aneuploidy, homologous recombination deficiency, copy number variation burden and intratumor heterogeneity. Low immune status was significantly associated with high aneuploidy score, high HRD and high CNA burden but had no relevance to ITH (Figures 2B–E). For immunogenicity, the differences in the tumor mutational burden and SNV-related neoantigen burden were not as large as those in homologous recombination deficiency and copy number variation burden (Figures 2F, G). In general, HLA expression and DNA damage were significantly associated with the immune class.

Figure 2

Prognostic Value of 40 Immune Components in HNSCC

We performed survival analysis among the three groups by different methods, and the results demonstrated that the immune classes were associated with prognosis. Kaplan–Meier analysis showed that the Immunity-H group had better survival than the Immunity-M group and Immunity-L group (P<0.01) (Figure 3A). The multivariate Cox proportional hazards model also revealed that the Immunity-M group and Immunity-L group independently predicted worse overall survival in HNSCC than the Immunity-H group (Table 1). In addition, we investigated the association between clinical factors and immune class. There were similar clinicopathologic characteristics among the three immune classes (Figure 3B). Considering that the immune microenvironment is dynamic, we compared the prognostic value of 40 immune components in the three different immune statuses (Figure 3C). In the whole HNSCC cohort, most immune components were related with survival. However, 40 immune components showed different prognostic effects with different immune statuses, for example, the Th1 cells, Th2 cells, and macrophages were associated with prognosis only in the Immunity-H group. Moreover, neutrophils predicted a worse prognosis in the Immunity-H group and a better prognosis in the Immunity-L group and are thus worth studying further.

Figure 3

Therapeutic Strategies According to Immune Class

We further explored the therapeutic effects of common treatments in HNSCC. In the Immunity-H group and the Immunity-M group, patients treated with radiotherapy exhibited better overall survival rates than patients without radiotherapy (p<0.05) (Figures 4A, B), whereas there were no significant differences in the Immunity-L group (p=0.2613) (Figure 4C). In addition, the proportion of radiosensitive patients in the Immunity-H group and the Immunity-M group was obviously greater than that in the Immunity-L group (Figure 4D). Charoentong et al (18). proposed an IPS for defining patients likely to respond to anti-PD-1 therapy. As illustrated in Figure 4E, the IPS score decreased as the immune status decreased (p<0.001); in addition, the immunotherapy markers PD-L1, PD1, and CTLA-4 also decreased as the immune status decreased (p<0.0001) (Figures 4F–H). We next estimated the chemosensitivity in subgroups, and the log-transformed IC50 values are shown in Figure 4I. Most drugs were most effective in the Immunity-H group, including 5-fluorouracil, belinostat, bexarotene, bicalutamide, idelalisib, lenalidomide, nilotinib, and ruxolitinib. Only cetuximab, bosutinib, and vinorelbine were effective in the Immunity-L group. As a result, the low immune status predicted worse efficacy of administered treatments than the moderate and high immune statuses.

Figure 4

PPI Network Construction Revealed That CXCR3, the Top Hub Gene in Immunoregulation, Was Associated With Prognosis

To explore potential therapeutic targets to reverse the immune status in the Immunity-L group, we identified 677 downregulated and 21 upregulated genes in the Immunity-L group compared with the Immunity-H group. As shown in Figure 5A, we constructed a PPI network composed of 177 nodes and 437 edges (Figure 5A). In addition, we identified the top 10 hub genes, including CXCR3, CXCR5, CCR2, CCR8, CCL1, CCL25, P2RY12, CNR2, PNOC, and GPR31 (Figure 5B). In terms of biological processes, these hub genes were significantly enriched in leukocyte migration and cell chemotaxis. In terms of molecular functions, these hub genes were enriched in G protein-coupled chemoattractant receptor activity, chemokine receptor activity, cytokine receptor activity, etc. (Figure 5C). KEGG pathway enrichment analysis showed that the hub genes were associated with viral protein interactions with cytokine and cytokine receptors, chemokine signalling pathways and cytokine-cytokine receptor interactions (Figure 5D).

Figure 5

We further analyzed the prognostic value of the 10 hub genes. Low expression of CCR8, CCL1, CCL25, CCR2, CXCR3, CNR2, CXCR5 and PNOC was associated with worse survival (p<0.05) (Figures 5E–L), whereas P2RY12 and GPR31 expression was not associated with survival. Moreover, CXCR3, the top hub gene, was downregulated in the Immunity-L group compared with the Immunity-H and Immunity-M groups (Figure 5M), and low CXCR3 expression predicted worse prognosis than high CXCR3 expression (P<0.05) (Figure 5N), which were validated in the GEO datasets. The correlation analysis implied that CXCR3 expression was positively correlated with the other 9 hub genes, including CXCR5 (R=0.3, P<0.05), CCR2 (R=0.76, P<0.05), CCR8 (R=0.58, P<0.05), CCL1 (R=0.15, P<0.05), CCL25 (R=0.12, P<0.05), CNR2 (R=0.52, P<0.05), PNOC (R=0.52, P<0.05), P2RY12 (R=0.65, P<0.05), and GPR31 (R=0.69, P<0.05), which further revealed that CXCR3 plays a core role in immune regulation. For a more comprehensive study, we accessed the CXCR3 expression using immunohistochemistry and CXCR3 protein was mainly expressed in the cytoplasm of cell (Figures 5O, P). CXCR3 positive expression was identified in 23 cases (43.39%) and positive CXCR3 status was associated with better prognosis (Figure 5Q).

LASSO regression analysis of the top 10 hub genes showed that the expression levels of 4 genes, P2RY12, CXCR3, PNOC, and CCR8, could be used to predict low immune status. In addition, there were significant associations between CXCR3 (p<0.001), PNOC (p=0.026), and CXCR8 expression and low immune status (p=0.005) in the multivariate logistic regression model (Table 2). When using the three genes to predict the risk of patients with low immune status, the accuracy was 94.2%.

CXCR3, CXCR5, and CCL1 Copy Number Variants Correlated With Low Immune Status

In addition to transcriptome analysis, we also investigated the CNV differences of the 10 hub genes between the Immunity-L group and Immunity-H group. The CNV events related to CXCR3 on chromosome X, CCL1 on chromosome 17 and CXCR5 on chromosome in the Immunity-L group were significantly more frequent than those in the Immunity-H group (p<0.05) (Figure 6A). However, CNV events related to the other seven hub genes were not associated with immune class. In addition, we investigated the association between CNV events related to the three genes and immune infiltration in the TIMMER database. CXCR3 copy number variation was significantly associated with decreased levels of six immune cells (B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells) (Figure 6B), and CCL1 copy number variation was significantly associated with decreased levels of 4 immune cells (B cells, CD8 + T cells, neutrophils, and dendritic cells) (Figure 6C), but CXCR5 copy number variation was only associated with CD8 cell level (Figure 6D). These results showed that CXCR3 and CCL1 copy number variations significantly regulate immunity, whereas CXCR5 copy number variations may have little relation with immunity.

Figure 6

Discussion

In recent years, based on emerging evidence, solid tumors are no longer thought to exist in isolation but rather in a complicated environment called the tumor microenvironment composed of multiple cell types, including neutrophils, macrophages, regulatory T cells, and more. These cells interact with each other via various cytokines, chemokines, and growth factors (30). The purpose of our study was to determine the heterogeneity of HNSCC samples within different immune statuses through analysis of immune cells, factors, and pathways in the TME.

Using the scores for 40 immune components calculated by ssGSEA, we divided HNSCC samples into three immune classes, the Immunity-H group; the Immunity-M group; and the Immunity-L group, and confirmed the utility of the classifications via different methods and different cohorts. In previously published literature (15, 28), we found the following characteristics of the classes that might contribute to low immune status: 1, defective HLA gene expression leading to defective presentation of antigens and immune system activation; 2, DNA damage driving endogenous immune deficiency; and 3, aberrant immunogenicity. The first two reasons were confirmed in this study; nevertheless, the third one could not be confirmed in this study.

We further analyzed the prognosis of the three immune classes. The patients in the Immunity-H group had the best survival compared with patients in the Immunity-M group and the Immunity-L group. Malignant tumors promote tumor progression by suppressing effective antitumor immunity, and immune cells, factors and pathways are indispensable parts of the immune response (31). These findings could explain the different prognoses among the three immune classes. Interestingly, the 40 immune components showed prognostic effects on prognosis in different immune statuses; for example, neutrophils predicted opposite prognoses in the Immunity-H group and the Immunity-L group. Regarding this interesting phenomenon, we hypothesized that neutrophils could be differentially polarized and driven by the tumor microenvironment to transition into different subtypes with different functions in the Immunity-H group and the Immunity-L group (32). In 2019, Fridlender and colleagues (33) first suggested the description of antitumorigenic and protumorigenic neutrophils, but no definitive surface marker were identified to distinguish them. Our results further illustrated that immune cells can vary their functions to adapt to different tumor environments.

Our research also provided a guide for tailoring therapeutic strategies. First, we found that radiotherapy improved the overall survival of the patients in the Immunity-H group and the Immunity-M group, and the percentage of radiosensitive patients increased with the improvement of immune status. Similarly, most of the chemotherapy drugs were also most effective in the Immunity-H group. Recently, an increasing number of studies have indicated that immunotherapy greatly improves the prognosis of patients with advanced disease. In our research, we predicted the immunotherapy response by the anti-PD-1 score proposed by Charoentong et al. and some related markers in previous studies. Anti-PD1 treatment could revive early-stage dysfunctional T cells but not late-stage dysfunctional T cells (34), which supported our results. We found that the IPSs and immune marker levels were the lowest in the Immunity-L group compared with the Immunity-M group and the Immunity-H group, which showed that the patients in the Immunity-L group suffered from late-stage immune dysregulation and presented a low immune status because of long-term immune disorders. In general, monotherapy of radiotherapy, chemotherapy, or immunotherapy was efficacious in the Immunity-H group, and combined radiotherapy and chemotherapy or immunotherapy was suitable in the Immunity-M group; however, there were no satisfactory effects from existing treatments methods for the patients in the Immunity-L group.

To explore some potential targets to improve the prognosis of patients in the Immunity-L group, we obtained 697 DEGs via a comparison between the Immunity-H group and the Immunity-L group and 10 hub genes by PPI network analysis of these DEGs. The enrichment analysis of the 10 hub genes showed that the genes mainly participated in functions and pathways related to chemokines and cytokines. The 10 genes, CXCR3, CXCR5, CCR2, CCR8, CCL1, CCL25, CNR2, and PNOC, for which low expression was associated with worse survival were potential targets, and CXCR3 and CCL1 copy number variations were significantly associated with low immunity. Through validation analysis in the GEO database, we found that the chemokine receptor CXCR3 was the most valuable target. HNSCC samples with lower expression of CXCR3 had significantly shorter OS than HNSCC samples with higher expression of CXCR3, and this relation was validated in the GEO datasets. Similarly, Fangfang Chen et al. indicated that the overexpression of CXCR3 was associated with increased dendritic cell and tumor-infiltrating lymphocyte infiltration and improved OS in gastric cancer (35). However, some researchers found that higher CXCR3 expression was associated with more distant metastasis and shorter OS than lower CXCR3 expression (36, 37). Although the prognostic effects were different in different solid tumors, these studies all demonstrated the important role of CXCR3 in cancer development. To better understand these results, we further summarized the role of CXCR3 in previous studies (25, 35). CXCR3 expression was obviously associated with immune infiltration. The overexpression of CXCR3 might decrease proportion of Th2 cells and IL-4 level, reducing M2 macrophage infiltration and increased the dendritic cell and tumor-infiltrating lymphocyte. We also found the CXCR3 had two different variants which played distinct biological functions (38, 39). The tumor-driven changes influenced the expression of the CXCR3 variants and their ligands promote cancer progression, which were worthy of further investigation.

Our research had some limitations. First, there are more than 40 immune components in the tumor microenvironment. In addition, the immunoregulatory mechanisms could not be fully explained in our study.

Conclusions

In conclusion, our study proposed and verified three novel immune classes (the Immunity-H group, the Immunity-M group, and the Immunity-L group) of HNSCC according to 40 immune components in the tumor environment. Patients in the Immunity-H group responded to monotherapy of radiotherapy, chemotherapy or immunotherapy, and combination radiotherapy and chemotherapy or immunotherapy benefited patients in the Immunity-M group. Moreover, CXCR3 was found to play a significant role in immunoregulation, and CXCR3-targeted therapy may be an ideal candidate treatment to improve the prognosis of patients in the Immunity-L group.

Funding

This research was supported by the Key Research and Development Plan of Liao Ning province (grant no.2017225023), Wu Jieping Medical Foundation (320.6750.2020-05-8) and CSCO-HENGRUI Oncology Research Foundation (Y-HR2019-0326).

References

• 1

  BrayFFerlayJSoerjomataramISiegelRLTorreLAJemalA. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492

  • CrossRef
  • Google Scholar

• 2

  MarurSForastiereAA. Head and Neck Cancer: Changing Epidemiology, Diagnosis, and Treatment. Mayo Clinic Proc (2008) 83(4):489–501. doi: 10.4065/83.4.489

  • CrossRef
  • Google Scholar

• 3

  HammermanPSHayesDNGrandisJR. Therapeutic insights from genomic studies of head and neck squamous cell carcinomas. Cancer Discov (2015) 5(3):239–44. doi: 10.1158/2159-8290.CD-14-1205

  • CrossRef
  • Google Scholar

• 4

  VermorkenJBTrigoJHittRKoralewskiPDiaz-RubioERollandFet al. Open-label, uncontrolled, multicenter phase II study to evaluate the efficacy and toxicity of cetuximab as a single agent in patients with recurrent and/or metastatic squamous cell carcinoma of the head and neck who failed to respond to platinum-based therapy. J Clin Oncol (2007) 25(16):2171–7. doi: 10.1200/JCO.2006.06.7447

  • CrossRef
  • Google Scholar

• 5

  GavrielatouNDoumasSEconomopoulouPFoukasPGPsyrriA. Biomarkers for immunotherapy response in head and neck cancer. Cancer Treat Rev (2020) 84:101977. doi: 10.1016/j.ctrv.2020.101977

  • CrossRef
  • Google Scholar

• 6

  WeberCEKuoPC. The tumor microenvironment. Surg Oncol (2012) 21(3):172–7. doi: 10.1016/j.suronc.2011.09.001

  • CrossRef
  • Google Scholar

• 7

  WuTDaiY. Tumor microenvironment and therapeutic response. Cancer Lett (2017) 387:61–8. doi: 10.1016/j.canlet.2016.01.043

  • CrossRef
  • Google Scholar

• 8

  ZhangZLiuSZhangBQiaoLZhangYZhangY. T Cell Dysfunction and Exhaustion in Cancer. Front Cell Dev Biol (2020) 8:17. doi: 10.3389/fcell.2020.00017

  • CrossRef
  • Google Scholar

• 9

  BindeaGMlecnikBTosoliniMKirilovskyAWaldnerMObenaufACet al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity (2013) 39(4):782–95. doi: 10.1016/j.immuni.2013.10.003

  • CrossRef
  • Google Scholar

• 10

  BhattacharyaSAndorfSGomesLDunnPSchaeferHPontiusJet al. ImmPort: disseminating data to the public for the future of immunology. Immunol Res (2014) 58(2–3):234–9. doi: 10.1007/s12026-014-8516-1

  • CrossRef
  • Google Scholar

• 11

  HänzelmannSCasteloRGuinneyJ. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

  • CrossRef
  • Google Scholar

• 12

  YoshiharaKShahmoradgoliMMartínezEVegesnaRKimHTorres-GarciaWet al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

  • CrossRef
  • Google Scholar

• 13

  LiTFuJZengZCohenDLiJChenQet al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res (2020) 48(W1):W509–14. doi: 10.1093/nar/gkaa407

  • CrossRef
  • Google Scholar

• 14

  SturmGFinotelloFPetitprezFZhangJDBaumbachJFridmanWHet al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics (2019) 35(14):i436–45. doi: 10.1093/bioinformatics/btz363

  • CrossRef
  • Google Scholar

• 15

  ThorssonVGibbsDLBrownSDWolfDBortoneDSYangT-HOet al. The Immune Landscape of Cancer. Immunity (2018) 48(4):812–30.e14. doi: 10.1016/j.immuni.2018.03.023

  • CrossRef
  • Google Scholar

• 16

  SchumacherTNKesmirCvan BuurenMM. Biomarkers in Cancer Immunotherapy. Cancer Cell (2015) 27(1):12–4. doi: 10.1016/j.ccell.2014.12.004

  • CrossRef
  • Google Scholar

• 17

  KimHSKimSCKimSJParkCHJeungH-CKimYBet al. Identification of a radiosensitivity signature using integrative metaanalysis of published microarray data for NCI-60 cancer cells. BMC Genomics (2012) 13:348. doi: 10.1186/1471-2164-13-348

  • CrossRef
  • Google Scholar

• 18

  CharoentongPFinotelloFAngelovaMMayerCEfremovaMRiederDet al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

  • CrossRef
  • Google Scholar

• 19

  WangYWangZXuJLiJLiSZhangMet al. Systematic identification of non-coding pharmacogenomic landscape in cancer. Nat Commun (2018) 9(1):3192. doi: 10.1038/s41467-018-05495-9

  • CrossRef
  • Google Scholar

• 20

  RobinsonMDMcCarthyDJSmythGK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (2010) 26(1):139–40. doi: 10.1093/bioinformatics/btp616

  • CrossRef
  • Google Scholar

• 21

  SzklarczykDFranceschiniAWyderSForslundKHellerDHuerta-CepasJet al. STRING v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res (2015) 43(Database issue):D447–52. doi: 10.1093/nar/gku1003

  • CrossRef
  • Google Scholar

• 22

  ChinC-HChenS-HWuH-HHoC-WKoM-TLinC-Y. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol (2014) 8(Suppl 4):S11. doi: 10.1186/1752-0509-8-S4-S11

  • CrossRef
  • Google Scholar

• 23

  ShannonPMarkielAOzierOBaligaNSWangJTRamageDet al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

  • CrossRef
  • Google Scholar

• 24

  YuGWangL-GHanYHeQ-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 25

  HuMLiKMaskeyNXuZYuFPengCet al. Overexpression of the chemokine receptor CXCR3 and its correlation with favorable prognosis in gastric cancer. Hum Pathol (2015) 46(12):1872–80. doi: 10.1016/j.humpath.2015.08.004

  • CrossRef
  • Google Scholar

• 26

  TangZLiCKangBGaoGLiCZhangZ. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res (2017) 45(W1):W98–102. doi: 10.1093/nar/gkx247

  • CrossRef
  • Google Scholar

• 27

  LiTFanJWangBTraughNChenQLiuJS. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. - PubMed - NCBICancer Res. 77(21):e108–10. doi: 10.1158/0008-5472.CAN-17-0307

  • CrossRef
  • Google Scholar

• 28

  SeligerB. Molecular mechanisms of HLA class I-mediated immune evasion of human tumors and their role in resistance to immunotherapies. HLA (2016) 88(5):213–20. doi: 10.1111/tan.12898

  • CrossRef
  • Google Scholar

• 29

  HalpertMMKonduriVLiangDVazquez-PerezJHofferekCJWeldonSAet al. MHC class I and II peptide homology regulates the cellular immune response. FASEB J (2020) 34(6):8082–101. doi: 10.1096/fj.201903002R

  • CrossRef
  • Google Scholar

• 30

  PeltanovaBRaudenskaMMasarikM. Effect of tumor microenvironment on pathogenesis of the head and neck squamous cell carcinoma: a systematic review. Mol Cancer (2019) 18(1):63. doi: 10.1186/s12943-019-0983-5

  • CrossRef
  • Google Scholar

• 31

  TakeuchiYNishikawaH. Roles of regulatory T cells in cancer immunity. Int Immunol (2016) 28(8):401–9. doi: 10.1093/intimm/dxw025

  • CrossRef
  • Google Scholar

• 32

  PatelSFuSMastioJDominguezGAPurohitAKossenkovAet al. Unique pattern of neutrophil migration and function during tumor progression. Nat Immunol (2018) 19(11):1236–47. doi: 10.1038/s41590-018-0229-5

  • CrossRef
  • Google Scholar

• 33

  FridlenderZGSunJKimSKapoorVChengGLingLet al. Polarization of tumor-associated neutrophil phenotype by TGF-beta: ‘N1’ versus ‘N2’ TAN. Cancer Cell (2009) 16(3):183–94. doi: 10.1016/j.ccr.2009.06.017

  • CrossRef
  • Google Scholar

• 34

  PhilipMFairchildLSunLHorsteELCamaraSShakibaMet al. Chromatin states define tumour-specific T cell dysfunction and reprogramming. Nature (2017) 545(7655):452–6. doi: 10.1038/nature22367

  • CrossRef
  • Google Scholar

• 35

  ChenFYinSNiuLLuoJWangBXuZet al. Expression of the Chemokine Receptor CXCR3 Correlates with Dendritic Cell Recruitment and Prognosis in Gastric Cancer. Genet Test Mol Biomarkers (2018) 22(1):35–42. doi: 10.1089/gtmb.2017.0125

  • CrossRef
  • Google Scholar

• 36

  RenYKanYZKongLF. [Study on the effects of target-silencing CXCR3 expression on malignant proliferation of hepatocellular carcinoma]. Zhonghua Gan Zang Bing Za Zhi (2018) 26(7):508–12. doi: 10.3760/cma.j.issn.1007-3418.2018.07.006

  • CrossRef
  • Google Scholar

• 37

  ZhangYXuLPengM. CXCR3 is a prognostic marker and a potential target for patients with solid tumors: a meta-analysis. Onco Targets Ther (2018) 11:1045–54. doi: 10.2147/OTT.S157421

  • CrossRef
  • Google Scholar

• 38

  WuQDhirRWellsA. Altered CXCR3 isoform expression regulates prostate cancer cell migration and invasion. Mol Cancer (2012) 11(1):3. doi: 10.1186/1476-4598-11-3

  • CrossRef
  • Google Scholar

• 39

  ReyndersNAbboudDBaragliANomanMZRogisterBNiclouSPet al. The Distinct Roles of CXCR3 Variants and Their Ligands in the Tumor Microenvironment. Cells (2019) 8(6):613. doi: 10.3390/cells8060613

  • CrossRef
  • Google Scholar