Patterns of Immune Infiltration in HNC and Their Clinical Implications: A Gene Expression-Based Study

Background: Immune infiltration of head and neck cancer (HNC) highly correlated with the patient's prognosis. However, previous studies failed to explain the diversity of different cell types that make up the function of the immune response system. The aim of the study was to uncover the differences in immune phenotypes of the tumor microenvironment (TME) between HNC adjacent tumor tissues and tumor tissues using CIBERSORT method and explore their therapeutic implications.

Method: In current work, we employed the CIBERSORT method to evaluate the relative proportions of immune cell profiling in 11 paired HNC and adjacent samples, and analyzed the correlation between immune cell infiltration and clinical information. The tumor-infiltrating immune cells of TCGA HNC cohort was analyzed for the first time. The fractions of LM22 immune cells were imputed to determine the correlation between each immune cell subpopulation and survival and response to chemotherapy. Three types of molecular classification were identified via “CancerSubtypes” R-package. The functional enrichment was analyzed in each subtype.

Results: The profiles of immune infiltration in TCGA HNC cohort significantly vary between paired cancer and para-cancerous tissue and the variation could reflect the individual difference. Total Macrophage, Macrophages M0 and NK cells resting were elevated in HNC tissues, while total T cells, total B cells, T cells CD8, B cell navie, T cell follicular helper, NK cells activated, Monocyte and Mast cells resting were decreased when compared to paracancerous tissues. Among each cell immune subtype, T cells regulatory Tregs, B cells naïve, T cells follicular helper, and T cells CD4 memory activated was significantly associated with HNC survival. Three clusters were observed via Cancer Subtypes R-package. Each cancer subtype has a specific molecular classification and subtype-specific immune cell characterization.

Conclusions: Our data suggest a difference in immune response may be an important driver of HNC progression and response to treatment. The deconvolution algorithm of gene expression microarray data by CIBERSOFT provides useful information about the immune cell composition of HNC patients.

Introduction

Head and neck cancer (HNC) remain the primary cause of cancer-related mortality in the world. Oral squamous cell carcinoma (OSCC) is the most important type of HNC, accounting for about 90%, characterized by advanced diagnosis (1), poor prognosis (2), low overall survival (3), and high recurrence (4). Head and neck cancer (HNC) are regarded as the tenth most common cancer in the world and the seventh most common cause of cancer deaths. There are ~400,000 oral and pharyngeal diseases, 160,000 laryngeal cancers and 300,000 deaths worldwide each year (5–7). Head and neck cancer is a common heterogeneous tumor that exists in the oral, pharyngeal, and larynx lesions (8). Contributing factors may include fewer mutant oncogenes, the presence of extensive p53 and lower levels of tumor hypoxia, but the mechanisms remain to be explored in detail. The development of cancer genomics has been widely used to reveal subgroups of patients with different prognoses and different responses to treatment. However, the phenotype of cancer is defined not only by the intrinsic activity of tumor cells but also by immune cells recruited and attracted in the tumor-associated microenvironment. Until now, the roles of immune cells in the tumor-related microenvironment during tumor development remain unclear.

There is increasing evidence that tumor cell immune cells (TIICs), which including B cells, T cells, dendritic cells, macrophages, neutrophils, monocyte, and mast cells, can regulate cancer progression and stimulate attractive therapeutic targets. In addition, the effect of TIICs on cancer development has been largely documented (9–12). In the work, the CIBERSOFT method was applied to define the 22 TICCs subsets of the immune response in HNC in order to examine the mutual relationship with molecular subsets and clinical characteristics. Our findings also revealed distinct immune phenotypes for molecular HNC subclasses. Moreover, the present study provides a novel insight into immunotherapy strategy for HNC.

Materials and Methods

Data Source and Data Processing

