Ruijie Ming
Xiangrui Li
Enhao Wang
Jiahui Wei
Bo Liu
Peng Zhou
Wenting Yu
Shimin Zong
Hongjun Xiao- Department of Otorhinolaryngology, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
Purpose: This study aimed to construct a prognostic signature consisting of immune-related RNA-binding proteins (RBPs) to predict the prognosis of patients with head and neck squamous cell carcinoma (HNSCC) effectively.
Methods: The transcriptome and clinical data of HNSCC were downloaded from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. First, we ascertained the immunological differences in HNSCC, through single-sample gene set enrichment analysis, stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE), and cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) deconvolution algorithm. Then we used univariate proportional hazards (Cox) regression analysis and least absolute shrinkage and selection operator (LASSO) Cox regression analysis to screen immune-related RBPs and acquire the risk score of each sample. Subsequently, we further investigated the difference in prognosis, immune status, and tumor mutation burden in high- and low-risk groups. Finally, the efficacy of immunotherapy was measured by the tumor immune dysfunction and exclusion (TIDE) score.
Results: We derived 15 immune-related RBPs, including FRMD4A, ASNS, RAB11FIP1, FAM120C, CFLAR, CTTN, PLEKHO1, SELENBP1, CHCHD2, NPM3, ATP2A3, CFDP1, IGF2BP2, NQO1, and DENND2D. There were significant differences in the prognoses of patients in the high- and low-risk groups in the training set (p < 0.001) and the validation set (p < 0.01). Furthermore, there were statistical differences between the high-risk group and low-risk group in immune cell infiltration and pathway and tumor mutation load (p < 0.001). In the end, we found that patients in the low-risk group were more sensitive to immunotherapy (p < 0.001), and then we screened 14 small-molecule chemotherapeutics with higher sensitivity to the high-risk group (p < 0.001).
Conclusion: The study constructed a prognostic signature of HNSCC, which might guide clinical immunotherapy in the future.
Introduction
Head and neck squamous carcinoma (HNSCC), which has a mortality rate of 50.5%, is one of the most common tumors, accounting for 3.6% of malignant tumors (1). HNSCC is a histologically and genetically heterogeneous disease that originates from a variety of anatomical parts, including the oral cavity, tongue, salivary glands, nasopharynx, and larynx (2). Smoking, drinking, and human papillomavirus infection are the main causes of HNSCC (3). Patients with HNSCC often experience cervical lymph node metastasis, local recurrence, and resistance to radiotherapy and chemotherapy (4).
At present, the treatment strategy for HNSCC patients is still based on tumor location and disease stage, not tumor biology. Many biomolecular markers, such as proteins, DNA, RNA, and microRNA, have been proposed to detect primary and secondary malignancies in the initial stages of the disease, but the above indicators are still very limited in terms of prognostic assessment and optimization of treatment options. In order to improve the treatment outcome of HNSCC, a clinically useful method is urgently needed to identify the risk of HNSCC and judge the effectiveness of adjuvant therapy.
The tumor microenvironment (TME) plays a vital role in the occurrence, progression, and treatment response of tumors. TME includes proliferating tumor cells, tumor stroma, blood vessels, cancer-related fibroblasts, infiltrating inflammatory cells, and various related signal molecules (5, 6). In the microenvironment of HNSCC, immune cells and mesenchymal cells, as the two main non-tumor components, have caused a large number of inflammatory reactions (7). Since HNSCC is an immunosuppressive disease, immune checkpoint inhibitors have emerged as a new treatment option (8). The basic principle of immunotherapy is to block the immunosuppressive effect of immune checkpoints while activating the endogenous immune system, thus increasing the number and cytotoxicity of T cells, which is beneficial to attack tumor cells (9). Consequently, it would be valuable to investigate the role of immune cells and their regulators in the TME of HNSCC.
From the nucleus to the peripheral cytoplasm, RNA-binding proteins (RBPs) play a vital role in the post-transcriptional regulation of genes (10). RBPs are able to affect pre-mRNA processing, transport and localization, mRNA stability/degradation, and translation (11). In a variety of tumors, some RBPs were found to be dysfunctional and aberrantly regulated (12, 13). Meanwhile, RBPs are important components of the immune system, which respond quickly to inflammatory mediators and in modulating inflammatory responses (14). Considering the important role of RBPs in immunity, it is necessary to explore the relationship between RBPs and HNSCC.
This study aimed to develop a prognostic prediction model for HNSCC based on immune-related RBPs. First, we classified HNSCC patients into two immune phenotypes based on the enrichment fraction of immune cells, then screened for differentially expressed RBPs in two immune phenotypes, and defined them as immune-related RBPs. Through univariate proportional hazards (Cox) regression analysis and least absolute shrinkage and selection operator (LASSO) Cox regression analysis, we identified immune-related RBPs related to prognosis and then constructed a risk model for patients with HNSCC. Based on the validation of the prognostic relevance and predictive capacity of the risk model, we further analyzed the infiltrating immune cells and immune-related pathways, somatic mutations, copy number variations (CNVs), the efficacy of immunotherapy, and sensitivity of chemotherapeutic agents in patients with HNSCC. The results showed that the risk model consisting of immune-related RBPs can effectively differentiate the clinical outcomes and show superiority in predicting the prognosis of patients with HNSCC.
Methods
Data Access
The transcriptome data in the fragment per kilobase million (FPKM) format and clinical data of 499 patients with HNSCC were downloaded from The Cancer Genome Atlas (TCGA) as the training set (https://portal.gdc.cancer.gov) (15) and downloaded the transcriptome data and clinical data of 97 HNSCC samples from the GSE41613 dataset of the Gene Expression Omnibus (GEO) database for validation (https://www.ncbi.nlm.nih.gov/geo/) (16). The data of somatic mutation and CNVs of patients with HNSCC were downloaded from UCSC (http://xena.ucsc.edu/) (17). The gene list of RBPs was collected from Gerstberger (10), SONAR (18), GO: RNA binding (19), poly(A) RBPs (20–24), CARIC (25), and XRNAX (26).
Immunophenotyping Based on Single-Sample Gene Set Enrichment Analysis
Single-sample gene set enrichment analysis (ssGSEA) is an algorithm based on rank ordering, which can calculate the degree of enrichment of a single sample in a given gene set (27). On this basis, the enrichment scores of immune cells and some related immune processes were calculated through the GSEA program (28, 29) and then quantified through the default parameters of the “Gene Set Variation Analysis (GSVA)” R package (30). Subsequently, the “ConsensusClusterPlus” R package was used to co-cluster the infiltration levels of 23 types of immune cells in HNSCC samples from TCGA to identify and distinguish immune subtypes (31). In the cumulative distribution function (CDF), the K value with the largest area under the curve was selected as 2, and so the HNSCC samples were divided into two types (31). The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm was utilized to calculate the immune score, stromal score, ESTIMATE score, and tumor purity (32). The immune cell infiltration calculated by the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) deconvolution algorithm was used to verify the immune difference between the two types (33). Finally, the GSEA program was used to compare the differences in pathway enrichment between the above immunotypes from the Kyoto Encyclopedia of Genes and Genomes (KEGG) (34).
Construction and Validation of Risk Model
The “limma” R package was used to distinguish RBPs with different expressions between immunotypes. With a 1.4-fold difference and corrected p less than 0.05 as the screening conditions, 238 immune-related RBPs were obtained. Subsequently, 47 immune-related RBPs associated with prognosis were obtained through univariate proportional hazards regression (p < 0.05). The “glmnet” package was then utilized to perform LASSO Cox regression analysis (35). After 1,000 times of cross-validation, 15 immune-related RBPs and the correlation coefficients of the corresponding risk genes were obtained to construct a risk model at the same time. , in which Expi is the expression of each risk gene and Coefi is its correlation coefficient. All patients were divided into a high-risk group and a low-risk group characterized by the median risk score of patients with HNSCC in the training set. The Kaplan–Meier curves were used to compare the overall survival (OS) difference of patients in the high- and low-risk groups. Receiver operating characteristic (ROC) curves were generated to evaluate the effectiveness and accuracy of the risk score in predicting the prognosis of patients with HNSCC. Next, the “ggExtra” R package was used to calculate the correlation between the risk score and the OS of patients with HNSCC. The independent correlation between the risk score and the prognosis of patients with HNSCC was then evaluated by univariate and multivariate proportional hazards regression analyses. Subsequently, a nomogram that could predict the prognosis of individual patients with HNSCC was constructed based on the stage, T stage, N stage, and risk group of patients with HNSCC through the “rms” R package (36). The C index was then used to assess the ability of the nomogram to distinguish prognosis, and a calibration chart was drawn to evaluate the accuracy of the nomogram. In addition, GSEA and gene set variation analysis (GSVA) were used to compare the differences in KEGG pathway enrichment between risk groups.
Analysis of Somatic Mutation and Copy Number Variations
The tumor mutation burden (TMB) of HNSCC samples from TCGA was analyzed through the “maftools” R package (37). The differences in TMB between the high- and low-risk groups were compared and showed the top 20 genes with the highest mutation rate and their mutation types in the high- and low-risk groups. Then the impact of TMB on the OS of patients with HNSCC was evaluated through the Kaplan–Meier survival curves. After that, gistic 2.0 was used to detect significant copy number amplification or deletion (38). In the end, the CNVs of 22 pairs of autosomes between the high- and low-risk groups were compared and showed the top 20 genes with most CNVs and their variation types.
Prediction of the Curative Effect of Immunotherapy and Chemotherapy
Tumor immune dysfunction and exclusion (TIDE) (http://tide.dfci.harvard.edu/) was used to calculate the TIDE score, which was reported to be able to predict the response of patients with a malignant tumor to immunotherapy (39). On the other hand, the “pRRophetic” R package was used to compare the half-maximal inhibitory concentration (IC50) differences of some common small-molecule chemotherapeutics between the high- and low-risk groups and screened out chemotherapeutics that may have better efficacy for patients in the high-risk group (40).
Statistical Analysis
All statistical analyses were based on R 4.0.4 software (https://www.r-project.org/). Categorical variables were tested by the chi-square test or Fisher’s exact test. The t-test or Wilcoxon test was performed on continuous variables. p < 0.05 was deemed statistically significant.
Results
Development and Validation of the Prognostic Model Based on Immunophenotyping of Head and Neck Squamous Cell Carcinoma
The flowchart of this research is shown in Figure 1. First, we obtained patient data from TCGA database and divided the patients into two groups according to differences in immune cells. The CIBERSORT deconvolution and ESTIMATE algorithm confirmed the difference in the immune microenvironment between the Sub1 and Sub2 groups (Figure 1A). After differential expression analysis, 238 immune-related RBPs were identified. Through univariate and LASSO Cox regression analysis, 15 immune-related RBPs related to prognosis were selected, and then the Kaplan–Meier curves showed the difference between the high- and low-expression immune-related RBPs groups (Figure 1B). Subsequently, we found that the risk score was significantly related to the OS of patients with HNSCC in training and validation sets, respectively (Figure 1C). In addition, the differences in immune cells and pathways between the high- and low-risk groups are further elaborated (Figure 1D). In terms of genes, we showed the differences in somatic mutation and CNVs (Figure 1E). In addition, we have also produced a nomogram combining the stage, T stage, N stage, and risk group to predict the prognosis (Figure 1F). In the end, the efficacy of immunotherapy was analyzed through the TIDE score, and the sensitivity of different risk groups to small-molecule chemotherapeutics was also revealed (Figure 1G).
Figure 1 Flowchart of this study. Two immune subtypes identified by single-sample gene set enrichment analysis (ssGSEA) and co-clustering analysis, and difference of infiltrating immune cells assessed by CIBERSORT deconvolution algorithm and ESTIMATE algorithm (A). Fifteen immune-related RNA-binding proteins (RBPs) screened out through “limma” package, univariate and least absolute shrinkage and selection operator (LASSO) Cox analysis, and the Kaplan–Meier curves for high- and low-expression immune-related RBP groups (B). Validation of the risk model composed of immune-related RBPs for prognosis in The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) database (C). Immune cell infiltration and pathways in high- and low-risk groups (D). Somatic mutation and copy number variations (CNVs) in high- and low-risk groups (E). Construction and calibration of prognosis nomogram (F). The differences of tumor immune dysfunction and exclusion (TIDE) score and sensitivity to chemotherapeutics of patients with head and neck squamous cell carcinoma (HNSCC) in high- and low-risk groups (G). *p < 0.05; **p < 0.01; ***p < 0.001.
Based on the transcriptome data of TCGA HNSCC, we evaluated and quantify 23 kinds of immune cells and 15 kinds of immune processes by ssGSEA. After that, co-clustering analysis was used to distinguish the infiltration of 23 immune cells in the HNSCC samples in TCGA. When K = 2, the CDF curve had the largest area under the curve, so all samples were divided into two types (Sub1 and Sub2) (Figure S1A). Among them, there were 271 cases in the Sub1 group and 228 cases in the Sub2 group. It was worth mentioning that the immune cells and pathways were more enriched in the Sub2 group than the Sub1 group (Figure 2A). Compared with the Sub2 group, the Sub2 group had lower immune score (Figure 2B, p < 0.001), lower stromal score (Figure 2C, p < 0.001), lower ESTIMATE score (Figure 2D, p < 0.001), and higher tumor purity (Figure 2E, p < 0.001). For the purpose of authenticating the difference between the two types, we used the CIBERSORT deconvolution algorithm and the ESTIMATE algorithm to calculate the infiltration of immune cells. Among the Sub1 group, M0 macrophages, activated dendritic cells, and mast cells infiltrated more, while in the Sub2 group, primitive B cells, plasma cells, CD8 T cells, activated CD4 memory T cells, follicular helper T cells, Treg cells, M1 macrophages, resting mast cells, and eosinophils infiltrated more (Figure 2F, p < 0.05). As far as the human leukocyte antigen (HLA) family is concerned, the expression of the Sub1 group is lower (Figure 2G, p < 0.001). Considering the rise of immune checkpoint inhibitor therapy, we also analyzed the differences between immune checkpoints. The expressions of checkpoint LAG3, PDCD1, HAVCR2, CTLA4, and CD274 in the Sub2 group are extremely higher than that in the Sub1 group (Figure 2H, p < 0.001). In addition, as the result of pathway enrichment shows, there was more immune-related pathway enrichment in the Sub2 group, such as cytokine receptor interaction, chemokine signaling pathway, JAK-STAT signaling pathway, cell adhesion molecules cams, toll-like receptor signaling pathway, and natural killer cell-mediated cytotoxicity (Figure 2I, p < 0.001). It was worth noting that the Kaplan–Meier curves showed a better prognosis of the Sub2 group than that of the Sub1 group (Figure 2J, p = 0.007).
Figure 2 Immune subtypes of head and neck squamous cell carcinoma (HNSCC) were identified based on the tumor-infiltrating immune cells. Heatmap of single-sample gene set enrichment analysis (ssGSEA) scores for Sub1 group (n = 271) and Sub2 group (n = 228) (A). Comparison of immune score (B), stromal score (C), ESTIMATE score (D), and tumor purity (E) between Sub1 and Sub2 groups. Difference of immune cell infiltration between Sub1 and Sub2 groups (F). The expressions of HLA family genes in Sub1 and Sub2 groups (G). The discrepancy of immune checkpoint genes between Sub1 and Sub2 groups, including LAG3, PDCD1, HAVCR2, CTLA4, and CD274 (H). The divergence of enrichment pathways between Sub1 and Sub2 groups (I). Kaplan–Meier curves of Sub1 and Sub2 groups (J). ***p < 0.001.
Construction and Validation of Risk Model
We screened 238 immune-related RBPs through the “limma” R package (Figure 3A). Among these 238 immune-related RBPs, most of them were highly expressed in the Sub2 group, and the others were highly expressed in the Sub1 group (Figure 3B). Subsequently, 47 prognostic-related immune-related RBPs were obtained through univariate proportional hazards regression (Figure 3C, p < 0.05). In order to avoid overfitting, we then used LASSO Cox regression analysis and cross-validated 1,000 times to obtain 15 immune-related RBPs and the correlation coefficients of their corresponding risk genes (Figures 3D, E, Table S1). The risk model was thus constructed:
Figure 3 Construction of risk model for prognosis in patients with head and neck squamous cell carcinoma (HNSCC). Volcano plot exhibiting the differentially expressed immune-related RNA-binding proteins (RBPs) between Sub1 group (n = 271) and Sub2 group (n = 228) in HNSCC (A). Heatmap of differentially expressed immune-related RBPs in Sub1 and Sub2 groups (B). The result of univariate Cox analysis (C) and least absolute shrinkage and selection operator (LASSO) Cox analysis (D, E).
The negative correlation coefficient indicated that the expression of the gene was beneficial to the prognosis, and the positive value indicated no benefit or even hindrance.
We assigned TCGA data as the training set and GEO data as the validation set. According to the median risk score of patients with HNSCC in TCGA, all patients were divided into the high-risk group and low-risk group (Figures 4A, H). In the training and validation sets, the mortality of patients in the high-risk group was higher than that in the low-risk group (Figures 4B, I). In the high-risk group, ASNS, CTTN, CHCHD2, NPM3, CFDP1, IGF2BP2, and NQO1 were expressed higher, while in the low-risk group, there were higher expressions of FRMD4A, RAB11FIP1, FAM120C, CFLAR, PLEKHO1, SELENBP1, ATP2A3, and DENND2D (Figures 4C, J). The OS was negatively correlated with the risk score, which meant the OS of patients with HNSCC gradually decreases as the risk score increased (Figures 4D, K). The area under the ROC (AUC) of the risk score of the training set was 0.60 (1 year), 0.70 (3 years), and 0.64 (5 years) (Figure 4E). In contrast, the AUC of the validation set was 0.63 (1 year), 0.63 (3 years), and 0.64 (5 years) (Figure 4L). The Kaplan–Meier curves also indicated that the high-risk group had a poor prognosis (Figures 4F, G, p < 0.01).
Figure 4 Application and validation of the risk model for prognosis. Samples in The Cancer Genome Atlas (TCGA) dataset were designated as training set, and samples in Gene Expression Omnibus (GEO) dataset were designated as validation set. On basis of the mean risk score of samples in training set, patients were divided into high-risk (red dot) and low-risk (green dot) groups. Distribution of the risk scores of the patients in training set (A). Distribution of survival time of patients in training set (B). The heatmap depicting the expression difference of 15 immune-related RNA-binding proteins (RBPs) between the high-risk group and the low-risk group in training set (C). Correlation between overall survival and risk score in training set (D). ROC curves of risk score for predicting 1, 3, and 5 years of overall survival in training set (E). Kaplan–Meier curves of high- and low-risk groups in training set (F). Distribution of the risk scores of the samples in validation set (H). Distribution of survival time of samples in validation set (I). The heatmap showing the expression patterns of 15 immune-related RBPs between the high- and low-risk groups in validation set (J). Correlation between overall survival and risk score in validation set (K). Receiver operating characteristic (ROC) curves of risk score for predicting 1, 3, and 5 years of overall survival in validation set (L). Kaplan–Meier curves of high- and low-risk groups in validation set (G).
In order to verify the validity and independence of the risk score, we combined the clinical characteristics and pathological staging data from TCGA database to perform univariate and multivariate Cox regression analyses. Univariate Cox analysis showed that age (p < 0.05), stage (p < 0.001), T (p < 0.01), N (p < 0.001), and risk score (p < 0.001) were significantly related to the prognosis (Figure 5A, Table S2). Multivariate analysis indicated that age (p < 0.01), N (p < 0.05), and risk score (p < 0.001) were significantly correlated with the prognosis (Figure 5B, Table S3). This implied that our risk model based on immune-related RBPs could be used as independent and effective indicators for the prognosis of patients with HNSCC.
Figure 5 Independence of risk score and construction of nomogram consisting of risk score and clinicopathological characteristics. Univariate Cox regression analysis was used to validate whether age, gender, grade, stage, T, N, and risk score had an independent influence on prognosis (A). Multivariate Cox regression analysis was used to validate whether age, gender, grade, stage, T, N, and risk score had independent influence on prognosis (B). Construction of integrated nomogram to predict survival in head and neck squamous cell carcinoma (HNSCC) (C). Calibration curve for predicting 1, 3, and 5 years of overall survival (D). *p < 0.05, **p < 0.01.
In addition, we combined the stage (I–II and III–IV), T stage (T1–2 and T3–4), N stage (N0 and N1–3), and risk group (low and high) to construct 1-, 3-, and 5-year prognostic nomogram models (Figure 5C), which could guide clinical judgment more conveniently and effectively. For example, when an 80-year-old patient in a low-risk group is stage III–IV, T3–4, and N1–3, he would get a score of 375, which means that the probability of his survival time at less than 1 year, less than 3 years, and less than 5 years is 0.268, 0.574, and 0.698, respectively. The following calibration chart showed the difference between the OS predicted by this nomogram and the actual OS from TCGA database and suggested that the nomogram had certain accuracy (Figure 5D).
Finally, we evaluated the relationship between each of the 15 immune-related RBP genes in the model and the OS of patients with HNSCC. Patients with high expressions of ASNS, IGF2BP2, CFDP1, CHCHD2, CTTN, NPM3, and NQO1 have poor OS, while patients with high expressions of FRMD4A, FAM120C, ATP2A3, PLEKHO1, RAB11FIP1, DENND2D, CFLAR, and SELENBP1 have a better OS (Figures 6A–O, p < 0.05).
Figure 6 Validation of each immune-related RNA-binding protein (RBP) in the risk model. Kaplan–Meier curves showing the differences of overall survival in high- and low-expression immune-related RBPs ASNS (A), IGF2BP2 (B), CFDP1 (C), CHCHD2 (D), CTTN (E), NPM3 (F), NQO1 (G), FRMD4A (H), FAM120C (I), ATP2A3 (J), PLEKHO1 (K), RAB11FIP1 (L), DENND2D (M), CFLAR (N), and SELENBP1 (O) between high-expression (blue) group and low-expression (yellow) group.
Exploration of the Immune Microenvironment
The established risk model was based on immune-related RBPs, so it was necessary to confirm whether the model was related to the immune microenvironment of HNSCC. CIRBERSORT results showed that the 15 immune-related RBPs in the model all had associated immune cells (Figure 7A). Through the ESTIMATE algorithm, we found that the immune score, stromal score, and ESTIMATE score were lower and that the tumor purity was higher in the high-risk group (Figures 7B–E, p < 0.001). Then, we compared the expressions of the HLA family, and most of them were lower in the high-risk group (Figure 7F, p < 0.05). Subsequently, the checkpoint expressions of PDCD1, CD274, CTLA4, HAVCR2, and LAG3 in the low-risk group were relatively high (Figures 7G–K, p < 0.001). Every immune checkpoint is negatively correlated with the risk score (Figures 7L–P, p < 0.001).
Figure 7 Immune landscape of patients with head and neck squamous cell carcinoma (HNSCC) in high- and low-risk groups. Correlation matrix of 15 immune-related RNA-binding proteins (RBPs) and infiltrating immune cells (A). Comparison of immune score (B), stromal score (C), ESTIMATE score (D), and tumor purity (E). The differential expressions of HLA family genes in patients with HNSCC in high- and low-risk groups (F). The expression level of immune checkpoint genes PDCD1 (G), CD274 (H), CTLA4 (I), HAVCR2 (J), and LAG3 (K) in low-risk group and high-risk group. The correlation between risk score and immune checkpoints PDCD1 (L), CD274 (M), CTLA4 (N), HAVCR2 (O), and LAG3 (P). “ns” means p ≥ 0.05, *p < 0.05, **p < 0.01, ***p < 0.001.
Analysis of Somatic Mutation and Copy Number Variations
In the high-risk group and the low-risk group, the genes with the highest mutations are TP53, TTN, FAT1, and CDKN2A. Moreover, there are more mutations of TP53, FAT1, CDKN2A, NOTCH1, SYNE1, and NSD1 in the high-risk group, and the mutation rate of PIK3CA is higher in the low-risk group (Figures 8A, B). TMB is higher in the high-risk group (Figure 8C, p < 0.001). The prognosis of patients with high TMB was significantly worse than that of patients with low TMB (Figure 8D). Considering that the risk score was an independent prognostic factor, we evaluated the superimposed influence of TMB and risk score. The prognosis in descending order is the low-mutation and low-risk group, the high-mutation and low-risk group, the low-mutation and high-risk group, and the high-mutation and high-risk group (Figure 8E, p < 0.001).
Figure 8 Somatic mutation and copy number variations (CNVs) in high- and low-risk groups. Heatmap of somatic mutations in high-risk group (A) and low-risk group (B). The difference of tumor mutation burden between high- and low-risk groups (C). Kaplan–Meier curves showing the differences in high- and low-tumor mutation burden (TMB) groups (D). Kaplan–Meier curves revealing the differences in high-TMB and high-risk group, high-TMB and low-risk group, low-TMB and high-risk group, and low-TMB and low-risk group (E). Amplification and deletion of copy number in the high-risk group (inner) and low-risk group (outer) (F). The 20 genes with maximum CNVs in high-risk group, and the percentage meaning the proportion of patients with head and neck squamous cell carcinoma (HNSCC) who suffered gene deletion (blue) or amplification (red) in high-risk group (G). Top 20 genes with maximum CNVs in low-risk group, and the percentage representing the ratio of patients with HNSCC who suffered gene deletion (blue) or amplification (red) in low-risk group (H).
Extensive copy number amplification was detected in 22 pairs of autosomes in all two groups. In the low-risk group, high-frequency deletion regions were found on chromosomes 3 and 13, and high-frequency amplification regions were found on chromosome 8 (Figure 8F). In the high-risk group, CNVs analysis indicated the following most relevant genes: MIR7641-2|chr8, CASC19, CCAT1, CASC21, CASC8, CCAT2, POU5FIB, and CSMD1 (Figure 8G). Among them, gene CSMD1 had a significant copy number deletion (Figure 8G). On the other hand, the five most correlative genes in the low-risk group included NAALADL2, TP63, LINC01206, TPRG1, and TPRG1-AS2 (Figure 8H).
GSEA (Figure 9A) and GSVA (Figure 9B) revealed the differences in pathway enrichment between the high- and low-risk groups. Most of the pathways enriched in the low-risk group were associated with immune responses, which may be involved in immune-related RBPs, including Fc gamma R-mediated phagocytosis, B-cell receptor signaling pathway, T-cell receptor signaling pathway, autoimmune thyroid disease, cell adhesion molecules cams, cytokine–cytokine receptor interaction, leukocyte transendothelial migration, and natural killer cell-mediated cytotoxicity.
Figure 9 Enrichment signaling pathways of different risk groups. The pathway enrichment of gene set variation analysis (GSVA) between the low- and high-risk groups (A). The pathway enrichment of gene set enrichment analysis (GSEA) between the low- and high-risk groups (B).
Prediction of the Efficacy of Immunotherapy and Chemotherapy
We used the TIDE score to predict the immunotherapy response of patients with HNSCC to immunotherapy. It could be briefly described that the higher the TIDE score, the higher the likelihood of immune dysfunction or evasion, and the less likely the patient will benefit from immune checkpoint inhibitors. As a result, the TIDE score of the high-risk group was significantly higher than that of the low-risk group, which means that immunotherapy is less effective in the high-risk group (Figure 10A, p < 0.001).
Figure 10 The value of the risk model in predicting the efficacy of immunotherapy and chemotherapy. The score of tumor immune dysfunction and exclusion of patients with head and neck squamous cell carcinoma (HNSCC) in high- and low-risk groups (A). The box plots of the estimated IC50 for bosutinib (B), bryostatin.1 (C), camptothecin (D), cytarabine (E), docetaxel (F), doxorubicin (G), erlotinib (H), gefitinib (I), gemcitabine (J), lapatinib (K), paclitaxel (L), parthenolide (M), sorafenib (N), and thapsigargin (O).
In addition, we screened out 14 small-molecule chemotherapeutics that may be more effective for patients with HNSCC in the high-risk group. The IC50 represents the concentration of an inhibitor that is required for 50% inhibition of carcinoma cells. A lower IC50 value means better drug sensitivity. Patients in the high-risk group were more sensitive to bosutinib, bryostatin.1, camptothecin, cytarabine, docetaxel, doxorubicin, erlotinib, gefitinib, gemcitabine, lapatinib, paclitaxel, parthenolide, sorafenib, and thapsigargin (Figures 10B–O, p < 0.001).
Discussion
Immunotherapy has become an effective method for treating malignant tumors (41). Furthermore, immunosuppressant therapy has made important progress in the treatment of patients with HNSCC (42). Nevertheless, it cannot be ignored that only a limited one-third of patients respond to immunotherapy in most types of tumors (43). Further studies of immune-related RBPs in HNSCC may provide new ways to improve the clinical prognosis of patients. At present, there is an urgent need for an accurate and operational prognostic evaluation model for HNSCC in clinical practice. Based on TCGA and GEO databases and a variety of algorithms starting with ssGSEA, our study established a new model for predicting immune response, efficacy of conventional chemotherapy and immunotherapy, and individual outcome.
There are many kinds of myeloid immune cells in the HNSCC microenvironment that have a unique immune profile prior to treatment (44). In this study, we retrospectively analyzed the transcriptomic data of 499 HNSCC patients in TCGA database and further classified them into Sub1 and Sub2 on the basis of differences in immune cell infiltration. Regarding the infiltrating immune cells in the Sub1 group, M0 macrophages infiltrated more, while in the Sub2 group, there was more infiltration of naive B cells, plasma cells, T cells, and M1 macrophages. Compared with the Sub1 group, the Sub2 group had higher immune, stromal, and ESTIMATE scores but lower tumor purity, and its prognosis was significantly better than that of the Sub1 group. In addition, the expressions of the HLA family were significantly lower in the Sub1 group, which assisted tumor cells to escape the immune system (45). The immune-related pathways in the Sub2 group were more abundant.
Some RBPs are able to rapidly react to inflammatory mediators and regulate the reprogramming of immune cells to tumor-associated phenotypes (12). After recognizing the difference in RBP expressions between the Sub1 and Sub2 groups, we constructed a risk model containing 15 immune-related RBPs through univariate and LASSO Cox regression analysis. Then, according to the risk score calculated by the above model, patients with HNSCC were divided into low- and high-risk groups. Patients in the high-risk group had poorer clinical outcomes. The model even had good validity and stability in determining the prognosis at 1, 3, and 5 years, which was further confirmed in the GEO database. After confirming the risk score as an independent prognostic factor, we constructed a prognostic nomogram model according to the staging, T, N, and age of patients, which was also accurate in predicting OS at 1, 3, and 5 years.
The relationship between the prognostic characteristics of 15 immune-related RBPs and the immune microenvironment has also been investigated. In contrast with the high-risk group, the low-risk group had lower tumor purity and higher immune score, stromal score, and ESTIMATE score. In the high-risk group, M0 macrophages, activated dendritic cells, and mast cells infiltrated more, but the expression of the HLA family decreased. In the low-risk group, there were more M1 macrophages, naive B cells, CD4 memory and CD8 T cells, plasma cells, and eosinophils. This is similar to previous reports suggesting that exhausted immunity with lower survival is characterized by enrichment of stromal activation and anti-inflammatory M2 macrophage, whereas enhanced immunity associated with better prognosis is characterized by M1 macrophages providing stronger pro-inflammatory signaling, enhanced cytolytic activity, and massive lymphocyte infiltration (7). The activation of M1 macrophages is beneficial to patients because it can induce acute inflammation secreting tumor-killing molecules such as tumor necrosis factor α (TNFα) (46). On the other hand, if acute inflammation is not controlled, differentiation of M2 macrophages facilitates chronic inflammation, promoting tumor cell growth, angiogenesis, fibrosis, and immunosuppression (47), which is certainly harmful to patients. Both B cells and partial T cells also contribute to the prognosis of patients with HNSCC. As reported by Norouzian et al., the composition of B-cell subpopulations changes in TME of HNSCC, and the B cells with atypical memory and regulatory phenotype are significantly related to favorable prognostic (48). Notably, the high abundance of tumor-infiltrating lymphocyte B and high density of direct B-cell/CD8+ T-cell interactions predict a better outcome (49). Dense T-cell infiltration, especially cytotoxic CD8 T cells, represents superior antitumor ability (50, 51).
Based on the risk score, we further elaborated on TMB, somatic mutations, and CNVs. The high-risk group had a higher TMB, which implied a higher mortality rate. Mutations in TP53 were overwhelmingly predominant in both groups and were more frequent in the high-risk group than in the low-risk group (70% vs. 55%). As previously reported by Lawrence, TP53 mutations and CDKN2 inactivation are intimately involved in HNSCC (52). Remarkably, TP53 mutations are common and associated with a poor prognosis in patients with HNSCC (53). MIR7641 is highly expressed in the exosomes of metastatic tumor cells and can enhance the proliferation, migration, and invasion of recipient tumor cells (54, 55). Cub and Sushi Multiple Domains-1 (CSMD1) acts as a tumor suppressor, whose low expression promotes the invasion of HNSCC and gastric tumor (56, 57) and is also correlated with a poorer prognosis of HNSCC (58).
The efficacy of immunotherapy has been reported to be generally superior to that of conventional chemotherapy (42). Besides, the combined application of PD-1/PD-L1 inhibitors and platinum drugs also shows positive therapeutic potential (59). This implies that the exploration of medication regimens is potentially valuable. On the one hand, the immune checkpoint expressions of LAG3, PDCD1, HAVCR2, CTLA4, and CD274 increased in the low-risk group. The efficacy of corresponding immune checkpoint inhibitors is better for the low-risk group but the opposite for the high-risk group. On the other hand, the low-risk group has a lower TIDE score, which means that the lower TIDE score is related to a better curative effect. TIDE can be used to identify two mechanisms of tumor immune escape: inducing T-cell immunotherapeutic dysfunction in tumors with high infiltration of cytotoxic T lymphocytes (CTLs) and blocking T-cell infiltration in tumors with low CTL in TME (39). In our study, the low-risk group had more infiltration of CTLs, so they would respond better to immunotherapy, due to better recovery from T-cell dysfunction. The high-risk group had less CTL infiltration, so they would benefit less from immunotherapy, which may be due to T-cell repulsion. In short, the low-risk group will benefit more from immunotherapy. Furthermore, we screened out 15 chemotherapy drugs that are more suitable for the high-risk group. The new model constructed by immune-related RBPs could serve as a new marker to help guide the selection of chemotherapeutic drugs and distinguish who would benefit more from antitumor immunotherapy.
Some of the RBPs in this article have been reported to have a practical relationship with tumorigenesis and progression. Asparagine synthetase (ASNS) catalyzes the synthesis of the nonessential amino acid asparagine, while ASNS knockdown significantly hinders cell proliferation (60). In other words, stable ASNS gene expression guarantees the growth of tumor cells. Cortactin (CTTN) gene encodes a protein, cortacn, which plays an essential role in the migration of oral carcinoma cells by regulating filamentous actin and prominent structures on cell membranes (61). The high expression of CTTN was related to a poorer OS rate (62). Coiled-coil-helix-coiled-coil-helix domain-containing protein 2 (CHCHD2) as a small mitochondrial protein can regulate mitochondrial outer membrane permeabilization and is one of the negative regulators that mediate apoptosis (63). CHCHD2 indicates a poor prognosis and is overexpressed in hepatocellular carcinoma, breast tumor, non-small cell lung carcinoma, and renal cell carcinoma (64, 65). The loss of the human Cranio Facial Development Protein 1 (CFDP1) affects the dynamic changes of chromosomes and cell cycle progression (66). Moreover, some studies have confirmed that CFDP1 is a risk gene for pancreatic carcinoma (67, 68). High expression of insulin growth factor 2 mRNA binding protein 1 (IGF2BP1) is associated with a poor prognosis such as advanced clinical stage, increased tumor size, lymph node metastasis, and low survival rate of patients with HNSCC (69, 70). NAD(P)H quinone oxidoreductase (NQO1), a cytoplasmic enzyme that mediates the reduction of quinone substrates, is highly expressed in a multitude of tumors and can catalyze quinone drugs to poison tumor cells (71). NQO1 is considered a promising direct tumor target. For example, the drug β-lapachone, catalyzed by NQO1, triggers the innate perception of T cells in the TME, thereby enhancing antitumor capacity and even overcoming checkpoint blockade (72). Casein kinase 2-interacting protein-1 (CKIP-1, also known as PLEKHO1) inhibits tumor growth by causing inactivation of serine/threonine kinases and self-degradation of Smurf1, which is a potential oncogenic target in various tumor cells (73). Selenium binding protein 1(SELENBP1) is significantly downregulated in esophageal adenocarcinoma, ovarian tumor, and oral squamous cell carcinoma, but its overexpression can lead to incremental cellular senescence and apoptosis, as well as enhanced cytotoxicity of cisplatin (74–76). Three different genes (ATP2A1-3) encode the Ca2+-ATPases from the Sarco/endoplasmic reticulum (SERCA) to maintain calcium homeostasis between the cell cytoplasm and the endoplasmic reticulum, and they have been reported to downregulate transcription in gastric and colon tumors (77). In particular, ATP2A2 gene inactivation is closely related to oral squamous cell carcinoma (78). DENN/MADD domain-containing protein 2D (DENND2D) is less expressed in malignant tumors and is thought to contribute to the worsening prognosis and high recurrence rate (79–81). However, other RBPs may have a prospective regulatory impact on HNSCC. Some articles reported the relationship between genes FRMD4A and HNSCC. High expression of FRMD4A is associated with an increased risk of HNSCC recurrence, and the silencing of FRMD4A inhibits the growth and metastasis of human squamous cell carcinoma in skin and tongue metastases and reduces the proliferation and cell adhesion of squamous cell carcinoma (82, 83). Interestingly, in our study, patients with high expression of FRMD4A experienced a better prognosis (Figure 6H), which is worthy of further study. RAB proteins play the role of small GTPases in the regulation of vesicle and protein transport, membrane targeting, and fusion, and a group of them can actively or inversely regulate tumor cell generation, migration, and invasion (84). RAB11 affects the invasiveness of breast cancer cells (85). RAB11FIP1 is positively related to dendritic cells and CD4 T cells, and the low expression of RAB11FIP1 revealed a poor prognosis for lung adenocarcinoma (86). CASP8 and FADD-like apoptosis regulator (CFLAR), also known as c-FLICE-like inhibitory protein (c-FLIP), is a vital anti-apoptotic protein (87). Some studies have identified FLIP as an independent poor prognostic indicator for colorectal carcinoma, cervical carcinoma, and acute myeloid leukemia (88).
Although some studies have explored the association of RBPs with HNSCC (89, 90), our research has made further progress. On the basis of differentiated immunophenotyping, we take the lead in the screening of differentially expressed RBPs, which represents a more effective prognostic biomarker and a more accurate predictor of response to immunotherapy in different groups of patients. In general, the prognosis model system constructed based on the immune-related RBPs and clinical information of patients with HNSCC drew the landscape in the immune microenvironment of HNSCC and could effectively predict the prognosis of patients with HNSCC in the high- and low-risk groups. The nomogram based on this model is more helpful for predicting the clinical outcome of patients with HNSCC. Last but not least, the differences in immune checkpoints and TIDE scores between the high- and low-risk groups provide new ideas for the immunotherapy of patients with HNSCC.
Our study still has some limitations. First, we only used public databases to construct and verify the prognostic risk model, and we need to validate this model in subsequent clinical trials. Second, how the immune-related RBPs regulate immunity still needs to be verified by experiments in vitro and in vivo. Eventually, human papillomavirus is an independent prognostic factor for HNSCC, which is worth further stratified analysis.
Conclusion
In summary, the signature constructed by 15 immune-related RBPs could effectively predict the clinical outcome of patients with HNSCC. Subsequently, we demonstrated the immune landscape, TMB, CNVs, and efficacy of immunotherapy in different risk groups, which might guide clinical therapy.
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.
Author Contributions
RM and XL designed the study, performed the experiments, and plotted the data. EW collected original data. JW, BL, WY, and PZ drafted and edited the manuscript. HX and SZ reviewed the manuscript. HX and SZ supervised the project. HX and SZ funded the experiments for the study. All authors contributed to the article and approved the submitted version.
Funding
This study was supported by grants from the National Natural Science Foundation of China (grant numbers 81771002, 82071057).
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.
Acknowledgments
The authors gratefully acknowledge the data generated by all public datasets utilized in this research. The authors very much appreciate the efforts of all researchers in developing R package tools. The authors are very grateful to these authors for their selfless dedication to TCGA database and standardizing these data work from TCGA. The authors cherish these precious public database resources very much.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.795781/full#supplementary-material
References
1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492
2. Takeuchi T, Kawasaki H, Luce A, Cossu AM, Misso G, Scrima M, et al. Insight Toward the Microrna Profiling of Laryngeal Cancers: Biological Role and Clinical Impact. Int J Mol Sci (2020) 21(10):3693. doi: 10.3390/ijms21103693
3. Galot R, Le Tourneau C, Guigay J, Licitra L, Tinhofer I, Kong A, et al. Personalized Biomarker-Based Treatment Strategy for Patients With Squamous Cell Carcinoma of the Head and Neck: EORTC Position and Approach. Ann Oncol (2018) 29(12):2313–27. doi: 10.1093/annonc/mdy452
4. Li Q, Xu L, Li Y, Yang R, Qiao Q, Wang Y, et al. P2RY14 Is a Potential Biomarker of Tumor Microenvironment Immunomodulation and Favorable Prognosis in Patients With Head and Neck Cancer. Front Genet (2021) 12:670746. doi: 10.3389/fgene.2021.670746
5. Quail DF, Joyce JA. Microenvironmental Regulation of Tumor Progression and Metastasis. Nat Med (2013) 19(11):1423–37. doi: 10.1038/nm.3394
6. New J, Arnold L, Ananth M, Alvi S, Thornton M, Werner L, et al. Secretory Autophagy in Cancer-Associated Fibroblasts Promotes Head and Neck Cancer Progression and Offers a Novel Therapeutic Target. Cancer Res (2017) 77(23):6679–91. doi: 10.1158/0008-5472.CAN-17-1077
7. Chen YP, Wang YQ, Lv JW, Li YQ, Chua MLK, Le QT, et al. Identification and Validation of Novel Microenvironment-Based Immune Molecular Subgroups of Head and Neck Squamous Cell Carcinoma: Implications for Immunotherapy. Ann Oncol (2019) 30(1):68–75. doi: 10.1093/annonc/mdy470
8. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The Immune Landscape of Cancer. Immunity (2018) 48(4):812–830 e14. doi: 10.1016/j.immuni.2018.03.023
9. Ghanizada M, Jakobsen KK, Gronhoj C, von Buchwald C. The Effects of Checkpoint Inhibition on Head and Neck Squamous Cell Carcinoma: A Systematic Review. Oral Oncol (2019) 90:67–73. doi: 10.1016/j.oraloncology.2019.01.018
10. Gerstberger S, Hafner M, Tuschl T. A Census of Human RNA-Binding Proteins. Nat Rev Genet (2014) 15(12):829–45. doi: 10.1038/nrg3813
11. Dreyfuss G, Kim VN, Kataoka N. Messenger-RNA-Binding Proteins and the Messages They Carry. Nat Rev Mol Cell Biol (2002) 3(3):195–205. doi: 10.1038/nrm760
12. Kechavarzi B, Janga SC. Dissecting the Expression Landscape of RNA-Binding Proteins in Human Cancers. Genome Biol (2014) 15(1):R14. doi: 10.1186/gb-2014-15-1-r14
13. Sebestyen E, Singh B, Minana B, Pages A, Mateo F, Pujana MA, et al. Large-Scale Analysis of Genome and Transcriptome Alterations in Multiple Tumors Unveils Novel Cancer-Relevant Splicing Networks. Genome Res (2016) 26(6):732–44. doi: 10.1101/gr.199935.115
14. Turner M, Diaz-Munoz MD. RNA-Binding Proteins Control Gene Expression and Cell Fate in the Immune System. Nat Immunol (2018) 19(2):120–9. doi: 10.1038/s41590-017-0028-4
15. Mounir M, Lucchetta M, Silva TC, Olsen C, Bontempi G, Chen X, et al. New Functionalities in the Tcgabiolinks Package for the Study and Integration of Cancer Data From GDC and Gtex. PloS Comput Biol (2019) 15(3):e1006701. doi: 10.1371/journal.pcbi.1006701
16. Lohavanichbutr P, Mendez E, Holsinger FC, Rue TC, Zhang Y, Houck J, et al. A 13-Gene Signature Prognostic of HPV-Negative OSCC: Discovery and External Validation. Clin Cancer Res (2013) 19(5):1197–203. doi: 10.1158/1078-0432.CCR-12-2647
17. Goldman MJ, Craft B, Hastie M, Repecka K, McDade F, Kamath A, et al. Visualizing and Interpreting Cancer Genomics Data via the Xena Platform. Nat Biotechnol (2020) 38(6):675–8. doi: 10.1038/s41587-020-0546-8
18. Brannan KW, Jin W, Huelga SC, Banks CA, Gilmore JM, Florens L, et al. SONAR Discovers RNA-Binding Proteins From Analysis of Large-Scale Protein-Protein Interactomes. Mol Cell (2016) 64(2):282–93. doi: 10.1016/j.molcel.2016.09.003
19. Gene Ontology C, Blake JA, Dolan M, Drabkin H, Hill DP, Li N, et al. Gene Ontology Annotations and Resources. Nucleic Acids Res (2013) 41(Database issue):D530–5. doi: 10.1093/nar/gks1050
20. Conrad T, Albrecht AS, de Melo Costa VR, Sauer S, Meierhofer D, Orom UA. Serial Interactome Capture of the Human Cell Nucleus. Nat Commun (2016) 7:11212. doi: 10.1038/ncomms11212
21. Castello A, Fischer B, Frese CK, Horos R, Alleaume AM, Foehr S, et al. Comprehensive Identification of RNA-Binding Domains in Human Cells. Mol Cell (2016) 63(4):696–710. doi: 10.1016/j.molcel.2016.06.029
22. Castello A, Fischer B, Eichelbaum K, Horos R, Beckmann BM, Strein C, et al. Insights Into RNA Biology From an Atlas of Mammalian Mrna-Binding Proteins. Cell (2012) 149(6):1393–406. doi: 10.1016/j.cell.2012.04.031
23. Beckmann BM, Horos R, Fischer B, Castello A, Eichelbaum K, Alleaume AM, et al. The RNA-Binding Proteomes From Yeast to Man Harbour Conserved Enigmrbps. Nat Commun (2015) 6:10127. doi: 10.1038/ncomms10127
24. Baltz AG, Munschauer M, Schwanhausser B, Vasile A, Murakawa Y, Schueler M, et al. The Mrna-Bound Proteome and Its Global Occupancy Profile on Protein-Coding Transcripts. Mol Cell (2012) 46(5):674–90. doi: 10.1016/j.molcel.2012.05.021
25. Huang R, Han M, Meng L, Chen X. Transcriptome-Wide Discovery of Coding and Noncoding RNA-Binding Proteins. Proc Natl Acad Sci USA (2018) 115(17):E3879–87. doi: 10.1073/pnas.1718406115
26. Trendel J, Schwarzl T, Horos R, Prakash A, Bateman A, Hentze MW, et al. The Human RNA-Binding Proteome and Its Dynamics During Translational Arrest. Cell (2019) 176(1-2):391–403 e19. doi: 10.1016/j.cell.2018.11.004
27. Zuo S, Wei M, Wang S, Dong J, Wei J. Pan-Cancer Analysis of Immune Cell Infiltration Identifies a Prognostic Immune-Cell Characteristic Score (ICCS) in Lung Adenocarcinoma. Front Immunol (2020) 11:1218. doi: 10.3389/fimmu.2020.01218
28. Zhang L, Zhao Y, Dai Y, Cheng JN, Gong Z, Feng Y, et al. Immune Landscape of Colorectal Cancer Tumor Microenvironment From Different Primary Tumor Location. Front Immunol (2018) 9:1578. doi: 10.3389/fimmu.2018.01578
29. Finotello F, Trajanoski Z. Quantifying Tumor-Infiltrating Immune Cells From Transcriptomics Data. Cancer Immunol Immunother (2018) 67(7):1031–40. doi: 10.1007/s00262-018-2150-z
30. Hanzelmann S, Castelo R, Guinney J. GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7
31. Wilkerson MD, Hayes DN. Consensusclusterplus: A Class Discovery Tool With Confidence Assessments and Item Tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170
32. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
33. 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(5):453–7. doi: 10.1038/nmeth.3337
34. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102
35. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Software (2010) 33(1):1–22. doi: 10.18637/jss.v033.i01
36. Mao R, Chen Y, Xiong L, Liu Y, Zhang T. Identification of a Nomogram Based on an 8-Lncrna Signature as a Novel Diagnostic Biomarker for Head and Neck Squamous Cell Carcinoma. Aging (Albany NY) (2020) 12(20):20778–800. doi: 10.18632/aging.104014
37. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118
38. Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 Facilitates Sensitive and Confident Localization of the Targets of Focal Somatic Copy-Number Alteration in Human Cancers. Genome Biol (2011) 12(4):R41. doi: 10.1186/gb-2011-12-4-r41
39. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1
40. Geeleher P, Cox N, Huang RS. Prrophetic: An R Package for Prediction of Clinical Chemotherapeutic Response From Tumor Gene Expression Levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468
41. Emens LA, Ascierto PA, Darcy PK, Demaria S, Eggermont AMM, Redmond WL, et al. Cancer Immunotherapy: Opportunities and Challenges in the Rapidly Evolving Clinical Landscape. Eur J Cancer (2017) 81:116–29. doi: 10.1016/j.ejca.2017.01.035
42. Cohen EEW, Soulières D, Le Tourneau C, Dinis J, Licitra L, Ahn M-J, et al. Pembrolizumab Versus Methotrexate, Docetaxel, or Cetuximab for Recurrent or Metastatic Head-and-Neck Squamous Cell Carcinoma (KEYNOTE-040): A Randomised, Open-Label, Phase 3 Study. Lancet (2019) 393(10167):156–67. doi: 10.1016/s0140-6736(18)31999-8
43. Sharma P, Hu-Lieskovan S, Wargo JA, Ribas A. Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy. Cell (2017) 168(4):707–23. doi: 10.1016/j.cell.2017.01.017
44. 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(1):183–199 e9. doi: 10.1016/j.immuni.2019.11.014
45. Jhunjhunwala S, Hammer C, Delamarre L. Antigen Presentation in Cancer: Insights Into Tumour Immunogenicity and Immune Evasion. Nat Rev Cancer (2021) 21(5):298–312. doi: 10.1038/s41568-021-00339-z
46. Josephs DH, Bax HJ, Karagiannis SN. Tumour-Associated Macrophage Polarisation and Re-Education With Immunotherapy. Front Biosci (Elite Ed) (2015) 7:293–308. doi: 10.2741/E735
47. Ruffell B, Coussens LM. Macrophages and Therapeutic Resistance in Cancer. Cancer Cell (2015) 27(4):462–72. doi: 10.1016/j.ccell.2015.02.015
48. Norouzian M, Mehdipour F, Balouchi Anaraki S, Ashraf MJ, Khademi B, Ghaderi A. Atypical Memory and Regulatory B Cell Subsets in Tumor Draining Lymph Nodes of Head and Neck Squamous Cell Carcinoma Correlate With Good Prognostic Factors. Head Neck Pathol (2020) 14(3):645–56. doi: 10.1007/s12105-019-01095-1
49. Hladikova K, Koucky V, Boucek J, Laco J, Grega M, Hodek M, et al. Tumor-Infiltrating B Cells Affect the Progression of Oropharyngeal Squamous Cell Carcinoma via Cell-to-Cell Interactions With CD8(+) T Cells. J Immunother Cancer (2019) 7(1):261. doi: 10.1186/s40425-019-0726-6
50. Fridman WH, Zitvogel L, Sautes-Fridman C, Kroemer G. The Immune Contexture in Cancer Prognosis and Treatment. Nat Rev Clin Oncol (2017) 14(12):717–34. doi: 10.1038/nrclinonc.2017.101
51. Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The Prognostic Landscape of Genes and Infiltrating Immune Cells Across Human Cancers. Nat Med (2015) 21(8):938–45. doi: 10.1038/nm.3909
52. Cancer Genome Atlas N. Comprehensive Genomic Characterization of Head and Neck Squamous Cell Carcinomas. Nature (2015) 517(7536):576–82. doi: 10.1038/nature14129
53. Zhou G, Liu Z, Myers JN. TP53 Mutations in Head and Neck Squamous Cell Carcinoma and Their Impact on Disease Progression and Treatment Response. J Cell Biochem (2016) 117(12):2682–92. doi: 10.1002/jcb.25592
54. Liu Y, Hua F, Zhan Y, Yang Y, Xie J, Cheng Y, et al. Carcinoma Associated Fibroblasts Small Extracellular Vesicles With Low Mir-7641 Promotes Breast Cancer Stemness and Glycolysis by HIF-1alpha. Cell Death Discov (2021) 7(1):176. doi: 10.1038/s41420-021-00524-x
55. Shen S, Song Y, Zhao B, Xu Y, Ren X, Zhou Y, et al. Cancer-Derived Exosomal Mir-7641 Promotes Breast Cancer Progression and Metastasis. Cell Commun Signal (2021) 19(1):20. doi: 10.1186/s12964-020-00700-z
56. Chen XL, Hong LL, Wang KL, Liu X, Wang JL, Lei L, et al. Deregulation of CSMD1 Targeted by Microrna-10b Drives Gastric Cancer Progression Through the NF-Kappab Pathway. Int J Biol Sci (2019) 15(10):2075–86. doi: 10.7150/ijbs.23802
57. Veeramachaneni R, Walker T, Revil T, Weck A, Badescu D, O’Sullivan J, et al. Analysis of Head and Neck Carcinoma Progression Reveals Novel and Relevant Stage-Specific Changes Associated With Immortalisation and Malignancy. Sci Rep (2019) 9(1):11992. doi: 10.1038/s41598-019-48229-7
58. Jung AR, Eun YG, Lee YC, Noh JK, Kwon KH. Clinical Significance of CUB and Sushi Multiple Domains 1 Inactivation in Head and Neck Squamous Cell Carcinoma. Int J Mol Sci (2018) 19(12):3996. doi: 10.3390/ijms19123996
59. Xue Y, Gao S, Gou J, Yin T, He H, Wang Y, et al. Platinum-Based Chemotherapy in Combination With PD-1/PD-L1 Inhibitors: Preclinical and Clinical Studies and Mechanism of Action. Expert Opin Drug Deliv (2021) 18(2):187–203. doi: 10.1080/17425247.2021.1825376
60. Li H, Ning S, Ghandi M, Kryukov GV, Gopal S, Deik A, et al. The Landscape of Cancer Cell Line Metabolism. Nat Med (2019) 25(5):850–60. doi: 10.1038/s41591-019-0404-8
61. Ramos-Garcia P, Gonzalez-Moles MA, Gonzalez-Ruiz L, Ayen A, Ruiz-Avila I, Navarro-Trivino FJ, et al. An Update of Knowledge on Cortactin as a Metastatic Driver and Potential Therapeutic Target in Oral Squamous Cell Carcinoma. Oral Dis (2019) 25(4):949–71. doi: 10.1111/odi.12913
62. Horn D, Gross M, Dyckhoff G, Fuchs J, Grabe N, Weichert W, et al. Cortactin Expression: Association With Disease Progression and Survival in Oral Squamous Cell Carcinoma. Head Neck (2018) 40(12):2685–94. doi: 10.1002/hed.25515
63. Liu Y, Clegg HV, Leslie PL, Di J, Tollini LA, He Y, et al. CHCHD2 Inhibits Apoptosis by Interacting With Bcl-X L to Regulate Bax Activation. Cell Death Differ (2015) 22(6):1035–46. doi: 10.1038/cdd.2014.194
64. Yao Y, Su J, Zhao L, Li R, Liu K, Wang S. CHCHD2 Promotes Hepatocellular Carcinoma and Indicates Poor Prognosis of Hepatocellular Carcinoma Patients. J Cancer (2019) 10(27):6822–8. doi: 10.7150/jca.31158
65. Gundamaraju R, Lu W, Manikam R. Chchd2: The Power House’s Potential Prognostic Factor for Cancer? Front Cell Dev Biol (2020) 8:620816. doi: 10.3389/fcell.2020.620816
66. Messina G, Atterrato MT, Prozzillo Y, Piacentini L, Losada A, Dimitri P. The Human Cranio Facial Development Protein 1 (Cfdp1) Gene Encodes a Protein Required for the Maintenance of Higher-Order Chromatin Organization. Sci Rep (2017) 7:45022. doi: 10.1038/srep45022
67. Liu D, Zhou D, Sun Y, Zhu J, Ghoneim D, Wu C, et al. A Transcriptome-Wide Association Study Identifies Candidate Susceptibility Genes for Pancreatic Cancer Risk. Cancer Res (2020) 80(20):4346–54. doi: 10.1158/0008-5472.Can-20-1353
68. Zhong J, Jermusyk A, Wu L, Hoskins JW, Collins I, Mocci E, et al. A Transcriptome-Wide Association Study Identifies Novel Candidate Susceptibility Genes for Pancreatic Cancer. J Natl Cancer Inst (2020) 112(10):1003–12. doi: 10.1093/jnci/djz246
69. Chou CH, Chang CY, Lu HJ, Hsin MC, Chen MK, Huang HC, et al. IGF2BP2 Polymorphisms Are Associated With Clinical Characteristics and Development of Oral Cancer. Int J Mol Sci (2020) 21(16). doi: 10.3390/ijms21165662
70. Paramasivam A, George R, Priyadharsini JV. Genomic and Transcriptomic Alterations in M6a Regulatory Genes are Associated With Tumorigenesis and Poor Prognosis in Head and Neck Squamous Cell Carcinoma. Am J Cancer Res (2021) 11(7):3688–97.
71. Oh ET, Park HJ. Implications of NQO1 in Cancer Therapy. BMB Rep (2015) 48(11):609–17. doi: 10.5483/bmbrep.2015.48.11.190
72. Li X, Liu Z, Zhang A, Han C, Shen A, Jiang L, et al. NQO1 Targeting Prodrug Triggers Innate Sensing to Overcome Checkpoint Blockade Resistance. Nat Commun (2019) 10(1):3251. doi: 10.1038/s41467-019-11238-1
73. Fu L, Zhang L. Physiological Functions of CKIP-1: From Molecular Mechanisms to Therapy Implications. Ageing Res Rev (2019) 53:100908. doi: 10.1016/j.arr.2019.05.002
74. Silvers AL, Lin L, Bass AJ, Chen G, Wang Z, Thomas DG, et al. Decreased Selenium-Binding Protein 1 in Esophageal Adenocarcinoma Results From Posttranscriptional and Epigenetic Regulation and Affects Chemosensitivity. Clin Cancer Res (2010) 16(7):2009–21. doi: 10.1158/1078-0432.CCR-09-2801
75. Huang KC, Park DC, Ng SK, Lee JY, Ni X, Ng WC, et al. Selenium Binding Protein 1 in Ovarian Cancer. Int J Cancer (2006) 118(10):2433–40. doi: 10.1002/ijc.21671
76. Zeng H, Zhao X, Tang C. Downregulation of SELENBP1 Enhances Oral Squamous Cell Carcinoma Chemoresistance Through KEAP1-NRF2 Signaling. Cancer Chemother Pharmacol (2021) 88(2):223–33. doi: 10.1007/s00280-021-04284-4
77. Meneses-Morales I, Izquierdo-Torres E, Flores-Peredo L, Rodriguez G, Hernandez-Oliveras A, Zarain-Herzberg A. Epigenetic Regulation of the Human ATP2A3 Gene Promoter in Gastric and Colon Cancer Cell Lines. Mol Carcinog (2019) 58(6):887–97. doi: 10.1002/mc.22978
78. Endo Y, Uzawa K, Mochida Y, Shiiba M, Bukawa H, Yokoe H, et al. Sarcoendoplasmic Reticulum Ca(2+) Atpase Type 2 Downregulated in Human Oral Squamous Cell Carcinoma. Int J Cancer (2004) 110(2):225–31. doi: 10.1002/ijc.20118
79. Kanda M, Nomoto S, Oya H, Takami H, Hibino S, Hishida M, et al. Downregulation of DENND2D by Promoter Hypermethylation Is Associated With Early Recurrence of Hepatocellular Carcinoma. Int J Oncol (2014) 44(1):44–52. doi: 10.3892/ijo.2013.2165
80. Kanda M, Shimizu D, Nomoto S, Takami H, Hibino S, Oya H, et al. Prognostic Impact of Expression and Methylation Status of DENN/MADD Domain-Containing Protein 2D in Gastric Cancer. Gastric Cancer (2015) 18(2):288–96. doi: 10.1007/s10120-014-0372-0
81. Ling B, Zheng H, Fu G, Yuan J, Shi T, Chen S, et al. Suppression of Non-Small Cell Lung Cancer Proliferation and Tumorigenicity by DENND2D. Lung Cancer (2013) 79(2):104–10. doi: 10.1016/j.lungcan.2012.10.012
82. Zheng X, Jia B, Lin X, Han J, Qiu X, Chu H, et al. FRMD4A: A Potential Therapeutic Target for the Treatment of Tongue Squamous Cell Carcinoma. Int J Mol Med (2016) 38(5):1443–9. doi: 10.3892/ijmm.2016.2745
83. Goldie SJ, Mulder KW, Tan DW, Lyons SK, Sims AH, Watt FM. FRMD4A Upregulation in Human Squamous Cell Carcinoma Promotes Tumor Growth and Metastasis and Is Associated With Poor Prognosis. Cancer Res (2012) 72(13):3424–36. doi: 10.1158/0008-5472.CAN-12-0423
84. Tzeng HT, Wang YC. Rab-Mediated Vesicle Trafficking in Cancer. J BioMed Sci (2016) 23(1):70. doi: 10.1186/s12929-016-0287-7
85. Yoon SO, Shin S, Mercurio AM. Hypoxia Stimulates Carcinoma Invasion by Stabilizing Microtubules and Promoting the Rab11 Trafficking of the Alpha6beta4 Integrin. Cancer Res (2005) 65(7):2761–9. doi: 10.1158/0008-5472.CAN-04-4122
86. Zhang W, Chen T, Liu J, Yu S, Liu L, Zheng M, et al. RAB11FIP1: An Indicator for Tumor Immune Microenvironment and Prognosis of Lung Adenocarcinoma From a Comprehensive Analysis of Bioinformatics. Front Genet (2021) 12:757169. doi: 10.3389/fgene.2021.757169
87. Giogha C, Lung TW, Pearson JS, Hartland EL. Inhibition of Death Receptor Signaling by Bacterial Gut Pathogens. Cytokine Growth Factor Rev (2014) 25(2):235–43. doi: 10.1016/j.cytogfr.2013.12.012
88. Humphreys L, Espona-Fiedler M, Longley DB. FLIP as a Therapeutic Target in Cancer. FEBS J (2018) 285(22):4104–23. doi: 10.1111/febs.14523
89. Duan X, Cheng X, Yin X, Ke Z, Song J. Systematic Analysis of the Function and Prognostic Value of RNA Binding Protein in Head and Neck Squamous Cell Carcinoma. Eur Arch Otorhinolaryngol (2021) 279(3):1535–47. doi: 10.1007/s00405-021-06929-9
Keywords: head and neck squamous cell carcinoma, RNA binding protein, prognostic, immune microenvironment, tumor mutation burden, copy number variations, immunotherapy, chemotherapeutic
Citation: Ming R, Li X, Wang E, Wei J, Liu B, Zhou P, Yu W, Zong S and Xiao H (2022) The Prognostic Signature of Head and Neck Squamous Cell Carcinoma Constructed by Immune-Related RNA-Binding Proteins. Front. Oncol. 12:795781. doi: 10.3389/fonc.2022.795781
Received: 15 October 2021; Accepted: 04 March 2022;
Published: 05 April 2022.
Edited by:
Heming Lu, People’s Hospital of Guangxi Zhuang Autonomous Region, ChinaReviewed by:
Jessica Dal Col, University of Salerno, ItalyJohn Morton, University of Colorado, United States
Copyright © 2022 Ming, Li, Wang, Wei, Liu, Zhou, Yu, Zong and Xiao. 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: Hongjun Xiao, xhjent_whxh@hust.edu.cn; Shimin Zong, 2018XH0090@hust.edu.cn
†These authors share first authorship