A Tumor-Infiltration CD8+ T Cell-Based Gene Signature for Facilitating the Prognosis and Estimation of Immunization Responses in HPV+ Head and Neck Squamous Cell Cancer

Background: CD8+ T cells, which play a vital role in response to adaptive immunity, are closely related to the immunization responses to kill tumor cells. Understanding the effects exerted by tumor-infiltrated CD8+ T cells in HPV+ and HPV- head and neck squamous cell carcinoma (HNSCC) patients is critical for predicting their prognosis as well as their responses towards immunization-related therapy.

Materials and Methods: HNSCC single cell transcriptome was used to screen for differentially expressed genes (DEGs) based on CD8+ T cells. A gene signature associated with CD8+ T cells was built and verified with the cancer genome atlas dataset with a view to predicting the prognosis of HNSCC patients. Risk scores were calculated for HNSCC cases and categorized into either high- or low-risk cohorts. The prognosis-correlated data of the risk scores were analyzed by using Kaplan-Meier survival curves and multi-variate Cox regression plots. In addition, the possibility of using the genetic profiles to predict responses toward immunization-related therapy was explored.

Results: From the DEGs screened from the sequencing of single-cell RNA, a gene signature of 4 genes (ACAP1, ANKRD28, C12orf75, and M6PR) were identified. It was seen that these genes could predict overall survival in HPV+ HNSCC patients. In addition, high- and low-risk HPV+ HNSCC patients showed marked differences in their CD8+ T-cell infiltration due to immunization when clinical characteristics were taken into consideration. This correlated with their immunization therapy responses.

Conclusions: Our work provides insights into explaining the restricted responses of current immunization checkpoint inhibiting substances in HPV+ HNSCC patients. A novel genetic signature to predict the prognosis and immunization-correlated therapeutic responses is presented. This will provide potential new therapeutic opportunities for HPV+ HNSCC patients.

Introduction

Head and neck squamous cell carcinoma (HNSCC) refers to cancer of the oral cavity, oropharynx and larynx, and is the sixth commonest carcinoma globally (1). More than 550,000 new patients are reported each year, accounting for approximately 4% of carcinomas worldwide (2). Because most cases present with locally advanced disease, HNSCC is correlated with poor prognosis and results in high mortality (3). Conventional treatments exhibit limited effectiveness, and new therapeutic strategies capable of broadening the existing treatment options for HNSCC are urgently required (4). Recent clinical trials have demonstrated that programmed death ligand 1 (PD-L1) or programmed death 1 (PD-1) blockade exhibit clinically meaningful anti-tumor activity together with an acceptable safety profile when used in the treatment of HNSCC patients (5, 6). However, despite this progress, only about 20-30% of HNSCC cases survived after anti-PD-1/PD-L1 therapies, and the response towards PD-1/PD-L1 blockade is still far from satisfactory. Therefore, there is an urgent need to further understand the immunization state of the cell during this disease and identify features correlated to the response ability of existing immunization therapies, thereby paving the way to the development of new single and multi-drug immunization therapies.

HNSCC, is a type of carcinoma that can arise due to genetic alteration caused by either exposure to carcinogens (such as alcohol and/or tobacco) or via malignant conversion due to HPV infection (7, 8). There is an alarming growth of HPV + HNSCC in western countries, with up to half of all HNSCC patients residing in the US, as the presence of HPV infections is considered to increase the risk of the disease (9). However, HPV-correlated HNSCC is suggested to exhibit distinctive biological and clinically-associated characteristics, with the presence of HPV conferring a survival advantage when compared to its absence (10). Distinct tumor-infiltration immunization populations were identified in HNSCC patients, with a greater proportion of dysfunctional CD8+ T cells seen in HPV- HNSCC (11). A higher rate of response towards PD-1/PD-L1 was also identified in HPV+ patients compared to those who were HPV-. Thus, the presence of HPV may well be a factor that can be used to classify HNSCC. However, the underlying mechanisms and potential associations between the HPV state and the tumor immunization environment still needs to be characterized.