Firstly, the RNA-seq data (FPKM values) of HNC cohort were downloaded from the Cancer Genome Atlas (TCGA) website (https://cancergenome.nih.gov/), including four files (gdc_download_20190201_030206.452249.tar.gz, gdc_manifest_20190201_025835.txt, metadata.cart.2019-02-01.json and clinical.cart.2019-02-01.json). Secondly, we merged the data using Perl software (mRNA_merge.pl). RNA-seq profiles were normalized using Voom standardized method (variance modeling at the observational level) (13).

Quantification of TIICs Using the CIBERSOFT Algorithm

The CIBERSOFT method is a gene-based deconvolution algorithm that infers 22 human immune cell types and uses the characteristics of 547 marker genes to quantify the relative scores for each cell type. To enhance the robustness of the results, CIBERSOFT algorithm uses Monte Carlo sampling to derive the deconvoluted P value for each sample. The standardized processed dataset of gene expression is uploaded to the CIBERSOFT website (https://cibersort.stanford.edu/index.php), which runs using 1000 aligned default signature matrices (14).

Total T cells were calculated as a sum of CD8+ T cells, CD4+ naïve T cells, CD4+ memory resting T cells, CD4+ memory activated T cells, follicular helper T cells, regulatory T cells (Tregs), and T cells gamma delta fractions. Total macrophage fraction was imputed as a sum of M0, M1, and M2 macrophage fractions. And total B cells was estimated as a sum of B cells memory and B cells naïve.

Survival Analysis

We conducted univariate Cox analysis and Kaplan-Meier survival analysis between LM22 immune cell subsets and OS and screened the prognostic 22 human immune cell phenotypes. The univariate Cox analysis and Kaplan-Meier survival analysis was performed by the survfit function in “survival” package in R software. Then, the multivariate Cox regression analysis was used to further validate prognostic 22 human immune cell phenotypes as prognosis factors. We also explore the association between LM22 immune cell and other clinical information, such as TNM stage, Radiation, Grade, Age, and HPV status.

Identification of Molecular Subclassification

In order to further to explore different TME cell infiltration patterns, a consensus cluster algorithm was applied to determine the number of clusters. The process was conducted using “CancerSubtypes” R-package (15). To uncover the potential difference in TME cell infiltration models among the clusters, differentially expressed genes (DEGs) and different immune cell types were identified using unpaired Student t tests. The data set with |log2 fold change| ≥ 0.2 and P –value less than 0.05 was considered selection criteria for subsequent analysis.

Functional and Pathway Enrichment Analysis

To uncover the potential biological significance of DEGs among TME subtypes, Gene Ontology (GO) Biological Process term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were conducted using “ClusterProfiler” R package (16). GO enrichment analysis was based on the threshold of P-value < 0.05 and false discovery rate (FDR) of < 0.05. Gene Set Variation Analysis (GSVA) was conducted to unveil an overall pathway of gene-set activity score for each sample (17). The Gene sets using the c2 curated signatures were downloaded from the Molecular Signature Database (MSigDB) of Broad Institute. The significant enrichment pathway was identified based on the |logFoldChange| ≥0.2 and adjust P-value < 0.05.

Statistical Analysis

Patients with a CIBERSOFT P value less than 0.05 were included in the main study. After the initial screen, eligible samples were divided into two groups: HNC tumor samples and adjacent samples. Immune cell profiles were estimated for each sample. In order to explore the association among 21 immune cell subtypes in three groups, Pearson correlation analysis was calculated using R software. The wilcox.test was employed to test the statistical significance between the two groups. For comparisons of more than two groups, Kruskal-Wallis tests were employed to examine the difference among groups (18). Cox regression analysis was applied to determine the association between inferred proportions of immune cell types and survival. The log-rank statistic was used to test the differences in OS survival between groups using Kaplan-Meier plots. In the univariate Cox regression analysis, variables with a P < 0.05 were considered as independent prognostic overall survival (OS) factors, and the included prognostic factors were used to build the multivariate Cox regression model for OS. Clinical variables, such as age, sex, HPV status, lymph node metastasis, distant metastasis, grade, and TNM stage, were included in the multivariate Cox regression model. To evaluate the relationship between different immune cell subtypes and response to radiation, the wilcox.test was conducted.

A heatmap was produced with the R package “ComplexHeatmap” (19). The R package “pROC” was used to visualize operating characteristic (ROC) curves to impute the area under the curve (AUC) and confidence intervals to evaluate the diagnostic accuracy of LM 22 immune cell (20). Statistical analysis was performed using R-Language (R-project.org) and packages obtained through the Bioconductor project (www.bioconductor.org). All P values were bilateral and a P value of < 0.05 was considered statistically significant.

Results

Overview of Data

A total of 546 samples, included 44 adjacent samples, and 502 tumor samples, were obtained from the TCGA. After performing CIBERSOFT algorithm, 454 patients (11 normal patients and 443 tumor patients) with a P value < 0.05 was considered in the study, including 41 paracancerous tissue, and 11 paired tumor tissue. Meanwhile, 547 TME corresponding gene expression profiles were also filtered for further analysis.

Profile of TME in HNC and Clinicopathological Characteristics of TME Subtypes

The landscape of TME cell infiltration models and MTE signatures was systematically evaluated by CIBERSOFT algorithm. Figure 1 summarizes the findings obtain from the 41 paired tumor samples and 11 paired adjacent samples. Obviously, the proportions of TME cells in HNC varies significantly between both intra- and intergroup. The fraction of total T cells and total B cells were higher in HNC adjacent tissue than tumor tissue. Total Macrophage was mainly observed in the adjacent tissue (Figure 2). We found that T cell follicular helper, T cells CD8, NK cells activated, Monocyte, Macrophages M0, NK cells resting, Mast cells resting, and B cell navie were significantly changed between normal and tumor group, while B cells memory, Plasma cells, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells regulatory (Tregs), T cells gamma delta, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated Eosinophils, and Neutrophils were not obviously altered between groups (Figure 3, Table 1).

FIGURE 1

Figure 1. The performance of CIBERSOST for estimating TIICs composition in HNC. (A) Relative proportions of 22 TIICs subpopulation in normal samples. (B) Relative proportions of 22 TIICs subpopulation in tumor samples. (C) Heat map of the 22 immune cell proportions.

FIGURE 2

Figure 2. The stacked histogram shows the distribution of 22 immune cell infiltration between HNC adjacent tissues and tumor tissues, including total immune cells (A), total T cells (B), total B cells (C), and total macrophages (D).

FIGURE 3

Figure 3. The Violin plot exhibits the difference between CIBERSOFT immune cell fractions between HNC adjacent tissues and tumor tissues.

TABLE 1

Table 1. Comparison of CIBERSORT immune fractions between HNC tumor tissues and adjacent tissues.

Form the Figure 4, we can see that several highly positive relationship among LM22 immune cells in the HNC paired adjacent tumor samples were found, while mutual relationship among LM22 immune cells reduced in the tumor samples, such as Macrophages M2 was highly positive associated with Monocytes in the immune phenotype profiles in the TMC form the HNC paired adjacent tumor samples, while T cells regulatory (Tregs) was highly negatively related to T cells CD8 in the HNC tumor group. Thus, we hypothesized that alterations in the TME cell infiltration ratio directly reflect the difference in immunity between the two groups.

FIGURE 4

Figure 4. Correlation matrix of all 22 immune proportions and immune cytolytic activity in the TCGA HNC cohort, including total samples (A), normal samples (B) and tumor samples (C).

In terms of clinical features, we can see that the fraction of T cells CD8, T cells follicular helper, T cells regulatory (Tregs) and T cells CD4 memory activated was higher in HNC HPV positive tissue than HNC HPV negative tissue, while NK cells resting, Macrophages M2, T cells CD4 memory resting, Macrophages M0 and Neutrophils were higher in the HNC HPV negative samples (Figure 5, Table 2). T cells follicular helper was mainly observed in the early stage (G1/G2) of HNC (Supplementary Figure 1, Supplementary Table 1). Neutrophils and Macrophages M1 was associated with radiation therapy (Supplementary Figure 2, Supplementary Table 2). Therefore, these results demonstrated that aberrant immune infiltration and it's heterogeneous in HNC as a tightly regulated process might play important roles in the development of tumor, it has important clinical meanings.

FIGURE 5

Figure 5. The difference between immune infiltration between HNC HPV positive and negative samples. (A) The total difference in immune cells infiltration. Box plot of the distribution of CIBERSORT P value for Macrophages M0 (B), Macrophages M2 (C), Neutrophils (D), NK cells resting (E), T cells CD4 memory activated (F), T cells CD4 memory resting (G), T cells CD8 (H), T cells follicular helper (I) and T cells regulatory (Tregs) (J), respectively.

TABLE 2

Table 2. Comparison of CIBERSORT immune fractions between HNC HPV positive patients and negative patients.

Identification of Prognostic LM22 Immune Cell Subtypes

The univariate Cox regression was conducted to screen the prognostic LM22 immune cell subsets in all tumor samples, and the results were shown that a total of 6 immune cell subsets (Eosinophils, T cells regulatory Tregs, Mast cells activated, B cells naïve and T cells follicular helper) was significantly correlated with OS when the P-value was less than 0.05 (Figure 6A). B cells naive, T cells regulatory (Tregs), T cells follicular helper were associated with improved outcome. Mast cells activated and Eosinophils were associated with poorer outcome. Then, we conducted a Kaplan-Meier curve plot and log-rank test for the above-mentioned immune cell subsets, and the results are shown in Figure 7. Eosinophils, Mast cells activated, B cells naïve and T cells follicular helper were significantly associated with OS with HNC patients.

FIGURE 6

Figure 6. The prognostic associations of subsets of immune cells in Cox regression analysis. (A) Univariate Cox regression analysis. (B) Multivariate Cox regression analysis.

FIGURE 7

Figure 7. The K-M survival analysis of four immune cells, including B cells naïve (A), Eosinophils (B), Mast cells activated (C), and T cells follicular helper (D).

Next, multivariate Cox proportional hazards regression analysis was further identified the prognostic LM22 immune cell subsets. The results indicated that T cells regulatory (Tregs) was associated with improved outcome (Figure 6B). The ROC curves were used to evaluate the prognostic power of prognostic LM22 immune cell subsets. The AUC for prognostic LM22 immune cell subsets biomarkers prognostic model was shown in Figure 8. In HNC, Dendritic cells activated, B cells naïve, Dendritic cells resting, Macrophages M0, Macrophages M2, Plasma cells, and T cells follicular helper had an AUC value of more than 0.550, and T cells follicular helper had the highest performance in the risk prediction of HNC patients. Detailed results are provided in Table 3.

FIGURE 8

Figure 8. ROC curves of prognostic predictors built by the 22 immune cell subsets.

TABLE 3

Table 3. The AUC value for each immune cell.

Subsequently, the association between clinical features and LM22 immune cell subsets was performed to confirm the mutual relationship. After analysis, clinical covariates of HPV status, and gender, grade, race, radiation, age were significantly correlated with immune cell subclasses. However, clinical covariates of alcohol consumption, tobacco consumption, TNM stage, and neoadjuvant treatment are not significantly correlated with LM22 immune cell subsets (Table 4, Figure 9). Among these clinical variables, radiation/grade/HPV was significantly associated with immune cell subclasses (Figure 10).

TABLE 4

Table 4. Differences in pathway activities scored per cell by GSVA between HNC tumor and normal tissues.

FIGURE 9

Figure 9. The correlation heat map of 22 immune cells and clinical features. The value in the space shows the correlation coefficient (A) and the P value (B).

FIGURE 10

Figure 10. The forest plots show the correlation between clinical characteristics and immune cell subsets. (A) Is Grade, (B) is for HPV status, and (C) is for Radiation.

Immune Cell Patterns in Molecular HNC Subclasses

In order to uncover different immune cell population infiltration of TME in HNC, molecular classification of HNC was identified by performing unsupervised consensus clustering of all tumor samples. The optimal number of clusters was determined by the K value. After assessing the relative change in the area under the cumulative distribution function (CDF) curve and consensus heatmap, we selected a three-cluster solution (K = 3), which has no appreciable increase in the area under the CDF curve (Figure 11). This finding resulted in 252 patients (59%) in cluster I, 38 patients (9%) in cluster II, and 138 patients (32%) in cluster III for the HNC cohort. The consensus matrix heatmap revealed the Cluster I and Cluster II, Cluster III well appeared individualized clusters in Figure 12A. The sample of each cluster was shown in Figure 12B. Moreover, clusters were associated with distinct patterns of survival in Figure 12C. Compared with Cluster I and Cluster III, the patients who were classified as Cluster II had a good prognosis (P < 0.001, log-rank test).

FIGURE 11

Figure 11. The cluster counts evaluated. (A) CDF curve of K = 2–7. (B) The relative change in area under the CDF curve of K = 2–7.

FIGURE 12

Figure 12. The cancer subtypes using SNFCC+ algorithm. (A) Log-rank test p-value for Kaplan–Meier survival analysis, (B) clustering heatmap visualizing the degree of the partitioning of the sample clusters, (C) average silhouette width representing the coherence of clusters.

Differentially Expressed Analysis of Genes/LM22 Immune Cell Fraction in Each HNC Subgroup

The Kruskal-Wallis test was used to identify quantitative genes/LM22 immune cell significantly associated with each subtype. As for differential infiltration of LM22 immune cell in HNC TME, Cluster 1 was defined by high level of Macrophages M1 and Macrophages M2, Cluster 2 was enriched by B cells memory, B cells naive, Monocytes, Plasma cells, T cells CD8, T cells gamma delta, T cells regulatory (Tregs), NK cells activated and T cells follicular helper. While Cluster 3 was defined by the high level of Macrophages M0, Mast cells activated, Neutrophils, NK cells resting, and T cells CD4 memory resting (Figure 13).

FIGURE 13

Figure 13. Box plot depicting the association between immune cell subset and three clusters, depicted P-values are from the Kruskal-Wallis tests. (A–O) Is for B cells memory, B cells naïve, Macrophages M0, Macrophages M1, Macrophages M2, Mast cells activated, Monocytes, Neutrophils, NK cells activated, NK cells resting, Plasma cells, T cells CD4 memory resting, T cells CD8, T cells gamma delta, and T cells regulatory (Tregs), respectively.

In order to seek the molecular differences between HNC molecular subtypes and deriving subtype-specific biomarkers, the unpaired student t-test was used to identify quantitative genes significantly associated with each subtype. Gene differentially expressed analysis for unmatched subgroups was executed with the threshold of absolute log-fold change cutoff greater than 0.2 and FDR = 0.05. Figure 14 was shown by DEGs in concentric circles radiating among 3 Clusters. In subgroup 1 compared to subgroups 2, a total of 227 mRNAs (192 upregulated genes and 35 downregulated genes) were differentially expressed genes; In subgroup 1 compared to subgroups 3, a total of 242 differentially expressed mRNAs (14 upregulated genes and 228 down-regulated genes) were detected; In subgroup 2 compared to subgroups 3, a total of 333 differentially expressed mRNAs (41 upregulated genes and 292 down-regulated genes) were observed.

FIGURE 14

Figure 14. Circular visualization of chromosomal positions, hierarchical clustering, connectivity and severity gradients (Left) and a volcano plot (Right) in HNC subgroups. (A–C) Is for cluster I vs. Cluster II, cluster I vs. cluster III and cluster II vs. Cluster III.

Gene Functional Annotation and Gene Ontology (GO) Enrichment Analysis, GSVA Analysis for DEGs to Identify Molecular Subtypes

In subgroup 1 compared to subgroup 2, a total of 533 GO terms of biological process, 8 GO terms of cellular component, and 53 GO terms of molecular function were identified to be significant (adjust P-value < 0.05). In subgroup 1 compared to subgroups 3, A total of 692 GO terms of biological process, 19 GO terms of cellular component, and 66 GO terms of molecular function were identified to be significant. In subgroup 2 compared to subgroups 3, A total of 676 GO terms of biological process, 19 GO terms of cellular component, and 55 GO terms of molecular function were identified to be significant. Top GO terms included cytokine activity, immune/inflammatory response, and chemokine activity (Figure 14). In addition, all the pathways that were yielded from the KEGG analysis were associated with immune response (Figure 15).

FIGURE 15

Figure 15. The GO and KEGG analysis for three HNC clusters. (A–F) Is for cluster I vs. Cluster II, cluster I vs. cluster III and cluster II vs. Cluster III.

Gene set variation analysis with three clusters was analyzed by GSVA package of R software. The number of enriched pathways progressively increased from Subtype 1 through Subtype 3. The most significant gene sets enriched were ordered by significance (false discovery rate (FDR) P- and adjusted P-values) and listed in Table 4. As shown in Figure 16, several hallmark gene-sets, including “TNFA_SIGNALING_VIA_NFKB”, “APICAL_JUNCTION”, “EPITHELIAL_MESENCHYMAL_TRANSITION (EMT),” and “KRAS_SIGNALING_UP”, “ANGIOGENESIS”, were observed.

FIGURE 16

Figure 16. The GSVA analysis for three HNC clusters.

Discussion

According to statistics from the National Central Cancer Registry of China (21) and the National Cancer Institute of USA (22), HSC is still the primary cancer disease worldwide. Especially in China, HNC is in the top ten causes of incidence rates and mortality in urban and rural regions in China (21). HNC is known for its rapid clinical progression and poor prognosis. The overall survival rate of HNC patients has not improved significantly over the past 20 years with the 5- year survival rate is ~45–50% (23). Therefore, HNC remains high-risk cancer affecting large numbers of people, threatening the lives of patients through local recurrence and metastatic cervical lymph nodes. Meanwhile, the potential molecular mechanism of HNC is not unknown and, more efforts should be used to this area to uncover the corresponding targeted therapy.

In current work, the CIBERSOFT was applied to evaluate the differential immune cell infiltration in paired HNC and adjacent normal tissue, and the results revealed that considerable difference in immune cell fraction occurred in both intra- and intergroup. Our work also uncovers the detail of infiltration of LM22 immune cell subsets in HNC that the proportions of macrophages account for more than 45%, in which 27% is M0, 9% is M2, and M1 make up 9%. In addition, T cells CD4 memory resting as the second proportion (12%) (Table 1). Compared to HNC and adjacent normal tissues, the proportions of Total Macrophage, Macrophages M0 and NK cells resting, total T cells, total B cells, T cells CD8, B cell navie, T cell follicular helper, NK cells activated, Monocyte, Mast cells resting produced statistical significance (P < 0.05). The total T cells, total B cells, T cells CD8, B cell navie, T cell follicular helper, NK cells activated, Monocyte, and Mast cells resting were increased in adjacent normal tissues when compared to the HNC tissues. In contrast, Macrophages M0 and NK cells resting were decreased in adjacent normal tissues. The mechanisms behind the activation of NK cells resting and Macrophages M0, deactivation of B cells naïve, T cells gamma delta, NK cells activated, and Mast cells resting in HNC remain unclear. Macrophages cells and NK cells have been detected in HNC (24, 25), we hypothesize that tumor cell-derived metabolites such as oxidized natural polyamines might be responsible for activation of Macrophages and NK cells and deactivation of Mast cells.

The prognostic importance of immune cell infiltration has been identified for various solid tumor types (26–28). In univariate Cox analysis, we found that B cells naive, T cells regulatory (Tregs), T cells follicular helper were associated with improved outcome, while Mast cells activated and Eosinophils were associated with poorer outcome. Eosinophils and Mast cells activated were associated with poorer outcome. To the best of our knowledge, the T- cell immune response is the primary event in antitumor immunity (29, 30), especially in HNC (31). The T and B cells are present in immune cell infiltrates of HNC and the degree of tumor-infiltrating T and B cells correlated with improved survival of cancers (32). Nevertheless, Eosinophils has been recognized to correlate with angiogenesis and metastasis (33) and increased Eosinophils proportions are related to poor prognosis. This immune phenotype was in line with a previous study (34). In addition, several immune cells are associated with clinical outcome. T cells follicular helper correlated with an increase in the level of GRADE stage. Total T cells and Macrophages cells, Neutrophils are associated with HPV status. Neutrophils and Macrophages M1 are related to radiation therapy. Our data further validated these findings from the previous studies that specific immune cells were related to predicting clinical outcome (35–37).

HNC can be robustly divided into three subtypes by SNFCC+ method. And the results demonstrated that three subtype classifications were significantly associated with patients' survival. Compared with Cluster I and Cluster III, the patients who were classified as Cluster II had a good prognosis. Compare with each Cluster of LM22 immune cell in HNC TME, Cluster I was defined by a high level of Macrophages cells infiltration, Cluster II was enriched by B cells infiltration and T cells infiltration. While Cluster III was defined by the high level of Mast cells infiltration, Neutrophils infiltration, and NK cells infiltration. Compared with each Cluster of DEGs, each cluster has its characteristic functional enrichment terms. Therefore, the three subtypes differed in overall survival, molecular characteristics.

In summary, our analysis of LM22 immune cell subsets using by CIBERSORT deconvolution algorithm provides the full information on immune cell landscape of HNC. Our finding has also uncovered a vital role in prediction for clinical outcome. In current work, this comprehensive assessment of LM22 immune cell infiltration models in TME shed light on how tumors respond to immunotherapy and might help clinicians to explore the development of the new drug.

Data Availability Statement

The datasets analyzed in this study are available in The Cancer Genome Atlas (TCGA) public repository (https://cancergenome.nih.gov/).

Author Contributions

JSo, DY, JZ, and ZD: wrote the main manuscript text. JSo and JL prepared Figures 1–16. JSu and JZ contributed to data analysis. All authors reviewed the manuscript.

Funding

This study was supported by the NSFC (Natural Science Foundation of China) (Grant No. 81873608), Mechanism of PGE2 activating ureteral smooth muscle cAMP signaling pathway involved in ureteral drainage; The Science and Technology Plan Project of Guizhou Province in 2019 (Grant No. [2019]2799), Miao Yao Ning Tetai promotes stone discharge by relaxing ureteral smooth muscle with PGE2/cAMP signaling pathway; The High-level innovative talent project of Guizhou Province in 2018 (Grant No. [2018]5639), Guizhou high-level innovative talent project (100 levels) and The Science and Technology Plan Project of Guiyang in 2019 (Grant No. [2019]2-15), and Application of Artificial Intelligence in the Background of Big Data in Accurate Diagnosis and Treatment of Prostate Cancer.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2019.01285/full#supplementary-material

Supplementary Figure 1. The difference of immune infiltration between HNC the early stage (G1/G2) and the late stage (G3/G4) samples. (A) The total difference in immune cells infiltration. (B) Box plot of the distribution of CIBERSORT P value for T cells follicular helper.

Supplementary Figure 2. The difference of immune infiltration between HNC with radiation therapy and without radiation therapy samples. (A) The total difference in immune cells infiltration. (B,C) Box plot of the distribution of CIBERSORT P value for Macrophages M2 and Neutrophils, respectively.

Supplementary Table 1. Comparison of CIBERSORT immune fractions between grade G1/G2 and G3/G4.

Supplementary Table 2. Comparison of CIBERSORT immune fractions between with and without radiation therapy.

References

1. Preciado DA, Matas A, Adams GL. Squamous cell carcinoma of the head and neck in solid organ transplant recipients. Head Neck. (2002) 24:319–25. doi: 10.1002/hed.10055

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ordonez R, Otero A, Jerez I, Medina JA, Lupianez-Perez Y, Gomez-Millan J. Role of radiotherapy in the treatment of metastatic head and neck cancer. Onco Targets Ther. (2019) 12:677–83. doi: 10.2147/OTT.S181697

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Patel BP, Shah PM, Rawal UM, Desai AA, Shah SV, Rawal RM, et al. Activation of MMP-2 and MMP-9 in patients with oral squamous cell carcinoma. J Surg Oncol. (2005) 90:81–8. doi: 10.1002/jso.20240

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Choi KK, Kim MJ, Yun PY, Lee JH, Moon HS, Lee TR, et al. Independent prognostic factors of 861 cases of oral squamous cell carcinoma in Korean adults. Oral Oncol. (2006) 42:208–17. doi: 10.1016/j.oraloncology.2005.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Folz BJ, Silver CE, Rinaldo A, Fagan JJ, Pratt LW, Weir N, et al. An outline of the history of head and neck oncology. Oral Oncol. (2008) 44:2–9. doi: 10.1016/j.oraloncology.2007.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Mehanna H, Paleri V, West CM, Nutting C. Head and neck cancer–Part 1: epidemiology, presentation, and prevention. BMJ. (2010) 341:c4684. doi: 10.1136/bmj.c4684

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin. (2018) 68:7–30. doi: 10.3322/caac.21442

CrossRef Full Text | Google Scholar

8. Arantes LM, de Carvalho AC, Melendez ME, Carvalho AL, Goloni-Bertollo EM. Methylation as a biomarker for head and neck cancer. Oral Oncol. (2014) 50:587–92. doi: 10.1016/j.oraloncology.2014.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Shibutani M, Maeda K, Nagahara H, Fukuoka T, Iseki Y, Matsutani S, et al. Tumor-infiltrating lymphocytes predict the chemotherapeutic outcomes in patients with stage IV colorectal cancer. In Vivo. (2018) 32:151–8. doi: 10.21873/invivo.11218

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Zhang S, Zhang E, Long J, Hu Z, Peng J, Liu L, et al. Immune infiltration in renal cell carcinoma. Cancer Sci. (2019) 110:1564–72. doi: 10.1111/cas.13996

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Aponte-Lopez A, Fuentes-Panana EM, Cortes-Munoz D, Munoz-Cruz S. Mast Cell, the neglected member of the tumor microenvironment: role in breast cancer. J Immunol Res. (2018) 2018:2584243. doi: 10.1155/2018/2584243

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Cai DL, Jin LP. Immune cell population in ovarian tumor microenvironment. J Cancer. (2017) 8:2915–23. doi: 10.7150/jca.20314

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Law CW, Chen Y, Shi W, Smyth GK. voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. (2014) 15:R29. doi: 10.1186/gb-2014-15-2-r29

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xu T, Le TD, Liu L, Su N, Wang R, Sun B, et al. CancerSubtypes: an R/Bioconductor package for molecular cancer subtype identification, validation and visualization. Bioinformatics. (2017) 33:3131–3. doi: 10.1093/bioinformatics/btx378

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. (2013) 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Hazra A, Gogtay N. Biostatistics series module 3: comparing groups: numerical variables. Indian J Dermatol. (2016) 61:251–60. doi: 10.4103/0019-5154.182416

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. (2016) 32:2847–9. doi: 10.1093/bioinformatics/btw313

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinform. (2011) 12:77. doi: 10.1186/1471-2105-12-77

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Pan R, Zhu M, Yu C, Lv J, Guo Y, Bian Z, et al. Kadoorie biobank collaborative: cancer incidence and mortality: a cohort study in China, 2008–2013. Int J Cancer. (2017) 141:1315–23. doi: 10.1002/ijc.30825

CrossRef Full Text | Google Scholar

22. Tealab SH, Sedhom NFH, Hassouna A, Gouda I, Ismail H. Prevalence of human papilloma virus in oropharyngeal, tongue and lip squamous cell carcinoma: an experience from the Egyptian National Cancer Institute. J Investig Med. (2019) 67:1061–6. doi: 10.1136/jim-2018-000968

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Omar EA. The outline of prognosis and new advances in diagnosis of oral squamous cell carcinoma (OSCC): review of the literature. J Oral Oncol. (2013) 2013:519312. doi: 10.1155/2013/519312

CrossRef Full Text | Google Scholar

24. Chen X, Fu E, Lou H, Mao X, Yan B, Tong F, et al. IL-6 induced M1 type macrophage polarization increases radiosensitivity in HPV positive head and neck cancer. Cancer Lett. (2019) 456:69–79. doi: 10.1016/j.canlet.2019.04.032

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Friedman J, Padget M, Lee J, Schlom J, Hodge J, Allen C. Direct and antibody-dependent cell-mediated cytotoxicity of head and neck squamous cell carcinoma cells by high-affinity natural killer cells. Oral Oncol. (2019) 90:38–44. doi: 10.1016/j.oraloncology.2019.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Zhang J, Wang J, Qian Z, Han Y. CCR5 is associated with immune cell infiltration and prognosis of lung cancer. J Thorac Oncol. (2019) 14:e102–3. doi: 10.1016/j.jtho.2018.12.037

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Miksch RC, Schoenberg MB, Weniger M, Bosch F, Ormanns S, Mayer B, et al. Prognostic impact of tumor-infiltrating lymphocytes and neutrophils on survival of patients with upfront resection of pancreatic cancer. Cancers. (2019) 11:39. doi: 10.3390/cancers11010039

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Zhou R, Zhang J, Zeng D, Sun H, Rong X, Shi M, et al. Immune cell infiltration as a biomarker for the diagnosis and prognosis of stage I-III colon cancer. Cancer Immunol Immunother. (2019) 68:433–42. doi: 10.1007/s00262-018-2289-7

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Hermans IF, Ritchie DS, Yang J, Roberts JM, Ronchese F. CD8+ T cell-dependent elimination of dendritic cells in vivo limits the induction of antitumor immunity. J Immunol. (2000) 164:3095–101. doi: 10.4049/jimmunol.164.6.3095

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhang Y, Bush X, Yan B, Chen JA. Gemcitabine nanoparticles promote antitumor immunity against melanoma. Biomaterials. (2019) 189:48–59. doi: 10.1016/j.biomaterials.2018.10.022

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Alhamarneh O, Amarnath SM, Stafford ND, Greenman J. Regulatory T cells: what role do they play in antitumor immunity in patients with head and neck cancer? Head Neck. (2008) 30:251–61. doi: 10.1002/hed.20739

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Prizment AE, Vierkant RA, Smyrk TC, Tillmans LS, Nelson HH, Lynch CF, et al. Cytotoxic T cells and granzyme B associated with improved colorectal cancer survival in a prospective cohort of older women. Cancer Epidemiol Biomarkers Prev. (2017) 26:622–31. doi: 10.1158/1055-9965.EPI-16-0641

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Nissim Ben Efraim H, Levi-Schaffer F. Roles of eosinophils in the modulation of angiogenesis. Chem Immunol Allergy. (2014) 99:138–54. doi: 10.1159/000353251

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Varricchi G, Galdiero MR, Loffredo S, Lucarini V, Marone G, Mattei F, et al. Eosinophils: the unsung heroes in cancer? Oncoimmunology. (2018) 7:e1393134. doi: 10.1080/2162402X.2017.1393134

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Angell HK, Lee J, Kim KM, Kim K, Kim ST, Park SH, et al. PD-L1 and immune infiltrates are differentially expressed in distinct subgroups of gastric cancer. Oncoimmunology. (2019) 8:e1544442. doi: 10.1080/2162402X.2018.1544442

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Falk AT, Yazbeck N, Guibert N, Chamorey E, Paquet A, Ribeyre L, et al. Effect of mutant variants of the KRAS gene on PD-L1 expression and on the immune microenvironment and association with clinical outcome in lung adenocarcinoma patients. Lung Cancer. (2018) 121:70–5. doi: 10.1016/j.lungcan.2018.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Falkenius J, Johansson H, Tuominen R, Frostvik Stolt M, Hansson J, Egyhazi Brage S. Presence of immune cells, low tumor proliferation and wild type BRAF mutation status is associated with a favourable clinical outcome in stage III cutaneous melanoma. BMC Cancer. (2017) 17:584. doi: 10.1186/s12885-017-3577-x

PubMed Abstract | CrossRef Full Text | Google Scholar