ORIGINAL RESEARCH article
Discovery and Verification of an Immune-Related Gene Pairs Signature for Predicting Prognosis in Head and Neck Squamous Cell Carcinoma
- 1Department of Hand and Microsurgery, Xiangya Hospital of Central South University, Changsha, China
- 2Department of Neurosurgery, Huashan Hospital, Fudan University, Shanghai, China
- 3Department of Anesthesiology, The First Affiliated Hospital of USTC, University of Science and Technology of China, Hefei, China
The study of IRGPs to construct the prognostic signature in head and neck squamous cell carcinoma (HNSCC) has not yet elucidated. The objective of this study was to explore a novel model to predict the prognosis of HNSCC patients. The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets were set as training and validation cohorts, respectively. The least absolute shrinkage and selection operator (LASSO) and time-dependent ROC were employed to screen the highest frequency immune-related gene pairs (IRGPs) and their best cut-off value. Survival analysis, Cox regression analysis were applied to discover the effects of selected IRGPs signature on survival outcomes. The immune cell proportions were deconvoluted by the CIBERSORT method. After a couple of filtering, we obtained 22 highest frequency IRGPs. The overall survival time of HNSCC patients with a high score of IRGPs was shorter as compared to the ones with a low score in two independent datasets (P < 0.001). Six kinds of immune cells were found to be differentially distributed in the two different risk groups of HNSCC patients (P < 0.001). GO and GSEA analysis showed these differentially expressed genes enriched in multiple molecular functions. The new IRGPs signature probably confers a new insight into the prognosis prediction of HNSCC patients.
Head and neck squamous cell carcinoma (HNSCC) is one of the most common cancers globally, and more than 600,000 cases are diagnosed annually (Leemans et al., 2011). Smoking, heavy alcohol, and HPV infection are the main potential causes of HNSCC development. Local or distant metastasis is the major factor for cancer-associated death (Huang et al., 2018). Despite multimodal treatments, such as surgery, chemoradiotherapy, and targeted therapies, the 5-year survival rate of HNSCC patients remains poor (Li et al., 2020; Zhang et al., 2020). Interestingly, even patients with the same clinical symptoms, TNM stage, and treatments would have different prognoses (Karam and Raben, 2019). Therefore, the inaccurate effects of the TNM stage on predicting the prognosis of HNSCC patients lead to a clinically meaningful puzzle to explore and develop some novel and better molecular models.
Immunity is a disease defense system composed of a series of biological structures and processes in an organism (Montecino-Rodriguez et al., 2013). It is composed of various proteins, cells, organs, and tissues. Disorders of the immune system can lead to a variety of diseases (O’Byrne and Dalgleish, 2001; de Visser et al., 2006). Many studies have confirmed that the immune system plays an important role in the occurrence and development of tumors. Charlotte Lecerf et al. (2019) recently found OX40L, PDGFRB, and PD-1 act as the potential prognostic biomarkers and major therapeutic targets. Moreover, blockade to the PD1/PD-L1 axis has been approved by the FDA and shows groundbreaking progress in improving HNSCC patients’ survival time (Tsangaris et al., 1969). Discovering and understanding the relationship between immune-related genes and HNSCC could be crucial.
The tumor microenvironment (TME) includes various cell types and extracellular components, such as vasculature, cancer-associated fibroblasts, immune cells, tumor-associated endothelial cells, and extracellular matrix (Wu and Dai, 2017). Accumulated studies show that TME influences therapeutic response and clinical outcomes. Regulation of TME by CTX/L-NIL could improve the sensitivity of tumor cells of chemoradiotherapy resistance (Hanoteau et al., 2019). The expression level of PD-L1 in exosomes of HNSCC correlated with the clinical outcome of patients. The blocking of PD-L1+exosomes weakens immune suppression in HNSCC (Theodoraki et al., 2018).
Previous studies showed many gene-expression-based models for predicting cancer patient prognosis (Resteghini et al., 2018). Unfortunately, none of them could be implemented in routine clinical practice because of some issues like overfitting on small discovery data sets and lack of sufficient validation. With more and more large-scale gene expression profiles published online, it could identify potentially more reliable prognosis-related biomarkers in HNSCC. However, these diverse large-scale data also represent a daunting challenge. For instance, the traditional method to normalize gene expression is a difficult task given the potential biological heterogeneity among data sets and technical biases across measurement platforms. Instead, the relative ranking of gene expression may be a good way to eliminate the requirement for data scaling and normalization (Li et al., 2017), producing a robust result in various applications.
Considering that the important roles of immune-related genes, there are emerging studies that calculate the score of immune-related gene pairs (IRGPs) based on its relative raking and report a significant correlation between IRGPs and the prognosis of cancer, including hepatocellular carcinoma, colorectal cancer, gastric cancer, lung adenocarcinoma and cervical cancer (Wu et al., 2019; Zhang et al., 2019; Nie et al., 2020; Sun et al., 2020; Zhao et al., 2021). To our best knowledge, no study applying IRGPs to construct the prognostic signature in HNSCC. Here, we used TCGA datasets, GEO datasets, and immune-related genes from the ImmPort database to establish and validate the immune prognostic signature, which was constructed by 22 immune-related gene pairs consisting of 28 unique immune-related genes. The relationships of IRGPs with clinical features and underlying molecular function enrichments and immune cell proportions were further discovered in our study.
Materials and Methods
Acquisition of Relevant Public Data
The RNA-seq expression data and clinical data of 525 HNSCC patients in the training group were downloaded from the TCGA dataset1, which was used as the training cohort. Two independent datasets were used as the validation group were downloaded from the GEO dataset in this study, including GSE65858 (N = 270, Illumina HumanHT-12 V4.0 expression bead chip) (Wichmann et al., 2015) and GSE2379 (N = 34, Affymetrix Human Genome U95 Version 2 Array and Affymetrix Human Genome U95A Array) (Cromer et al., 2004). For each dataset, the probe ID was converted to the corresponding gene symbol according to its annotation file without making further standardization. If more than one probe matches the same gene, the average value is calculated as the expression value of the gene. The immune-related gene list was download from the ImmPort database2.
Establishment and Validation of an Immune Gene-Related Prognostic Model
All the 2,498 immune-related genes from the ImmPort database were selected for the candidates’ immune-related genes to be intersected with the TCGA dataset and GEO dataset, respectively. The immune-related genes obtained from the above three cohorts are intersected to obtain a matrix that included 145 co-IRGs. During screen immune-related genes (IRG), the median absolute deviation must be greater than 0.5. Next, each gene will have a chance to be paired up with another gene in the matrix. The estimated score of the gene pair is up to the reduction from the gene counts. If the gene pair i is made up of genei and genei+1. The score will be 1 if the expression count of genei is higher than genei+1; otherwise, its score will be 0. More importantly, we should discard the gene if the percentage of the score was 1 or 0 was less than 20%. We used LASSO to get the maximum information but the parsimonious model (iteration = 1,000). The IRGPs’ signature was calculated as follows: Risk score = (Exprgenepair–1 × Coefgenepair–1) + (Exprgenepair–2 × Coefgenepair–2) + … + (Exprgenepair–n × Coefgenepair–n). Expr means the expression value of gene pairs, and Coef means the regression coefficient derived. Then we used the TCGA survival data, risk score, and time-dependent ROC at 3 years in the TCGA dataset to screen the highest frequency IRGPs and calculate the cut-off value to determine the patient with high risk or low score. The capability of this model to predict overall survival (OS) was visualized and compared by Kaplan-Meier, univariate and multivariate Cox proportional hazard analysis. The external dataset GSE65858 and GSE22379 were used for the validations.
Correlations Between IRGPs and Tumor-Infiltrating Immune Cells Proportions
The CIBERSORT algorithm, a method for characterizing cell composition from gene expression profiles of complex tissues (Newman et al., 2015), was used to explore the relationship between tumor immune cells and IRGPs. It should enable large-scale analysis for inferring immune cells, cellular biomarkers in tumor transcriptomes (Gentles et al., 2015). The proportions of 22 immune cell types were estimated in the TCGA cohort. A radar map is drawn for comparing these immune cell distributions between the patients with a high score and the ones with a low score according to the results of the CIBERSORT algorithm.
Gene Enrichment Analysis
Gene ontology (GO) analysis was used to discover the enriched signaling pathway based on the genes in the IRGPs prognostic signature in HNSCC. Gene set enrichment analysis (GSEA) was used to understand the potential biological mechanisms these genes in IRGPs prognostic signature may involve.
All statistical analyses were performed using R (version 4.0.23). The log-rank test was used in the Kaplan-Meier analysis. The ROC curve was performed by the “survivalROC” package. The “survival” package was used to perform survival analysis curves, univariate Cox survival analysis and multivariate Cox survival analysis. The radar map was drawn by the “fmsb” package. The GSEA curve and the bubble chart of GO enrichment analysis were performed by the “fgsea” and “ggplot2” packages. For all tests, the significant differences were considered at P < 0.05.
Construction and Validation of IRGPs Signature in HNSCC
The flow chart of this study is shown in Figure 1. The characterizations of the training and validation cohort are summarized in Table 1. As showing in Table 2, 22 highest frequency IRGPs were obtained from 2,498 immune-related genes in the ImmPort database after a couple of filtering mentioned above. The best cut-off value for the determination of high or low scores was calculated by ROC analysis in the TCGA dataset of HNSCC patients (sensitivity = 0.701, specificity = 0.812) (Figure 2A). The Kaplan Merrier analysis showed that patients with a high score of IRGPs had been predicted a poorer outcome as opposed to a low score of IRGPs (P < 0.001, Figure 2B), which was further validated in the GSE65858 cohort (P < 0.05, Figure 2C) and GSE2379 (P = 0.032, Figure 2D).
Figure 2. Construction and verification of IRGPs signature. (A) ROC curve analysis was used to get the best cut-off value of the score to divided patients into high score group or low score group. Kaplan-Meier survival curves assessed the effect between high score group and low score group in the TCGA dataset (B), validated in GSE65858 (C) dataset and GSE22379 dataset (D) (P < 0.05).
To test whether the effect of IRGPs on HNSCC prognosis is independent of clinical parameters like age, clinical stage, lymph node status, etc., univariable and multiple variable Cox proportional regression analysis were performed. The death hazard ratio of univariate regression analyses of a high score is independently and significantly higher over a low score in the training TCGA cohort (HR = 3.967, 95% CI 3.026–5.200, P < 0.001, Figure 3A), which was validated both in the GSE65858 cohort (HR = 1.404, 95% CI 1.009–1.952, P = 0.021, Figure 3C) and GSE2379 cohort (HR = 9.572, 95% CI 1.430–64.051, P = 0.02, Figure 3E). The same results in multivariate regression analysis. TCGA cohort (HR = 3.844, 95% CI 2.911–5.076, P < 0.001, Figure 3B), GSE65858 (HR = 1.609, 95% CI 1.058–2.361, P = 0.0115, Figure 3D) and GSE2379 (HR = 23.838, 95% CI 1.607–67.926, P = 0.023, Figure 3F).
Figure 3. Explored and validation of independent clinical prognostic factors related to the signature of IRGPs in HNSCC. Univariate analysis and multivariate analysis were used to explored independent clinical prognostic factors in the TCGA dataset (A,B), validated in GSE65858 dataset (C,D) and GSE2379 dataset (E,F).
The Differential Tumor-Infiltrating Immune Cells Between HNSCC Patients With Different IRGPs
There are many different proportions of immune cells found in HNSCC patients with high-score over low-score of IRGPs (Figure 4). As shown in the radar map (Figure 4A), we found that naïve B cells, Plasma cells, Resting mast cells, CD8+ T cells, and Treg cells were significantly enriched in patients with a low score of IRGPs than the high score of IRGPs (P < 0.001) (Figures 4B–F). However, M0 macrophages were significantly higher in the high score group than the low score group (Figure 4G). The comparisons about the other immune cells in the two groups of patients were shown in Supplementary Figure 1. Furthermore, we explored the potential associations between the IRGPs’ signature and 22 immune cells by Pearson correlation analysis. We found that the signature risk score was significantly and negatively correlated to resting mast cells (R = –0.28), CD8+ T cells (R = –0.36), naïve B cells (R = –0.31), Treg cells (R = –0.19), and plasma cells (R = –0.27) (Figures 5A–E, all P < 0.01). However, the signature risk score was positively correlated to M0 macrophages (Figure 5F, R = 0.35, P < 0.01).
Figure 4. Immune cells correlation of IRGPs signature in HNSCC. (A) Using CIBERSORT calculated the abundance of 22 different immune cells in high score group and low score group of HNSCC. The abundance distribution of (B) Resting mast cells, (C) CD8+ T cells, (D) Naïve B cells, (E) Treg cells, (F) Plasma cells, and (G) M2 macrophages within the two groups (P < 0.05).
Figure 5. Relationship determined by Pearson correlation analysis between IRGPs and infiltration relative abundances of immune cells. (A) Resting mast cells, (B) CD8+ T cells, (C) Naïve B cells, (D) Treg cells, (E) T follicular helper cells, (F) Plasma cells.
Performance of IRGPs on Predicting HNSCC Prognosis
We applied the ROC analysis to explore the predictive ability of the IRGPs’ signature on the 5-year overall survival in the two validated HNSCC cohorts. The value of under the curve (AUC) in GSE65858 and GSE2379 was 0.765 and 0.905, respectively (Figures 6A,B).
Figure 6. The potential signaling pathways of the signature of IRGPs were explored by GEO enrichment analysis. (A,B) The receiver operating characteristic (ROC) curve illustrated the predictive ability of the 22 IRGPs signature. The areas under the curves for 5 years’ survival were 0.765 and 0.905 in two validation cohorts, respectively. (C) The signature of IRGPs in the TCGA cohort of HNSCC was involved in the positive regulation of cell activation, regulation of lymphocyte activation, regulation of immune effector process, humoral immune response, plasma membrane, and phagocytosis.
Enrichment Analysis Explored Potential Signaling Pathways
Go analysis showed that the signature of IRGPs in the TCGA cohort of HNSCC was involved in the positive regulation of cell activation, regulation of lymphocyte activation, regulation of immune effector process, humoral immune response, plasma membrane and phagocytosis (Figure 6C). According to the GSEA analysis, we revealed positive regulation of cell activation, regulation of lymphocyte activation, regulation of immune effector process, regulation of gene expression epigenetic, phagocytosis, mRNA binding, membrane signaling receptor complex, production of molecular mediator, and humoral immune response (Figure 7). Based on their normalized enrichment score, a lot of enriched signaling pathways were identified. It mainly including regulation of cell activation, regulation of lymphocyte activation, regulation of immune effector process, phagocytosis, and humoral immune response, which may be the immune-related pathways of significant enrichments in HNSCC. It provides possible immunologic concepts to support the prediction of IRGPs for HNSCC prognosis.
Figure 7. GSEA enrichment plots from gene set enrichment analysis. The results show that nine cancer hallmark gene set included (A) positive regulation of cell activation, (B) regulation of lymphocyte activation, (C) Regulation of immune effector process, (D) Regulation of gene expression epigenetic, (E) Phagocytosis, (F) mRNA binding, (G) Membrane signaling receptor complex, (H) Production of molecular mediator, (I) Humoral immune response are significant enrichment of low score of immune-related pathways in patients with HNSCC (P < 0.05).
HNSCC is the sixth most prevalent malignant tumor worldwide. The immune system is the main defense system of people against diseases, including cancers (Huse, 2017). Exploring the mechanisms between the immune-related genes and cancer prognosis may open a new perspective for predicting the outcomes of HNSCC patients. Currently, immunotherapy, receiving extensive attention all around the world, has been demonstrated as a groundbreaking remedy for cancer patients by blocking CTLA-4 or PD-1 pathways (Yang, 2015). Besides, researchers show that human peripheral blood T lymphocytes transduced by CD19-CAR could treat leukemia well (Brentjens et al., 2003). Therapeutic cancer vaccines and tumor neoantigens are other options for joining the armies of immunotherapy. Due to a lack of ideal vaccine design and the existence of an immunosuppressive tumor microenvironment, therapeutic vaccination remains suboptimal (Melief et al., 2015).
Continued efforts to understanding the role of the immune microenvironment will provide better methods for the treatment of HNSCC. So, we applied 22 immunes relate gene pairs to construct a prognostic signature and further explore its potential role in HNSCC prognosis. A previous study has reported that immune-related genes provide a model for predicting outcomes in HNSCC (She et al., 2020). However, a limitation is the lack of robust normalization to reduce batch effects if different platforms or sources are involved. Therefore, the first step needs to standardize all genes to reduce the errors caused by different samples and platforms. To overcome this limitation, we adopted a method named IRGPs (validated in a variety of cancers; Peng et al., 2018) to assess the score of the sample by using the comparable relative score of a pair of gene expressions in the same sample, which could remove the biases from the different platforms and samples.
In this study, 515 HNSCC patients from TCGA datasets act as the training group, the data of 270 patients from GSE65858 and 34 patients from GSE2379 act as the validation groups. Lasso penalized Cox regression was used to construct a 22-IRGPs signature and ROC curve analysis was used to calculate the best cut-off value. The HNSCC patients with a high or low score of IRGPs have a significant difference in survival analysis. Univariate and multivariate Cox regression analyses showed the high score was an independent prognostic factor. Besides, we used two external cohorts to confirm the conclusions, which remain the same as in the training group. Consequently, our signature has high accuracy and could be a potential prognostic factor in HNSCC. We found that naïve B cells, Plasma cells, Resting mast cells, CD8+ T cell, and Treg cells were significantly enriched in patients with a low score of IRGPs by immunocyte infiltration analysis, which may play an inhibitory role in the development of HNSCC. It has report showed that the proportions of naïve B cells, resting mast cells, and CD8+ T cells in non-cancer tissues were significantly higher than HNSC tissues (Krishna et al., 2018). It is further confirmed the accuracy of the results of our study. Furthermore, some research also found that plasma cells were associated with distant cancer metastasis, clinical stage. Naïve B cells and Tregs were significantly associated with pathological grade (Liang et al., 2020). These differential proportions of tumor-infiltrating cells may be reasons for the different overall survival observed and validated between a high and low score of IRPGs and associated with the IRGP interactions in HNSCC patients.
The 22 IRGPs in our model consisted of 28 immune-related genes, including Wnt5a, CCL1, CCL5, CCL9, HLA-DOB, IFN, and CD247, which mainly involved some important pathways and cellular processes such as proinflammatory, antimicrobial responses, antigen identification, and immune cell signaling. Wnt signaling regulates a series of signaling pathways, including cell proliferation, differentiation, adhesion, and motility (Wend et al., 2010). Wnt5a act as a non-canonical Wnt ligand that plays a critical role during growth and development (Oishi et al., 2003). The abnormal activation and inhibition of Wnt5a lead to the genesis in tumors. Increased Wnt5a level promotes the ability of invasion in metastasis melanoma (Weeraratna et al., 2002). In contrast, Wnt5a acts as a potential inhibitor. The increased Wnt5a level inhibits cell growth and migration in thyroid carcinoma (Kremenevskaja et al., 2005). Furthermore, several reports have already suggested that Wnt5a plays an important role in cancer-associated inflammation. Wnt5a induced gastric epithelial cells to produce IL-1, leading macrophages gathered in gastric mucosa inducing gastric inflammation (Asem et al., 2016). Most of the research reported that chemokines and their receptors play an important role in mediated cancer progress and metastasis. For instance, CCL1 derived from tumor-associated macrophages promotes breast cancer metastasis through activating NF-κB/SOX4 signaling (Wang et al., 2018). CXCL9 predominantly mediates lymphocytic infiltration to the focal sites and suppresses tumor growth (Tokunaga et al., 2018). The knockout of NLRX1 could significantly decrease IFN-I dependent T cell infiltration and inhibit tumor progression (Luo et al., 2020). HLA-DOB consists of a class II molecule and a beta chain and expresses on the antigen-presenting cells and B lymphocytes (Pruitt et al., 2009). The above descriptions indicated that immune response contributes to the development and progression of carcinomas in the interaction of the immune gene. It is consistent with our results that IRGPs were significant associated with some kinds of immune cells, including CD8+ T cells, Treg cells, and Resting mast cells (Li et al., 2021). The immune-related genes inside our IRGPs might play an important role in the tumor microenvironment of HNSCC. The GSEA and GO enrichment analysis implied that regulation of positive regulation of cell activation, regulation of lymphocyte activation, regulation of immune effector process, regulation of gene expression epigenetic, phagocytosis, mRNA binding, membrane signaling receptor complex, production of molecular mediator and humoral immune response might be the significant enrichment of related pathways in HNSCC, which give us more confidence in the capability of IRGPs to predict the HNSCC prognosis as we all know that the expression of gene reduced can promote most of diseases development. For example, mice with the knockout of EIF4E-BP1, a direct target of mTOR, formed more and larger lesions than its wild type mice (Wang et al., 2019). The humoral immune response is an important part of the immune system. Since tumors originate from containing self-antigens of autologous cells, the abnormal exposure of these antigens can activate the autoimmune response to promote the development of cancer (Zaenker et al., 2016). Therefore, we suggested that the 22 IRGPs signature was constructed in our research act as an important role in the development and prognosis of HNSCC. The reason that we chose two datasets GSE65858 and GSE22379 as validation cohorts is that these two databases of HNSCC contain the most complete clinical data. Although the second database contains a relatively small number of patients, it contains comprehensive patients’ information such as survival time, survival status, and TMN staging that we want. The clinical data of age in the multivariate COX regression analysis among three groups are not in the same range, P < 0.001, =0.036, =0.479, respectively, in Figures 3B,D,F. Therefore, we don’t think age can be used as a prognostic factor for HNSCC.
However, we are obliged to admit the limitations of our research. Firstly, although we used the validation group to prove the training group, the number of samples is not large enough. Secondly, it can be more persuasive if we further validate our findings through RT-QPC or IHC.
In conclusion, we established a new IRGPs signature for predicting the prognosis in HNSCC, which may provide a possible tool to help clinical physicians to give HNSCC patients consult about the prognosis.
We developed a new IRGPs signature prognostic model in HNSCC. It probably confers a new insight into the prognosis prediction of HNSCC patients, which may provide some clues for future personalized medicine decisions.
Data Availability Statement
The raw data of this study are derived from the TCGA database (https://portal.gdc.cancer.gov/) and GEO data portal (https://www.ncbi.nlm.nih.gov/geo/; Accession Nos.: GSE65858 and GSE2379), which are publicly available databases. The R codes will be available upon request.
MH: conception, revision for important intellectual content, and supervision. JH: interpretation or analysis of data. XF: preparation of the manuscript. All authors contributed to the article and approved the submitted version.
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.
We are greatly thankful to Prof. Gangcai Zhu and Zhexuan Li at the Xiangya Hospital of Central South University for friendly guidance.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.654657/full#supplementary-material
Supplementary Figure 1 | Any other immune cells enriched of IGRPs signature in HNSCC.
Brentjens, R. J., Latouche, J.-B., Santos, E., Marti, F., Gong, M. C., Lyddane, C., et al. (2003). Eradication of systemic B-cell tumors by genetically targeted human T lymphocytes co-stimulated by CD80 and interleukin-15. Nat. Med. 9, 279–286.
Cromer, A., Carles, A., Millon, R., Ganguli, G., Chalmel, F., Lemaire, F., et al. (2004). Identification of genes associated with tumorigenesis and metastatic potential of hypopharyngeal cancer by microarray analysis. Oncogene 23, 2484–2498.
Gentles, A. J., Newman, A. M., Liu, C. L., Bratman, S. V., Feng, W., Kim, D., et al. (2015). The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat. Med. 21, 938–945. doi: 10.1038/nm.3909
Hanoteau, A., Newton, J. M., Krupar, R., Huang, C., Liu, H.-C., Gaspero, A., et al. (2019). Tumor microenvironment modulation enhances immunologic benefit of chemoradiotherapy. J. Immunother. Cancer 7:10. doi: 10.1186/s40425-018-0485-9
Huang, D., Qiu, Y., Li, G., Liu, C., She, L., Zhang, D., et al. (2018). KDM5B overexpression predicts a poor prognosis in patients with squamous cell carcinoma of the head and neck. J. Cancer 9, 198–204. doi: 10.7150/jca.22145
Krishna, S., Ulrich, P., Wilson, E., Parikh, F., Narang, P., Yang, S., et al. (2018). Human papilloma virus specific immunogenicity and dysfunction of CD8 T cells in head and neck cancer. Cancer Res. 78, 6159–6170. doi: 10.1158/0008-5472.CAN-18-0163
Lecerf, C., Kamal, M., Vacher, S., Chemlali, W., Schnitzler, A., Morel, C., et al. (2019). Immune gene expression in head and neck squamous cell carcinoma patients. Eur. J. Cancer 121, 210–223. doi: 10.1016/j.ejca.2019.08.028
Li, B., Cui, Y., Diehn, M., and Li, R. (2017). Development and validation of an individualized immune prognostic signature in early-stage nonsquamous non-small cell lung cancer. JAMA Oncol. 3, 1529–1537. doi: 10.1001/jamaoncol.2017.1609
Li, Z., Chen, C., Wang, J., Wei, M., Liu, G., Qin, Y., et al. (2021). Overexpressed PLAU and its potential prognostic value in head and neck squamous cell carcinoma. PeerJ. 9:e10746. doi: 10.7717/peerj.10746
Li, Z., Chen, X., Wei, M., Liu, G., Tian, Y., Zhang, X., et al. (2020). Systemic analysis of RNA alternative splicing signals related to the prognosis for head and neck squamous cell carcinoma. Front. Oncol. 10:87. doi: 10.3389/fonc.2020.00087
Luo, X., Donnelly, C. R., Gong, W., Heath, B. R., Hao, Y., Donnelly, L. A., et al. (2020). HPV16 drives cancer immune escape via NLRX1-mediated degradation of STING. J. Clin. Invest. 130, 1635–1652. doi: 10.1172/JCI129497
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi: 10.1038/nmeth.3337
Oishi, I., Suzuki, H., Onishi, N., Takada, R., Kani, S., Ohkawara, B., et al. (2003). The receptor tyrosine kinase Ror2 is involved in non-canonical Wnt5a/JNK signalling pathway. Genes Cells 8, 645–654.
Peng, P.-L., Zhou, X.-Y., Yi, G.-D., Chen, P.-F., Wang, F., and Dong, W.-G. (2018). Identification of a novel gene pairs signature in the prognosis of gastric cancer. Cancer Med. 7, 344–350. doi: 10.1002/cam4.1303
She, Y., Kong, X., Ge, Y., Yin, P., Liu, Z., Chen, J., et al. (2020). Immune-related gene signature for predicting the prognosis of head and neck squamous cell carcinoma. Cancer Cell Int. 20:22. doi: 10.1186/s12935-020-1104-7
Sun, X.-Y., Yu, S.-Z., Zhang, H.-P., Li, J., Guo, W.-Z., and Zhang, S.-J. (2020). A signature of 33 immune-related gene pairs predicts clinical outcome in hepatocellular carcinoma. Cancer Med. 9, 2868–2878. doi: 10.1002/cam4.2921
Theodoraki, M.-N., Yerneni, S. S., Hoffmann, T. K., Gooding, W. E., and Whiteside, T. L. (2018). Clinical significance of PD-L1 exosomes in plasma of head and neck cancer patients. Clin. Cancer Res. 24, 896–905. doi: 10.1158/1078-0432.CCR-17-2664
Tokunaga, R., Zhang, W., Naseem, M., Puccini, A., Berger, M. D., Soni, S., et al. (2018). CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation - a target for novel cancer therapy. Cancer Treat. Rev. 63, 40–47. doi: 10.1016/j.ctrv.2017.11.007
Wang, N., Liu, W., Zheng, Y., Wang, S., Yang, B., Li, M., et al. (2018). CXCL1 derived from tumor-associated macrophages promotes breast cancer metastasis via activating NF-κB/SOX4 signaling. Cell Death Dis. 9:880. doi: 10.1038/s41419-018-0876-3
Wang, Z., Feng, X., Molinolo, A. A., Martin, D., Vitale-Cross, L., Nohata, N., et al. (2019). 4E-BP1 is a tumor suppressor protein reactivated by mTOR inhibition in head and neck cancer. Cancer Res. 79, 1438–1450. doi: 10.1158/0008-5472.CAN-18-1220
Weeraratna, A. T., Jiang, Y., Hostetter, G., Rosenblatt, K., Duray, P., Bittner, M., et al. (2002). Wnt5a signaling directly affects cell motility and invasion of metastatic melanoma. Cancer Cell 1, 279–288.
Wichmann, G., Rosolowski, M., Krohn, K., Kreuz, M., Boehm, A., Reiche, A., et al. (2015). The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer. Int. J. Cancer 137, 2846–2857. doi: 10.1002/ijc.29649
Wu, J., Zhao, Y., Zhang, J., Wu, Q., and Wang, W. (2019). Development and validation of an immune-related gene pairs signature in colorectal cancer. Oncoimmunology 8:1596715. doi: 10.1080/2162402X.2019.1596715
Zaenker, P., Gray, E. S., and Ziman, M. R. (2016). Autoantibody production in cancer–the humoral immune response toward autologous antigens in cancer patients. Autoimmun. Rev. 15, 477–483. doi: 10.1016/j.autrev.2016.01.017
Zhang, M., Zhu, K., Pu, H., Wang, Z., Zhao, H., Zhang, J., et al. (2019). An immune-related signature predicts survival in patients with lung adenocarcinoma. Front. Oncol. 9:1314. doi: 10.3389/fonc.2019.01314
Zhang, S., Li, G., Liu, C., Lu, S., Jing, Q., Chen, X., et al. (2020). miR-30e-5p represses angiogenesis and metastasis by directly targeting AEG-1 in squamous cell carcinoma of the head and neck. Cancer Sci. 111, 356–368. doi: 10.1111/cas.14259
Keywords: immune-related gene pairs, head and neck squamous cell carcinomas, the Cancer Genome Atlas, gene expression omnibus, prognosis, tumor immunology
Citation: He J, Fang X and Han M (2021) Discovery and Verification of an Immune-Related Gene Pairs Signature for Predicting Prognosis in Head and Neck Squamous Cell Carcinoma. Front. Genet. 12:654657. doi: 10.3389/fgene.2021.654657
Received: 17 January 2021; Accepted: 29 April 2021;
Published: 24 May 2021.
Edited by:Zhaolan Hu, Stanford University, United States
Reviewed by:Xi Chen, University of Maryland, United States
Yuan Zhou, George Washington University, United States
Yuanyuan Zhang, Johns Hopkins Medicine, United States
Copyright © 2021 He, Fang and Han. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mingming Han, firstname.lastname@example.org
†These authors have contributed equally to this work and share first authorship