In the present study, we attempted to elucidate the correlation between the HPV state and immunization environment-related factors by using a multicenter database. A single-cell RNA sequencing dataset was used to assess the various subpopulations of immune cells and particular genes that may differ in HPV+ and HPV- HNSCC patients. Using a combination of RNA-seq data from a large number of HNSCC cases and their corresponding clinical information, a gene signature for tumor-infiltrated CD8+ T cells was established using multiple machine learning algorithms. This risk-associated gene signature was verified using the gene expression profiles and clinically associated information from an independent cancer genome atlas (TCGA) provisional dataset. The genetic signature obtained may provide future targets for increasing our knowledge of the mechanisms that govern HPV+ and HPV- HNSCC. This study may also increase immunization checkpoint blockade therapy efficacy with respect to treatment of this disease.

Materials and Methods

CD8+ T Cell Estimation in HNSCC Patients

The TIMER2.0 database (http://timer.comp-genomics.org/) was utilized to explore the relationship between tumor infiltration of CD8+ T cells and the prognosis of HNSCC in patients (12). We analyzed the immunization infiltration CD8+ T cells in different carcinoma types by multiple immunization deconvolution approaches, by using Cox regression correlation and Kaplan-Meier survival curves. These data obtained were used to correlate the prognosis data of the relevant immunization infiltration data in a range of carcinoma categories.

Research of CD8+ T Cell-Correlated Immunization-Correlated Genes in HNSCC

Single cell transcriptomes include 130,721 cells from HPV- HNSCC and HPV+ HNSCC (10). Both cohorts consist of PBMC and tumor infiltrated leukocytes (TIL). In the current analysis we have only included 60,676 cells from the TIL.

Single-Cell RNA-Seq Data Analysis

Specific to the integrated investigation of single-cell data, these were normalized with the SCTransform approach and then analyzed by conducting a mutual principal unit investigation (PCA) (https://satijalab.org/seurat/v3.1/integration.html) (13). The PCA analysis was also conducted on the integrated datasets with the cluster analysis being performed with uniform manifold approximation and projection (UMAP). Cluster analysis of single-cells was performed using Seurat’s graph-based clustering approach [R software package Seurat (version 2.3.4)] with the FindClusters feature resolution set to 0.1. Subsequently, the clusters were visualized by using the UMAP (version 0.2.6.0) graph. For quality control, unique molecular identifier counts of less than 500 and double multiples were removed. Furthermore, cells with > 5% mitochondrial genes and > 50% ribosomal genes were filtered out.

Collection and Processing of the HNSCC RNAseq Dataset

The RNA sequencing dataset for HNSCC and the corresponding clinically related data originated from the TCGA database (https://portal.gdc.carcinoma.gov/), which consists of 279 samples. The validation cohort dataset is TCGA provisional database, consisting of 249 samples. The raw gene expression dataset was processed. Probe IDs received the annotation toward the gene from the corresponding platform annotation profile of the GDC website and the raw matrix data received the quantile normalization and log2 conversion. Samples with missing data were excluded.

Building a CD8+ T Cell-Correlated Gene Signature

Single-cell data was classified into specific cell types and divided according to their respective tissue sources. The corresponding transcriptome investigation data were compared in order to screen for DEGs. To increase the efficiency of the study the candidate DEGs were taken as min.pct > 0.25 and | Log2 (FC) | > 0.5.

The association between HNSCC tumor-infiltration CD8+ T cell-correlated DEGs and overall survival time in TCGA HNSCC cases was studied. Univariate Cox regression analysis was carried out for identifying the genes associated with survival (p value < 0.05). Subsequently, the significance of candidate genes was selected using variable importance in a randomized survival forest (RSF) algorithm. A risk score model with the selected DEGs was built using multi-variate Cox regression approaches. In addition, the Kaplan-Meier test was employed for a number of gene features and p-values (log) were determined. Receiver operating characteristic (ROC) analysis was carried out for 3- and 5-year overall survival rates and area under the curves (AUCs) were determined for assessing the specificity and sensitivity of the gene signature. In addition, for testing the robustness of the results, the HNSCC tumor-infiltration CD8+ T cell-correlated gene signature was further verified with the TCGA HNSCC dataset.

The Effects of Age, Gender, Alcohol and Smoking on HNSCC Patients

To assess the correlation of risk score distribution and clinically related characteristics, HNSCC patients were grouped according to age and gender as well as their status regarding alcohol consumption and smoking. In addition, risk scores were calculated in order to assess the patients’ prognosis for HNSCC in the presence and absence of HPV by using multi-variate Cox regression correlations.

Statistical Analysis

Statistics investigations were carried out with R software (version 3.6.0). Kaplan-Meier tests and ROC analysis were performed with the “survivor” and “survROC” software packages (14). Optimal cutoff data points were calculated using the “survminer” package (15). Single-variate and multi-variate Cox regression correlations were used to assess the prognosis-correlated factors of interest. Hazard ratios (HR) and 95% confidence intervals (95% CI) were presented for all the prognosis-correlated factors. In statistical tests, P<0.05 was considered statistically significant.

Results

A Comparison of CD8+ T Cells in HPV− and HPV+ HNSCC Patients

From the TIMER2.0 website, a number of immunization deconvolution approaches including “XCELL (16)”, “MCPCOUNTER (17)”, “QUANTISEQ (18)”, “CIBERSORT-ABS” and “CIBERSORT (19)” were employed for estimating immunization infiltration of CD8+ T cells in HPV− and HPV+ HNSCC patients. With the single-variate Cox proportional risk model, the found that the tumor-infiltrated CD8+ T cells were protective for cases with HPV+ HNSCC, but this was not seen in HPV− HNSCC patients (Figure 1A). According to the Kaplan-Meier curves obtained, the survival period of the high tumor-infiltration CD8+ T-cell cohort was significantly longer than that of the low CD8+ T tumor infiltration seen in HNSCC patients who also had HPV, irrespective of the deconvolution approach used (Figure 1B).

FIGURE 1

Figure 1 Prognosis-correlated data of CD8+ T cells from patients with HPV− and HPV+ HNSCC. (A) A heat map of data from the multiple-variate Cox proportional risk model in terms of CD8+ T cells from patients with HPV− versus HPV+ HNSCC. The z-scores represent the risk scores. (B) Kaplan-Meier survival analysis showing the levels of CD8+ T cells from patients with HPV− (right column) and HPV+ (left column) HNSCC by using MCPcounter, CIBERSORT, CIBERSORT-ABS, XCELL and QUANTISEQ approaches.

Single-Cell RNA-Seq-Based DEGs Identification

The RNA-seq based on single cells consisted of 130,721 immunization cells from the tumor infiltrated leukocyte (TIL) samples obtained from HNSCC patients. Using the UMAP algorithm, this mixture of 11 cell types including CD8 T, CD4 T, dendritic, CD8+ T, mast, natural killer (NK) and plasma cells as well as mon/macrophages were unambiguously clustered and annotated (Figure 2A). Single-cell clustering was also based on the presence and absence of HPV (Figure 2B). According to the pie chart, the number of CD8+ T cells was an important unit of the HNSCC TILs (Figure 2C). The bar graphs also indicated that CD8+ T cells accounted for the greatest proportion of infiltration of all the immunization cells into the HPV− and HPV+ tumors in HNSCC patients (Figure 2D). The gene expression of LAG3 showed that exhausting T cells account for majority of the CD8 T cells (Figure 2E).

FIGURE 2

Figure 2 Identification of HPV− and HPV+ HNSCC tumor-infiltration of CD8+ T cell-correlated genes. (A) View of single cell samples from HPV− and HPV+ TIL. The annotated UMAP plot identifying 11 distinctive cell types. (B) The annotated UMAP plot of HPV− and HPV+ HNSCC TIL. (C) A pie chart of the seven cell types that make up the TIL of HNSCC. (D) Bar graphs of the cell proportions in the TIL of HPV− and HPV+ HNSCC. (E) Violin plots illustrating the expression of CD8A and exhausting T cell marker LAG3 in different TIL cell types.

Building a CD8+ T Cell-Correlated Gene Signature

Subsequently, the HPV+ HNSCC tumor-infiltration CD8+ T cell-correlated DEGs were screened based on the selection criteria in the approaches used (Figure 3A). To screen for the crucial survival-related factors, the DEGs from the CD8+ T cells of HPV+ HNSCC were analyzed using single-variate Cox regression for the TCGA dataset, and a total of 21 DEGs were identified to be significantly correlated to survival in these patients (p < 0.05) (Figure 3B). Based on the random forest algorithm, the top 5 significant genes, ACAP1, ANKRD28, C12orf75, M6PR and RGCC, were screened (Figure 3C).

FIGURE 3

Figure 3 A gene signature using infiltrated CD8+ T cells. (A) A volcano plot of the differentially expressed genes (DEGs) between HPV+ and HPV− HNSCC tumor-infiltrated CD8+ T cells. (B) A volcano plot showing the DEGs obtained from Cox regression analysis of survival-related HPV+ HNSCC-infiltrated CD8+ T cells. (C) Forest plot lines of the top 4 genes screened by using random survival forest analysis of HNSCC patients.

Further, we studied the expression of the significant genes in different immune cell subsets and found the association of ACAP1, ANKRD28, C12orf75 and M6PR, but not RGCC, with CD8+ T cells (Figure 4A). Also, we checked the TCGA data using CIBERSORT deconvolution and TIMER2 website, which demonstrated a substantial fraction of CD8+ T cells and the 4 signature genes in HPV+ compared to HPV- HNSCC (Figures 4B, C). RGCC was detectable, and therefore excluded from the signature genes.

FIGURE 4

Figure 4 Expression of the signature genes. (A) Violin plots illustrating the expression of the survival related signature genes in different TIL cell types of scRNA data. (B) Bulk seq data were deconvoluted using CIBERSORT method. Bar graph illustrating the proportion of infiltrated immune cells of HPV+ and HPV− HNSCC. (C) Expression of signature genes were analyzed using TIMER2 website. Box charts showed the expression of the signature genes of the TCGA data. “**” refers p < 0.01, “***” refers p < 0.001 and “****” refers p < 0.0001.

The risk scoring system was then built using these 4 genes with multi-variate Cox analysis using the TCGA dataset. In accordance with the formula, a risk score was calculated for the respective cases. The HPV+HNSCC cases in the TCGA dataset were then divided into high-risk and low-risk cohorts with the optimal cutoff data for the risk scores. Kaplan-Meier curves showed that the high-risk group survived for longer periods in comparison with those patients in the low-risk cohort (Figure 5A). In comparison, there was no clear distinction in the high-risk cohort of HPV- HNSCC case survival (Figure 5B). This was further validated in the TCGA provisional dataset for HPV+ HNSCC (Figure 5C). To estimate the predictive power of genetic characteristics, ROC curve analysis of the HPV+ HNSCC cases were plotted and this showed an AUC of 1 and 0.739 for 3 year-survival (Figure 5D), while ROC survival curves for the HPV- HNSCC patients were less significant (Figure 5E). This was verified by using the TCGA provisional dataset which showed an AUC curve of 0.839 for 3 year-survival (Figure 5F).

FIGURE 5

Figure 5 Validation of prognosis gene labels for HNSCC cases. (A) Kaplan-Meier (KM) analysis of the risk group that were defined with CD8+ T cell-correlated gene tags in the TCGA dataset for HPV+ HNSCC. (B) KM analysis of the risk group that were defined with CD8+ T cell-correlated gene tags in the TCGA dataset for HPV- HNSCC. (C) KM analysis of the risk group that were defined with CD8+ T cell-correlated gene tags in the TCGA provisional dataset for HPV+ HNSCC. (D) Three- and five-year ROC survival curves from the TCGA dataset for HPV+ HNSCC. (E) Three- and five-year ROC survival curves from the HPV- HNSCC TCGA dataset. (F) Three- and five-year ROC survival curves from the TCGA provisional dataset for HPV+ HNSCC.

Correlation of Risk Score Distribution and Clinical Characteristics of HNSCC Patients

HPV+ HNSCC cases in the TCGA dataset were categorized according to high or low risk score cohorts with the best cutoff data obtained. Box plots showed that age, sex and alcohol consumption were not correlated with the risk score (Figure 6A). However, smoking (Figure 6B) did show a correlation with the risk score. Furthermore, risk scores of HPV- HNSCC cases were also evaluated according to their gene signatures as established with the scRNA and TCGA datasets (Figures 6C, D). No correlation of the risk scores were found when the cases were divided according to their age and sex as well as alcohol consumption and smoking status. This revealed a specificity of the current gene signature for the assessment of smoking in HPV+ HNSCC patients.

FIGURE 6

Figure 6 Correlation between risk scores and clinically related characteristics of HNSCC patients. (A) Distribution of risk scores obtained when patients’ data obtained from the HPV+ HNSCC TCGA dataset were separated by age and sex. (B) Risk score distributions for alcohol and smoking in the HPV+ HNSCC TCGA dataset. (C) Distribution of risk scores obtained when patients’ data obtained from the HPV- HNSCC TCGA dataset were separated by age and sex. (D) Risk score distributions for alcohol and smoking in the HPV- HNSCC TCGA dataset. (E) Multi-variate Cox regression forest plots of the risk scores and clinically related characteristics in the HPV+ HNSCC TCGA dataset. (F) Multi-variate Cox regression forest plots of risk scores and clinically related characteristics in the HPV- HNSCC TCGA dataset. Alco and Smok refer to alcohol consumer and smoker.

In order to compare the prognosis-correlated factors to those general factors, risk scores for genetic characteristics and clinically-related variables were analyzed by multi-variate Cox regression (Figures 6E, F). The forest plots did not show any significance that the current risks examined, thus revealing the significance of the current gene signature and risk scoring system used in this study.

Profiling the Gene Expression of the HPV+ HNSCC Risk Groups

With the signature genes, we separated the HPV+ HNSCC into high and low risk groups. We investigated the DEGs between high and low risk groups of HPV+ HNSCC TCGA samples (Figure 7A). The gene ontology (GO) and KEGG pathways studies showed that the risk DEGs are enriched in the T cell activation and differentiation (Figures 7B, C), indicating the involvement of the signature genes in the T cells function. We further investigated the correlation of the signature genes with the immune checkpoint genes and found a decent correlation with CTLA4, LAG3 and PDCD1, but not with the DAMP signal gene S100A8 (Figures 8A–D). Especially, a strong correlation was detected between ACAP1 and the checkpoint related genes (Figure 8E).

FIGURE 7

Figure 7 Profiling the gene expression of the HPV+ HNSCC risk groups. (A) A volcano plot of the DEGs between high and low risk groups of HPV+ HNSCC TCGA samples. (B) Bar graphs showing the enriched gene ontology (GO) Biological Process of the risk DEGs and (C) the enriched KEGG pathways.

FIGURE 8

Figure 8 Correlation of the signature genes with the checkpoint related genes. Correlation of the signature genes with (A) CTLA4 (B)LAG3 and (C) PDCD1. (D) S100A8 was used as control. (E) A corr plot showing the signature genes with the checkpoint related genes.

Discussion

Using a single-variate Cox proportional risk model, we identified that the tumor-infiltration of CD8+ T cells were protective for cases with HPV+ HNSCC, but not for HPV- HNSCC cases. Further, the tumor immunization environment was explored using single-cell sequencing and screened for CD8+ T cell-specific gene feature differences between HPV+ and HPV- HNSCC patients. In addition, we built a prognosis-correlated genetic signature that divided the overall survival of HPV+ HNSCC into two risk groups with the high-risk cases showing poorer prognoses. A prognosis-correlated gene signature consisting of 4 genes with low risk. Then the association between CD8+ T cell genetic traits and clinically related parameters were determined using TCGA dataset and verified in TCGA provisional dataset, to demonstrate the accuracy of genetic traits for prognosis-correlated prediction.

According to the single-cell data, there were significant differences between the immune profiles of HPV- versus HPV+ HNSCC patients, which is of significance when designing their immunotherapy regimens. Although the proportion infiltrating CD8 cells were similar in the two HNSCC types, the CD8 gene cell expression profiles were not identical, so more tailored therapies will be required in order to improve the survival rates of patients. The distinct immune profiles in the microenvironment of HPV+ and HPV- HNSCC patients may result from the presence of viral antigens throughout the carcinogenesis process, resulting in the early innate immune responses and the enhancement of the T cells adaptive immune response (20, 21). Further comparisons of the gene transcriptomes of cells with and without the presence of viral antigens throughout the carcinogenesis process is needed to enable the true cause of immune profile differences seen in the presence and absence of HPV in HNSCC patients.

In line with our current results, a recent study noted that between nearly 10% of the infiltrating T cells present in HNSCC were such CD8-positive T cells that target HPV and express PD-1 (22). One of these cells has stem cell characteristics and can expand in large numbers to treat head and neck cancer if cancer immunotherapy is used. Current conventional treatments for head and neck cancer include radiotherapy and chemotherapy, but they may affect the number of immune cells in the body. Therefore, better results may be achieved if immunotherapy is used first and then combined with conventional modalities.

The shortcoming of our current study is the small size of the test group, which, although showing a good trend, somehow lacks statistical significance. Further studies are needed to examine not only HPV+ samples but also specific epitopes of HPV infection.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics Statement

The studies involving human participants were reviewed and approved by Youjiang Medical University for Nationalities. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

JS and HL designed this study. YW, LM, KC, and JZ analyzed the scRNA-seq data. HL, SH, and JS analyzed the bulk sequencing data. QW, ZW, and SS contributed to scRNA-seq data and helped organize this manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research work was funded by the National Natural Science Foundation of China (#31970745).

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

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

References

1. Argiris A, Karamouzis MV, Raben D, Ferris RL. Head and Neck Cancer. Lancet (2008) 371:1695–709. doi: 10.1016/S0140-6736(08)60728-X

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Geng X, Zhang Y, Zeng Z, Zhu Z, Wang H, Yu W, et al. Molecular Characteristics, Prognostic Value, and Immune Characteristics of M(6)A Regulators Identified in Head and Neck Squamous Cell Carcinoma. Front Oncol (2021) 11:629718. doi: 10.3389/fonc.2021.629718

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Timar J, Csuka O, Remenar E, Repassy G, Kasler M. Progression of Head and Neck Squamous Cell Cancer. Cancer Metastasis Rev (2005) 24:107–27. doi: 10.1007/s10555-005-5051-5

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, Activity, and Immune Correlates of Anti-PD-1 Antibody in Cancer. N Engl J Med (2012) 366:2443–54. doi: 10.1056/NEJMoa1200690

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Burtness B, Harrington KJ, Greil R, Soulieres D, Tahara M, de Castro G Jr., et al. Pembrolizumab Alone or With Chemotherapy Versus Cetuximab With Chemotherapy for Recurrent or Metastatic Squamous Cell Carcinoma of the Head and Neck (KEYNOTE-048): A Randomised, Open-Label, Phase 3 Study. Lancet (2019) 394:1915–28. doi: 10.1016/S0140-6736(19)32591-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Ferris RL, Licitra L, Fayette J, Even C, Blumenschein G Jr, Harrington KJ, et al. Nivolumab in Patients With Recurrent or Metastatic Squamous Cell Carcinoma of the Head and Neck: Efficacy and Safety in CheckMate 141 by Prior Cetuximab Use. Clin Cancer Res (2019) 25:5221–30. doi: 10.1158/1078-0432.CCR-18-3944

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Ang KK, Harris J, Wheeler R, Weber R, Rosenthal DI, Nguyen-Tan PF, et al. Human Papillomavirus and Survival of Patients With Oropharyngeal Cancer. N Engl J Med (2010) 363:24–35. doi: 10.1056/NEJMoa0912217

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Carvalho AL, Nishimoto IN, Califano JA, Kowalski LP. Trends in Incidence and Prognosis for Head and Neck Cancer in the United States: A Site-Specific Analysis of the SEER Database. Int J Cancer (2005) 114:806–16. doi: 10.1002/ijc.20740

PubMed Abstract | CrossRef Full Text | Google Scholar

9. LeHew CW, Weatherspoon DJ, Peterson CE, Goben A, Reitmajer K, Sroussi H, et al. The Health System and Policy Implications of Changing Epidemiology for Oral Cavity and Oropharyngeal Cancers in the United States From 1995 to 2016. Epidemiol Rev (2017) 39:132–47. doi: 10.1093/epirev/mxw001

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Cillo AR, Kurten CHL, Tabib T, Qi Z, Onkar S, Wang T, et al. Immune Landscape of Viral- and Carcinogen-Driven Head and Neck Cancer. Immunity (2020) 52:183–99.e9. doi: 10.1016/j.immuni.2019.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Kansy BA, Concha-Benavente F, Srivastava RM, Jie HB, Shayan G, Lei Y, et al. PD-1 Status in CD8(+) T Cells Associates With Survival and Anti-PD-1 Therapeutic Outcomes in Head and Neck Cancer. Cancer Res (2017) 77:6353–64. doi: 10.1158/0008-5472.CAN-16-3167

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for Analysis of Tumor-Infiltrating Immune Cells. Nucleic Acids Res (2020) 48:W509–14. doi: 10.1093/nar/gkaa407

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, et al. Comprehensive Integration of Single-Cell Data. Cell (2019) 177:1888–1902 e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Therneau TM, Li H. Computing the Cox Model for Case Cohort Designs. Lifetime Data Anal (1999) 5:99–112. doi: 10.1023/A:1009691327335

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zeng L, Fan X, Wang X, Deng H, Zhang K, Zhang X, et al. Bioinformatics Analysis Based on Multiple Databases Identifies Hub Genes Associated With Hepatocellular Carcinoma. Curr Genomics (2019) 20:349–61. doi: 10.2174/1389202920666191011092410

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Aran D, Hu Z, Butte AJ. Xcell: Digitally Portraying the Tissue Cellular Heterogeneity Landscape. Genome Biol (2017) 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression. Genome Biol (2016) 17:218. doi: 10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Finotello F, Mayer C, Plattner C, Laschober G, Rieder D, Hackl H, et al. Molecular and Pharmacological Modulators of the Tumor Immune Contexture Revealed by Deconvolution of RNA-Seq Data. Genome Med (2019) 11:34. doi: 10.1186/s13073-019-0638-6

PubMed Abstract | CrossRef Full Text | Google Scholar

19. 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

20. Menezes FDS, Fernandes GA, Antunes JLF, Villa LL, Toporcov TN. Global Incidence Trends in Head and Neck Cancer for HPV-Related and -Unrelated Subsites: A Systematic Review of Population-Based Studies. Oral Oncol (2021) 115:105177. doi: 10.1016/j.oraloncology.2020.105177

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Kobayashi K, Hisamatsu K, Suzui N, Hara A, Tomita H, Miyazaki T. A Review of HPV-Related Head and Neck Cancer. J Clin Med (2018) 7(9):241. doi: 10.3390/jcm7090241

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Eberhardt CS, Kissick HT, Patel MR, Cardenas MA, Prokhnevska N, Obeng RC, et al. Functional HPV-Specific PD-1(+) Stem-Like CD8 T Cells in Head and Neck Cancer. Nature (2021) 597(7875):279–84. doi: 10.1038/s41586-021-03862-z

PubMed Abstract | CrossRef Full Text | Google Scholar