Yizhou Yang1,2,3
Zeman Cai1Kaichun Huang1Mei Li1
Xiao Wang1,2,3
Yinbing Lin1,2,3Sijie Chen1,2,3Zhining Yang1,2
Zhixiong Lin1,2,3*- 1Department of Radiation Oncology, Cancer Hospital of Shantou University Medical College, Shantou, Guangdong, China
- 2Nasopharyngeal Carcinoma Research Center of Shantou University Medical College, Shantou, Guangdong, China
- 3Shantou University Medical College, Shantou, Guangdong, China
Background: Immunotherapy has been demonstrated favorable in head and neck squamous cell carcinoma (HNSCC). Studies indicated that immune-related gene prognostic index (IRGPI) was a robust signature, and N6-methyladenosine (m6A) methylation had a significant impact on the tumor immune microenvironment (TIME) and immunotherapy of head and neck squamous cell carcinoma. Thus, combining indicated that immune-related gene prognostic index with m6A status should offer a better predictive power for immune responses.
Methods: Head and neck squamous cell carcinoma samples from the cancer genome atlas (TCGA, n = 498) and gene expression omnibus database (GSE65858, n = 270) were used in this study. Cox regression analysis was used to construct the indicated that immune-related gene prognostic index through immune-related hub genes which were identified by weighted gene co-expression network analysis (WGCNA). The m6A risk score was constructed by least absolute shrinkage and selection operator (LASSO) regression analysis. Principal component analysis was used to construct a composite score, and systematically correlate subgroups according to tumor immune microenvironment cell-infiltrating characteristics.
Results: A composite score was determined based on indicated that immune-related gene prognostic index and m6A risk score. Head and neck squamous cell carcinoma patients in the cancer genome atlas were divided into four subgroups: A (IRGPI-High&m6A-risk-High, n = 127), B (IRGPI-High&m6A-risk-Low, n = 99), C (IRGPI-Low&m6A-risk-High, n = 99), and D (IRGPI-Low&m6A-risk-Low, n = 128), and overall survival (OS) was significantly different between subgroups (p < 0.001). The characteristics of tumor immune microenvironment cell infiltration in the four subgroups were significantly different in subgroups (p < 0.05). The receiver operating characteristic (ROC) curves show the predictive value of composite score for overall survival was superior to other scores.
Conclusion: The composite score is a promising prognostic signature which might distinguish immune and molecular characteristics, predict prognosis, and guide more effective immunotherapeutic strategies for head and neck squamous cell carcinoma.
Introduction
Treatment of head and neck squamous cell carcinoma (HNSCC) remains challenging (Mody et al., 2021). Immunotherapy, including immune checkpoint inhibitor therapy (anti-PD-1/L1, anti-CTLA-4, etc.), chimeric antigen receptor macrophages (CAR-M), and anti-tumor-associated macrophages (anti-TAMs), has shown promise in treatment of HNSCC, but only a small percentage of HNSCC patients have sustained immune responses (Evrard et al., 2019; Malfitano et al., 2020; Cheng et al., 2021; Sloas et al., 2021; Fasano et al., 2022; Poulose and Kainickal, 2022). Many studies have shown that immune-related gene prognostic index (IRGPI) may be a potential signature of immunotherapy for a variety of tumors. However, the existing IRGPI is not detailed enough for tumor immune microenvironment (TIME) stratification, and the accuracy of prognosis prediction needs to be further studied (Chen et al., 2021; Yao et al., 2022). Methylation to form N6-methyladenosine (m6A), as the most common modification in eukaryotic RNA, plays an either positive or negative important role on RNA synthesis, transport, and translation, indirectly affects IRGPI, TIME and immunotherapy of HNSCC (Yi et al., 2020; Jing et al., 2021). However, incomplete evaluation of m6A regulators and inadequate detailed risk stratification remain issues to be solved (Yi et al., 2020; Chen et al., 2021; Elmusrati et al., 2021; Mody et al., 2021). Consequently, we need to identify more accurately prognostic signatures associated with therapeutic benefits to guide individualized immunotherapy for HNSCC patients.
In this study, we develop a composite HNSCC prognostic signature based on IRGPI and m6A risk score, and explore the molecular and immunological characteristics of diverse subgroups of IRGPI, m6A risk score and composite score. They then were compared with microsatellite instability (MSI), tumor immune dysfunction and exclusion (TIDE), tumor inflammation signature (TIS), and other signatures.
Materials and methods
Data acquisition
Samples (n = 498, 455 tumor samples and 43 para-cancer samples) of HNSCC patients (n = 454) with RNA-seq data (FPKM value), miRNA-seq data, copy number variation (CNV) and well-defined clinical stages were screened from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov). We downloaded simple nucleotide mutations in tumor samples from 451 patients. RNA-seq data and survival information of 270 HNSCC samples (GSE65858) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo). Lists of immune-related genes were downloaded from the InnateDB (https://www.innatedb.ca/) and ImmPort (https://www.immport.org/shared/home) databases.
Identification of differentially-expressed genes (DEGs)
The RNA-seq of 498 samples (454 tumors vs. 43 para-cancer samples) obtained from TCGA was analyzed for genetic differences by the Wilcoxon test, and miRNA-seq data was analyzed using the edgeR package of R. Those genes with |log2FC| > 1 and false discovery rate (FDR) < 0.05 were considered as DEGs. After extracting the DEGs of immune-related genes obtained from ImmPort and InnateDB databases, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on immune-related DEGs by using the clusterProfiler package of R.
Identification of immune-related hub genes
WGCNA package of R was used to perform weighted gene co-expression network analysis (WGCNA) for immune-related DEGs. In order to identify hub genes, 498 HNSCC samples were clustered to detect outliers based on the immune-related DEG expression matrix, and the similarity matrix was constructed after excluding outlier data (cut-off value 36,000). Then, the similarity matrix was converted into the adjacency matrix with soft threshold β = 7 and signed network type, and topological overlap measure (TOM) was used to calculate the degree of association between genes and the adjacency matrix was converted into TOM similarity matrix. After that, the TOM similarity matrix was transformed into a dissimilarity matrix based on 1-TOM, and a dynamic pruning tree was constructed to identify the modules. Finally, we set the minimum cluster size to 25 and the merge threshold function to 0.25 to separate the six modules. The igraph package of R was used to construct a gene co-expression network for genes with edge weights greater than 0.2 in the module of interest (turquoise module). After batch adjustment for biotechnological biases, using the ComBat algorithm of the sva package in TCGA and GEO cohorts, univariate Cox regression analysis was used to identify immune-related hub genes that were associated with overall survival (OS).
Construction and validation of the IRGPI
Genes with significant impact on OS were selected from immune-related hub genes by multivariate Cox regression analysis. The IRGPI of each sample was calculated by multiplying the expression values of specific genes by their weighting in the Cox model, and then summed. Survival analysis was used to evaluate the predicting prognostic ability of IRGPI in TCGA and GEO cohorts. Patients were divided into IRGPI-high and IRGPI-low groups based on the median.
m6A RNA methylation regulators
According to published studies, we collected 29 regulators of m6A RNA methylation, including 16 readers “YTHDC1, YTHDC2, IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF2, YTHDF3, HNRNPA2B1, HNRNPC, RBMX, FMR1, LRPPRC, ELAVL1, EIF3A, and PRRC2A”, 10 writers “METTL3, METTL5, METTL14, METTL16, WTAP, VIRMA, RBM15, RBM15B, ZC3H13, and CBLL1” and three erasers “FTO, ALKBH3, and ALKBH5”.
Calculation of m6A risk score
Based on the gene expression of 29 regulators of m6A RNA methylation, the survival package of R was used for univariate Cox regression analysis to construct a prognostic network map. The glmnet package of R was used for the least absolute shrinkage and selection operator (LASSO) regression analysis. After simulation 1,000 times and cross-validation, the m6A regulators that strongly correlated with prognosis and their coefficients were obtained. Then, the m6A risk score was calculated by multiplying each obtained coefficient by the corresponding m6A regulator expression and summing the total values. HNSCC patients were divided into m6A-risk-high and m6A-risk-low groups according to the median value, and the differences in OS between the two groups were assessed.
Classification and composite score
Four subgroups were classified by pairwise combination of IRGPI and m6A risk score. IRGPI and m6A risk score were standardized by using the descriptives function of SPSS, and their weights were calculated by the principal component method in factor analysis. The composite score equaled the sum of the normalized IRGPI and m6A risk score multiplied by their corresponding weights. Univariate and multivariate Cox regression analyses were performed to verify the independent prognostic value of IRGPI, m6A risk score and composite score, and their relevant gene signatures were constructed by principal component analysis (PCA), with principal components 1, 2, and 3 selected as signature scores. The chi-square test (Monte-Carlo algorithm) was applied to test the correlation between clinicopathological features and composite score subgroups. Based on the interaction between mRNA and miRNA in the ENCORI database (http://starbase.sysu.edu.cn/) (Li et al., 2014), we used Cytoscape (v3.8.2) to construct the target interaction network of differentially-expressed miRNAs and mRNAs that were expressed by the genes used to construct the three scores. To further explore the interrelationships between proteins expressed by the genes reflected in the scores, we used the STRING database (https://www.string-db.org/) to construct the protein-protein interaction (PPI) network.
Comprehensive analysis of the molecular immunological characteristics and immunotherapy in subgroups
Based on the DEGs, gene set enrichment analysis (GSEA) was applied to the composite score subgroups with the GO, KEGG and HALLMARK gene sets, which were downloaded from the GSEA database (http://www.gsea-msigdb.org/gsea/index.jsp) using the clusterProfiler package of R. The number and quality of genetic mutations in subgroup tumor samples was analyzed by using the maftools package of R. Correlation analysis of the composite score with PD-1/L1, CTLA-4 and tumor mutation burden (TMB) were performed. To identify the immune characteristics of 454 HNSCC samples, we analyzed their RNA-seq data by using the CIBERSORT tool (https://cibersort.stanford.edu/) and performed 1,000 iterations to estimate the relative proportions of 22 kinds of immune cells among subgroups. We analyzed the expression differences of 22 kinds of immune cells among subgroups, and selected differentially-expressed immune cells (p < 0.05) to analyze their influence on OS, after using the survminer package of R to calculate the optimal cutoff value of the corresponding immune cells for dividing tumor samples into two groups depending on the degree of infiltration. To further clarify immunological and molecular functions among subgroups, we performed single sample GSEA (ssGSEA) by using the GSVA package of R, and compared normalized ssGSEA scores among subgroups. Survival analysis was conducted for items with statistical differences (p < 0.05) according to the optimal cutoff value, which also calculated by the survminer package of R, of the items for dividing tumor samples into two groups. We evaluated the percentages of the three scores in immunized subtypes and final survival states. By using the time ROC package of R, we performed time-dependent receiver operating characteristic (ROC) curve analyses to obtain the area under the curve (AUC) and compared the prognostic value among composite score, m6A risk score, IRGPI, TIS, TIDE, interferon gamma (IFNG), MSI, Merck18, CD274, CD8, Dysfunction, Exclusion, myeloid-derived suppressor cells (MDSC), tumor-associated fibroblast (CAF), and TAM (M2) scores. The TIS score was calculated as the mean of normalized expression on a log2-scale of 18 characteristic genes (Ayers et al., 2017). Other scores were obtained online from the TIDE database (http://tide.dfci.harvard.edu/).
Statistical analysis
The correlation among variables involved in the study was elucidated by Pearson correlation analysis. An independent t-test was performed to compare continuous variables among subgroups. Categorical data was tested by using the chi-square test. The Wilcoxon test was used to compare differences of TIS and scores obtained from TIDE database between subgroups. Multivariate survival analysis was performed using the Cox regression model, while other survival analyses were performed by Kaplan-Meier survival analysis with the log-rank test by using the survival package of R, and p < 0.05 was considered to indicate statistical significance. All the above analysis processes were shown in Supplementary Figure S1. All data was processed in SPSS (v.25.0) and R (v.4.1.2) software.
Results
Outcomes of differential expression analysis
In differential expression analysis of 454 HNSCC patient samples (43 normal vs. 454 tumor samples), 4,855 DEGs were identified, which included 3,550 upregulated genes and 1,305 downregulated genes (Supplementary Figures S2A, B), among which 495 upregulated and 143 downregulated were immune-related DEGs (Supplementary Figures S2C, D). Furthermore, a total of 312 differentially-expressed miRNAs were observed, of which 186 were upregulated and 126 were downregulated (Supplementary Figures S2E, F). Enrichment analysis showed that the 638 immune-related DEGs were associated with 1,995 GO items and 115 KEGG pathways (q < 0.05; p < 0.05, Supplementary Tables S1, 2; Supplementary Figures S2G, H).
Immune-related hub genes and the IRGPI
To obtain immune-related hub genes, 498 samples were clustered to detect outliers (cutoff value 36,000) based on the 638 immune-related DEGs matrix, and the data was analyzed using WGCNA after deleting two outlier samples (Supplementary Figures S3A, B). The logarithm of nodes with connectivity of K (log(K)) was negatively correlated with the logarithm log[P(K)] of node probability, and the correlation coefficient was greater than 0.85. Based on a scale-free network, the optimal soft threshold power was 7 (Supplementary Figures S3C, D). The 638 immune-related DEGs were allocated to six modules, which were determined based on average linkage hierarchical clustering and optimal soft threshold power (Supplementary Figures S3E, F). According to the gene dendrogram and the Pearson correlation coefficient between module and sample feature for each module, the turquoise module, which was closely related to HNSCC and relatively concentrated genes distribution, was selected for further study (Figures 1A, B). There were 34 nodes and 115 edges of the turquoise module of the networks with a threshold weight >0.2 (Supplementary Figure S3G). Survival analysis showed that the expression of the 11 hub genes was closely related to the OS of HNSCC patients (p < 0.05, Supplementary Figures S4A–K). Turquoise module DEGs were analyzed by univariate Cox regression analysis, and 11 genes with statistical significance were selected as immune-related hub genes (p < 0.05, Figure 1C). We explored somatic mutations in the 11 hub genes and their clinicopathological information (Supplementary Table S3). Most of the 11 hub genes had missense mutations and non-sense mutations, but only the mutation rates of MAPT, EGF, PTX3, and SEMA3G were greater than 1% (Supplementary Figure S4L).
FIGURE 1. Identification of immune-related hub genes and comprehensive assessment of IRGPI. (A) Weighted gene coexpression network analysis (WGCNA) of immune-related differentially expressed genes with a soft threshold β = 7. (B) Gene modules related to HNSCC obtained by WGCNA. (C) Univariate Cox analysis of 11 immune-related hub genes. (D) K–M survival analysis of the IRGPI subgroups in TCGA cohort. (E) K–M survival analysis of the IRGPI subgroups in the GEO cohort. (F–M) The difference of PD-1, PD-L1, CTLA-4 and TMB in IRGPI subgroups (Wilcoxon test) and their correlation analysis with IRGPI.
Multivariate Cox regression analysis on the 11 hub genes identified and provided coefficients for seven independent prognostic hub genes (HTN3, CTSG, MAPT, CCL28, PTX3, SEMA3G, and USP2). A prognostic index for tumor samples could be calculated by the formula: IRGPI = (HTN3 expression * 0.096) + (CTSG expression * −0.152) + (MAPT expression * 0.122) + (CCL28 expression * −0.094) + (PTX3 expression * 0.138) + (SEMA3G expression * −0.140) + (USP2 expression * 0.107). With the median IRGPI as the threshold, patients with high IRGPI had worse OS than those with low IRGPI (p < 0.001, Figure 1D). This result was verified in the HNSCC dataset of GSE65858 (n = 270) (p < 0.05, Figure 1E).
Expression of PD-1 and CTLA-4 were high in the IRGPI-low subgroup and correlation analysis indicated slightly negatively correlated with the IRGPI (Figures 1F–I). However, TMB and expression of PD-L1 did not show significant differences and relationship with IRGPI (Figures 1J–M).
Landscape and risk score of m6A regulators in HNSCC
There were 29 m6A regulators in this study, including 10 writers, three erasers, and 16 readers. As shown in Figures 2A, B, 23 of them, including 12 readers, 8 writers, and all erasers, were differentially expressed between HNSCC tumor and para-cancer samples, and were upregulated in tumor samples except for the reader YTHDC2. Based on gene expression of regulators, we analyzed the correlation between 29 m6A regulators in pairs with a p < 0.001 signifying a correlating pair. The results showed that m6A regulators were generally positively correlated, but only the reader YTHDF3 and PRRC2A were highly correlated (r = 0.77, Supplementary Figure S5A).
FIGURE 2. Landscape of m6A regulators and m6A risk score in HNSCC. (A) Heatmap of 29 m6A regulators in TCGA cohort. (B) The upper and lower ends of the boxes represented interquartile range of values. The lines in the boxes represented median value, and dots showed outliers. The asterisks represented the statistical p-value (*p < 0.05; **p < 0.01; ***p < 0.001). (C) The CNV variation frequency of m6A regulators. The height of the column represented the alteration frequency. (D) The location of CNV alteration of m6A regulators on 23 chromosomes. (E) Univariate Cox analysis of 29 m6A regulators. (F) The interaction between m6A regulators in HNSCC. (G, H) Determination of the m6A regulators strongly correlated with prognosis by the LASSO analysis. (I) K–M survival analysis of the m6A risk score subgroups in TCGA cohort. (J–Q) The difference of PD-1, PD-L1, CTLA-4, and TMB in m6A risk score subgroups (Wilcoxon test) and their correlation analysis with m6A risk score.
We also explored mutations and clinicopathological information of the m6A regulators (Supplementary Figure S5B). In 451 HNSCC samples, m6A regulator mutations occurred in 68 cases, with a mutation frequency of 15.08%. The results showed that readers PRRC2A and YTHDC2 had the highest mutation frequency (about 2%) among the 29 regulators. The CNV frequency showed it was very common among regulators, with 15 undergoing copy number amplification, and 14 undergoing copy number deletion (Figure 2C). The chromosomal locations of m6A regulators CNV were shown in Figure 2D.
Univariate Cox regression demonstrated that IGF2BP2, YTHDF1, HNRNPA2B1, HNRNPC, RBMX, LRPPRC, METTL5, ALKBH3, and ALKBH5 were correlated with OS (Supplementary Figure S6), with hazard ratios all greater than 1 (p < 0.05, Figure 2E). Based on the results of univariate Cox regression analysis, we constructed an m6A regulator network (Figure 2F), which showed the comprehensive landscape of the correlation of regulators in pairs and their prognostic significance for HNSCC. We found that all regulators except YTHDC2 were prognostic risk factors, and there was a widespread positive correlation between the expression of regulators. The LASSO algorithm was used to obtain the coefficient of regulators with prognostic value (Figures 2G, H), IGF2BP2, HNRNPC, and METTL5 were selected to construct a prognostic signature using the cohort in TCGA. The formula was as follows: m6A risk score = (0.067 * IGF2BP2 expression) + (0.119 * HNRNPC expression) + (0.487 * METTL5 expression). According to the median risk score, HNSCC patients with a low m6A risk score had better OS than those with a high one (p < 0.001, Figure 2I).
The expression of PD-1 and CTLA-4 were high and TMB was low in the m6A-risk-low subgroup and correlation analysis indicated slightly positively correlated with m6A risk score (Figures 2J–O). Expression of PD-L1 did not show significant differences, but negative correlated with the score (Figures 2P, Q).
Classification and composite score
Based on the IRGPI and m6A risk score, HNSCC patients in TCGA were pairwise combined as subgroup A (IRGPI-High&m6A-risk-High, n = 127), subgroup B (IRGPI-High&m6A-risk-Low, n = 99), subgroup C (IRGPI-Low&m6A-risk-High, n = 99) and subgroup D (IRGPI-Low&m6A-risk-Low, n = 128). IRGPI and m6A risk score were standardized and multiplied by the weight coefficient obtained from PCA to obtain the composite score, using the formula: Composite score = 0.5*Z-score (IRGPI) + 0.5*Z-score (m6A risk score). Survival analysis showed from subgroup A to D, the OS was better than that of the previous subgroup (p < 0.001, Figure 3A). The PCA based on subgroup gene expressions is shown in Figures 3B–D. The correlation of the composite score and clinicopathological information in Figure 3E showed the correlation between the composite score and the primary site was statistically significant (p < 0.001). Univariate Cox regression analysis showed that age, stage, radiotherapy, IRGPI, m6A risk score and composite score were significantly correlated with the prognosis of HNSCC (Figure 3F). Multivariate Cox regression analysis confirmed that the IRGPI and m6A risk score were independent prognostic factors after adjusting for other clinicopathological factors (Figure 3G), and the composite score was also an independent prognostic factor after excluding the effects of IRGPI and m6A risk score (Figure 3H). Combining the 10 genes, used to construct composite score, with differentially-expressed miRNA to construct a gene targeted network with 24 nodes and 18 edges showed that miRNA was mainly upregulated to inhibit gene expression in the network (Figure 3I). The PPI network revealed gene fusion between HNRNPC and MAPT, adjacent of MAPT and SEMA3G, and other interrelationships obtained from text mining (Figure 3J; Supplementary Table S4).
FIGURE 3. Assessment of the composite score. (A) K–M survival analysis of the composite score subgroups in TCGA cohort. (B–D) Principal component analysis of the transcriptome profiles of the three score-defined subgroups. (E) Heatmap of clinical relevance for composite score subgroups (***p < 0.001). (F) Univariate Cox analysis of clinicopathologic factors, IRGPI, m6A risk score and the composite score. (G) Multivariate Cox analysis of the factors significant in the univariate Cox analysis (p < 0.05). (H) Multivariate Cox analysis of composite score and clinicopathologic factors significant in the univariate Cox analysis (p < 0.05). (I) The target interaction network of differentially expressed miRNA and the 10 mRNAs which are expressed by the genes used to construct the three scores. (J) The Protein-Protein Interaction (PPI) Network of 10 genes used to construct the three scores.
Molecular characteristics of different subgroups
GSEA was performed to determine the GO, KEGG and Hallmark gene sets enriched in the different subgroups (p < 0.05, Supplementary Table S5). The gene sets of subgroup A were markedly enriched in tumor metastasis and drug resistance (Supplementary Figure S7A), the gene sets of subgroup B were markedly enriched in tumor growth promotion and immune rejection (Supplementary Figure S7B), the gene sets of subgroup C were markedly enriched in immune rejection and DNA repair (Supplementary Figure S7C), and the gene sets of subgroup D were markedly enriched in immune response (Supplementary Figure S7D).
The top 20 highest mutation rate genes in TCGA cohort demonstrated that the most common mutation type was missense mutation, followed by non-sense mutation and frameshift deletion, as plotted in Figures 4A–D. The mutation frequencies of TP53, TTN, and FAT1 were all greater than 15% among the four subgroups, CDKN2A, MUC16, CSMD3, and PIK3CA were greater than 15% in subgroups A, C, and D, KMT2D was greater than 15% in subgroups A and C, SYNE1 was greater than 15% in subgroups C and D, PCLO was greater than 15% in subgroup A, NOTCH1 was greater than 15% in subgroup B, and LRP1B, NSD1, and DNAH5 were greater than 15% in subgroup C.
FIGURE 4. Characteristics of tumor somatic mutations, PD-1, PD-L1, CTLA-4, and TMB in composite subgroups. (A–D) Mutation profiles of the 20 genes with the highest mutation rates in HNSCC samples of the TCGA cohort in the four composite score subgroups. (E–H) The correlation analysis of PD-1, PD-L1, CTLA-4, and TMB with the composite score. (I–L) The difference of PD-1, PD-L1, CTLA-4, and TMB in composite score subgroups (Wilcoxon test).
We discovered that the composite score was moderately negative correlated with PD-1, slightly negative correlated with CTLA-4 and PD-L1, and slightly positive correlated with TMB (Figures 4E–H), and their differences across subgroups were shown in Figures 4I–L.
Immune characteristics of different subgroups
The abundance of 22 types of immune cells in different subgroups was analyzed with the Wilcoxon test to compare their distribution (Supplementary Table S6). The abundance of native B cells, plasma cells, CD8+ T cells, activated memory CD4+ T cells, follicular helper T cells, regulatory T cells (Tregs), resting natural killer (NK) cells, M0 macrophages, resting dendritic cells, resting mast cells, eosinophils and neutrophils were significantly different in subgroups (p < 0.05, Figures 5A,B). The abundance of M0 macrophages was the highest in all subgroups, and the patients with higher abundance of M0 macrophages had lower OS. In contrast, patients with the high abundance of native B-cell, plasma cells, CD8+ T cells, follicular helper T cells and Tregs had better OS. Based on the optimal cutoff value of the abundance of each immune cell, HNSCC patients were divided into two groups for survival analysis, which indicated that the more abundance of follicular helper T cells, Tregs, activated memory CD4+ T cells, CD8+ T cells, activated B cells, plasma cells, resting dendritic cells and resting mast cells, the better the OS. In contrast, a high abundance of neutrophils, activated mast cells, and M0 and M2 macrophages tended to imply poor OS (Supplementary Figure S8).
FIGURE 5. The TIME landscape and the immune subtype of different composite subgroups. (A) The proportions of TIME cells in composite score subgroups for patients in TCGA cohort. (B) The abundance of each TIME infiltrating cell in composite score subgroups. The upper and lower ends of the boxes represented interquartile range of values. The lines in the boxes represented median value, and dots showed outliers. The asterisks represented the statistical p-value (*p < 0.05; **p < 0.01; ***p < 0.001). (C) The molecular and immune-related function of composite score subgroups. The gene sets of molecular and immune-related function were analyzed by the single simple gene set enrichment analysis (ssGSEA) and then compared between different composite score subgroups. The upper and lower ends of the boxes represented interquartile range of values. The lines in the boxes represented median value, and dots showed outliers. The asterisks represented the statistical p-value (*p < 0.05; **p < 0.01; ***p < 0.001). (D) The Alluvial diagram shows changes in IRGPI, m6A risk score, composite score, immune subtype, and survival status of patients in the TCGA cohort.
Next, we applied specific gene signatures by ssGSEA to define immunity and molecular function among the subgroups (Supplementary Table S7). As the results illustrate, the better the OS of the subgroup, the more cytolytic activity, mast cells and NK cells (Figure 5C). Also, a variety of T cells and plasma cells, immune suppression and signals related to tumor metastasis had statistical significance, which were more obvious in subgroups A and D, and these immune signatures occurred most commonly in subgroup D with the best OS. Based on the optimal cutoff value of each immune and molecular function signature for the tumor samples in TCGA, it revealed the signatures had a large impact on survival. Patients with more macrophages and parainflammatory cells had poor prognoses, while patients with more T cells, mast cells, NK cells, B cells and other signatures has better prognosis (Supplementary Figure S9).
The benefit of immunotherapy in different subgroups
Thorsson et al. (2018) divided all TCGA tumor samples into six immune subtypes, which included wound healing (Immune C1), IFN-gamma dominant (Immune C2), inflammatory (Immune C3), lymphocyte-depleted (Immune C4), immunologically quiet (Immune C5) and TGF-beta-dominant (Immune C6), and showed patients who belonged to types C1, C2 and C3 could benefit from neoadjuvant immunotherapy following surgery. We explored the proportion of HNSCC samples in TCGA cohort and found that almost all samples belonged to the C1 and C2 immune subtypes. In each subgroup, C2 was predominant with similar proportion of C1 and C2 (Figure 5D; Supplementary Table S8).
We used TIDE, IFNG, MSI, Merck18, CD274, CD8, Dysfunction, Exclusion, MDSC, CAF, and TAM (M2) scores to evaluate the potential clinical efficacy of immunotherapy in subgroups categorized by IRGPI, m6A risk score and composite score. In the IRGPI-high subgroup, TAM (M2), Exclusion, CAF and MDSC scores were higher, while CD8 and IFNG scores were higher in the IRGPI-low subgroup (Supplementary Figure S10A). TAM (M2), Exclusion and MDSC scores were higher in the m6A risk-high subgroup, while TIDE, CD8, Dysfunction, IFNG, and Merck18 scores were higher in the m6A risk-low subgroup (Supplementary Figure S10B). There were no statistically significant differences in the IRGPI and m6A risk score subgroups for the remaining scores. In the subgroups of composite score, the scores of CD274 and MSI did not differ among subgroups, while TIDE, Dysfunction, Exclusion, TAM (M2), IFNG, CD8, CAF, MDSC, and Merck18 scores demonstrated significant differences between subgroup A and subgroup D, but their differences among other subgroups were not uniform (Figure 6A). In general, patients in subgroup A may have the best benefit from immunotherapy, while patients in subgroups B and C had better immunotherapy benefit than patients in subgroup D. However, the magnitude of immunotherapy benefit between subgroups B and C remains to be further delineated.
FIGURE 6. Prognostic value of the composite score in the TCGA cohort. (A) TIDE, Dysfunction, Exclusion, TAM (M2), IFNG, CD8, CAF, MDSC, and Merck18 scores in different composite score subgroups. The scores between the two subgroups were compared through the Wilcoxon test (ns, not significant, *p < 0.05; **p < 0.0s1; ***p < 0.001). (B) The predictive value of composite score in TCGA HNSCC cohort at 1, 2, 3, and 5 years. (C) ROC analysis of composite score, IRGPI, m6A risk score and the scores of TIS, TIDE, IFNG, MSI, Merck18, CD274, and CD8, Dysfunction, Exclusion, MDSC, CAF, and TAM (M2) on OS at 1-, 2-, three- and 5-year follow-up in TCGA HNSCC cohort.
Subsequently, we evaluated the OS of the composite score and other prognostic indicators at 1, 2, 3, and 5 years by ROC curves, and found that prediction of OS within 3 years by the composite score was better than the long-term survival rate (Figure 6B, AUC = 0.663 (1-year), AUC = 0.679 (2-year), AUC = 0.687 (3-year), AUC = 0.640 (5-year)). Moreover, the predictive value of composite score for OS was superior to any other scores (Figure 6C).
Discussion
Increasing evidence indicates that multiple m6A regulators perform a vital role in regulating the occurrence, development, metastasis, drug resistance, radiotherapy resistance and other aspects of HNSCC (Liu et al., 2020; Jing et al., 2021; Wu et al., 2021; Jin et al., 2022; Yu et al., 2022). At the same time, in-depth research has uncovered more and more m6A regulators. Currently, characteristics of the HNSCC TIME mediated by the combined action of multiple m6A regulators have not been sufficiently explained (Jing et al., 2021). On the other hand, immunotherapy is widely applied as an effective method to treat HNSCC, but the overall response rate remains low. After years of exploring HNSCC prognostic signatures, while some simple signatures can forecast immunotherapy and OS (Fasano et al., 2022; Poulose and Kainickal, 2022), single prognostic signatures cannot very well predict which patients can benefit from immunotherapy. Compound signatures constructed by combining multiple signatures can further differentiate the prognosis and immune response of patients with diverse OS, which will contribute to our understanding of anti-tumor immune response of TIME and guide more effective immunotherapy strategies. It is inevitable that future research will focus on identifying compound immunotherapy prognosis signatures.
Here, we used different methods to establish three different prognostic scores, all of which are valid independent prognostic signatures for HNSCC. Compared with IRGPI and m6A risk score, the composite score is more detailed for OS stratification of HNSCC patients, and the AUC value is greater than other scores, which means that the prediction of prognosis is also better.
In this research, based on 29 m6A regulators, we reveal the m6A methylation landscape of HNSCC in more detail. We revealed that mutation rate of the 29 m6A regulators was not significant. Interestingly, about half of them showed increasing of CNV which might indicated that m6A regulators did not affect HNSCC primarily through mutation and CNV. Except for YTHDC2, all regulators were highly expressed in HNSCC samples and were prognostic risk factors. However, publications on m6A regulators in HNSCC are limited, and further studies are warranted to clarify the effect of m6A on HNSCC.
On the basis of the composite score, we further studied gene mutations of subgroups. The greatest differences in gene mutations among subgroups were TP53 mutations, which were more frequent in subgroup A than in subgroup D (72% vs. 46%). TP53 mutation is not only the most common genetic event in cancer, but is also associated with more aggressive disease and poorer patient outcomes in many cancers, especially in HNSCC (Poeta et al., 2007; Leroy et al., 2014). The Notch signaling pathway is known to be an evolutionarily conserved intercellular signaling pathway, and NOTCH1 mutation is associated with HNSCC (Schmidl et al., 2022; Zhou et al., 2022). The mutation rate of NOTCH1 gene in subgroup B (23%) is higher than that in other groups, which may be the reason why HNSCC patients in subgroup B have worse OS than those in subgroup C. It also means that some HNSCC patients may benefit from anti-NOTCH1 therapy. The TTN mutation rate of subgroup C (48%) is higher than that of other groups, and is mainly comprised of multiple mutations, but its prognosis and biological significance of HNSCC is unknown. Subsequently, we explored the relationship between three scores and known immunotherapy predictive signatures PD-1, PD-L1, CTLA-4, and TMB, which was evaluated in prospective clinical trials as a potential signature for predicting immune checkpoint inhibitor therapy response in HNSCC (Cristescu et al., 2018; Chan et al., 2019). In addition, a number of clinical trials have made breakthroughs in immunotherapy using immune checkpoint inhibitors against HNSCC, and a number of combined anti-CTLA-4 and anti-PD-1/L1 clinical trials are also underway (Ferris et al., 2020; Maio et al., 2021; Masarwy et al., 2021). The expression of PD-1 and CTLA-4 was synochronous in each composite score subgroups, which might indicate that combination of anti-CTLA-4 and anti-PD-1 therapy would increase the anti-cancer effects, especially in subgroup D.
Understanding the TIME landscape of HNSCC can help discover new treatments or alter the TIME to improve the effectiveness of immunotherapy (Feng and Hess, 2021; Watermann et al., 2021). The immunotyping by Thorsson et al. (2018) revealed that almost all HNSCC patients should receive neoadjuvant immunotherapy. As we know, immune-excluded and immune-desert phenotypes can be considered non-inflammatory tumors (Bagaev et al., 2021). Subgroup A corresponding to the immune-desert phenotype is characterized by immunosuppression, subgroups B and C corresponding to the immune-excluded phenotype are characterized by the activation of innate immunity and stroma, and subgroup D corresponding to the immune-inflamed phenotype is characterized by activation of adaptive immunity, known as a hot tumor (31), which is characterized by a large infiltration of immune cells. After exploring the characteristics of immune cell infiltration in the TIME of different subgroups, although the degree of immune infiltration in subgroup B was similar to that in subgroup C, low NK and mast cell infiltration and low cytolytic activity might be the reason for the worse prognosis of subgroup B compared to subgroup C (Schantz et al., 1986; Attramadal et al., 2016; Gu et al., 2019; Charap et al., 2020; Watermann et al., 2021). However, the tumor immune cell infiltration in subgroup A was the least, indicating a possible worse result of immunotherapy comparing to any other subgroup. In the TIME of HNSCC, infiltration of undifferentiated macrophages (M0) was the highest among the 22 types of immune cells, and there was a significant difference of the composite score between the subgroups. TAM (M2) is associated with elevated tumor growth, aggressive phenotype and poor prognosis in solid tumors (Komohara et al., 2016; Ngambenjawong et al., 2017; Chen et al., 2019). Although there was no statistical difference in the expression of TAM (M2) between the composite score subgroups, the infiltration by TAM (M2) was high and statistically significant in the high-risk IRGPI or m6A risk score subgroups, which is consistent with the results of our survival analysis. In addition, tumor drug resistance is the main cause of chemotherapy failure, and studies have shown that TAMs are closely related to the drug resistance of tumors (Lopez-Yrigoyen et al., 2021; Pu and Ji, 2022). CAR-M have been demonstrated to enter solid tumors and survive in the tumor environment, triggering long-lasting adaptive immune responses (Komohara et al., 2016; Sloas et al., 2021). Accordingly, the excessive TAM infiltration indicates that HNSCC patients are expected to benefit from anti-TAM and CAR-M therapy, while the patients of subgroup A corresponding to the immune desert type may benefit the most, which also means that HNSCC patients who fail other immunotherapies or chemotherapy may benefit. In addition, MDSCs are the precursors of dendritic cells, macrophages and granulocytes, and have the ability to inhibit the immune response (Draghiciu et al., 2015; Greene et al., 2020; Grover et al., 2021). The results of our study indicate that subgroup A patients, who have worse prognosis, also have a higher MDSC score. By evaluating the MDSC score of subgroups, defined by IRGPI, m6A risk score and composite score, we found that the composite score was more favorable in identifying the MDSC-related risks of HNSCC. This helps identify HNSCC patients who may benefit from anti-MDSC-targeted therapy. Evidently, compared with IRGPI and m6A risk scores, the composite score can provide a more precise grouping of HNSCC patients, which is conducive to a more detailed understanding of the TIME of HNSCC and guide more effective immunotherapy strategies.
Finally, we compared the predictive prognostic ability of composite score with IRGPI, m6A risk score, TIS, TIDE, IFNG, MSI, Merck18, CD274, CD8, Dysfunction, Exclusion, MDSC, CAF, and TAM (M2) scores on the HNSCC cohort in TCGA. The results showed that the composite score was better than other scores in predicting both short-term and long-term survival. Although the AUC value of composite score declined slightly at 5 years, this might due to the limited number of patients with a 5-year OS.
In conclusion, the composite score, constructed based on IRGPI and m6A risk score, is an independent prognostic signature, which is more useful in identifying immune and molecular characteristics, predicting patient prognosis and guiding therapy as compared with single prognostic signatures. Further studies are still needed to clarify that the composite score is potential prognostic indicator for immunotherapy.
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 author.
Author contributions
All authors contributed to the study conception and design. ZY and ZL conceived this study as well as acquisition and analysis obtain data. ZY collect clinical information and writing the manuscript. ZC, XW, SC, YL, and KH prepared figure and provided comments. ZY, ML, KH, and ZL were responsible for the statistical analysis of the data. ZL and ZY were responsible for revisions to the manuscript. The entire team helped contribute to and accepted the final manuscript.
Funding
Science and Technology Special Fund of Guangdong Province of China (190829105556145).
Acknowledgments
YY, ZL, and ZC. A composite score based on immune-related gene prognostic index and m6A risk score of head and neck squamous cell carcinoma, 26 August 2022, PREPRINT (Version 1) available at Research Square (https://doi.org/10.21203/rs.3.rs-1987322/v1).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1061569/full#supplementary-material
References
Attramadal, C. G., Kumar, S., Gao, J., Boysen, M. E., Halstensen, T. S., and Bryne, M. (2016). Low mast cell density predicts poor prognosis in oral squamous cell carcinoma and reduces survival in head and neck squamous cell carcinoma. Anticancer Res. 36 (10), 5499–5506. Epub 2016/11/01. doi:10.21873/anticanres.11131
Ayers, M., Lunceford, J., Nebozhyn, M., Murphy, E., Loboda, A., Kaufman, D. R., et al. (2017). Ifn-gamma-related mrna profile predicts clinical response to Pd-1 blockade. J. Clin. Invest. 127 (8), 2930–2940. Epub 2017/06/27. doi:10.1172/JCI91190
Bagaev, A., Kotlov, N., Nomie, K., Svekolkin, V., Gafurov, A., Isaeva, O., et al. (2021). Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell 39 (6), 845–865.e7. Epub 2021/05/22. doi:10.1016/j.ccell.2021.04.014
Chan, T. A., Yarchoan, M., Jaffee, E., Swanton, C., Quezada, S. A., Stenzinger, A., et al. (2019). Development of tumor mutation burden as an immunotherapy biomarker: Utility for the oncology clinic. Ann. Oncol. 30 (1), 44–56. Epub 2018/11/06. doi:10.1093/annonc/mdy495
Charap, A. J., Enokida, T., Brody, R., Sfakianos, J., Miles, B., Bhardwaj, N., et al. (2020). Landscape of natural killer cell activity in head and neck squamous cell carcinoma. J. Immunother. Cancer 8 (2), e001523. Epub 2021/01/12. doi:10.1136/jitc-2020-001523
Chen, Y., Li, Z. Y., Zhou, G. Q., and Sun, Y. (2021). An immune-related gene prognostic index for head and neck squamous cell carcinoma. Clin. Cancer Res. 27 (1), 330–341. Epub 2020/10/25. doi:10.1158/1078-0432.CCR-20-2166
Chen, Y., Song, Y., Du, W., Gong, L., Chang, H., and Zou, Z. (2019). Tumor-associated macrophages: An accomplice in solid tumor progression. J. Biomed. Sci. 26 (1), 78. Epub 2019/10/21. doi:10.1186/s12929-019-0568-z
Cheng, N., Bai, X., Shu, Y., Ahmad, O., and Shen, P. (2021). Targeting tumor-associated macrophages as an antitumor strategy. Biochem. Pharmacol. 183, 114354. Epub 2020/12/07. doi:10.1016/j.bcp.2020.114354
Cristescu, R., Mogg, R., Ayers, M., Albright, A., Murphy, E., Yearley, J., et al. (2018). Pan-tumor genomic biomarkers for Pd-1 checkpoint blockade-based immunotherapy. Science 362 (6411), eaar3593. Epub 2018/10/13. doi:10.1126/science.aar3593
Draghiciu, O., Lubbers, J., Nijman, H. W., and Daemen, T. (2015). Myeloid derived suppressor cells-an overview of combat strategies to increase immunotherapy efficacy. Oncoimmunology 4 (1), e954829. Epub 2015/05/08. doi:10.4161/21624011.2014.954829
Elmusrati, A., Wang, J., and Wang, C. Y. (2021). Tumor microenvironment and immune evasion in head and neck squamous cell carcinoma. Int. J. Oral Sci. 13 (1), 24. Epub 2021/08/04. doi:10.1038/s41368-021-00131-7
Evrard, D., Szturz, P., Tijeras-Raballand, A., Astorgues-Xerri, L., Abitbol, C., Paradis, V., et al. (2019). Macrophages in the microenvironment of head and neck cancer: Potential targets for cancer therapy. Oral Oncol. 88, 29–38. Epub 2019/01/09. doi:10.1016/j.oraloncology.2018.10.040
Fasano, M., Corte, C. M. D., Liello, R. D., Viscardi, G., Sparano, F., Iacovino, M. L., et al. (2022). Immunotherapy for head and neck cancer: Present and future. Crit. Rev. Oncol. Hematol. 174, 103679. Epub 2022/04/09. doi:10.1016/j.critrevonc.2022.103679
Feng, B., and Hess, J. (2021). Immune-related mutational landscape and gene signatures: Prognostic value and therapeutic impact for head and neck cancer. Cancers (Basel) 13 (5), 1162. Epub 2021/04/04. doi:10.3390/cancers13051162
Ferris, R. L., Haddad, R., Even, C., Tahara, M., Dvorkin, M., Ciuleanu, T. E., et al. (2020). Durvalumab with or without tremelimumab in patients with recurrent or metastatic head and neck squamous cell carcinoma: Eagle, a randomized, open-label phase iii study. Ann. Oncol. 31 (7), 942–950. Epub 2020/04/16. doi:10.1016/j.annonc.2020.04.001
Greene, S., Robbins, Y., Mydlarz, W. K., Huynh, A. P., Schmitt, N. C., Friedman, J., et al. (2020). Inhibition of mdsc trafficking with sx-682, a cxcr1/2 inhibitor, enhances nk-cell immunotherapy in head and neck cancer models. Clin. Cancer Res. 26 (6), 1420–1431. Epub 2019/12/19. doi:10.1158/1078-0432.CCR-19-2625
Grover, A., Sanseviero, E., Timosenko, E., and Gabrilovich, D. I. (2021). Myeloid-derived suppressor cells: A propitious road to clinic. Cancer Discov. 11 (11), 2693–2706. Epub 2021/10/13. doi:10.1158/2159-8290.CD-21-0764
Gu, X., Boldrup, L., Coates, P. J., Fahraeus, R., Wang, L., Wilms, T., et al. (2019). High immune cytolytic activity in tumor-free tongue tissue confers better prognosis in patients with squamous cell carcinoma of the oral tongue. J. Pathol. Clin. Res. 5 (4), 240–247. Epub 2019/06/27. doi:10.1002/cjp2.138
Jin, S., Li, M., Chang, H., Wang, R., Zhang, Z., Zhang, J., et al. (2022). The m6A demethylase ALKBH5 promotes tumor progression by inhibiting RIG-I expression and interferon alpha production through the IKKε/TBK1/IRF3 pathway in head and neck squamous cell carcinoma. Mol. Cancer 21 (1), 97. Epub 2022/04/10. doi:10.1186/s12943-022-01572-2
Jing, F. Y., Zhou, L. M., Ning, Y. J., Wang, X. J., and Zhu, Y. M. (2021). The biological function, mechanism, and clinical significance of M6a rna modifications in head and neck carcinoma: A systematic review. Front. Cell Dev. Biol. 9, 683254. Epub 2021/06/18. doi:10.3389/fcell.2021.683254
Komohara, Y., Fujiwara, Y., Ohnishi, K., and Takeya, M. (2016). Tumor-associated macrophages: Potential therapeutic targets for anti-cancer therapy. Adv. Drug Deliv. Rev. 99, 180–185. Epub 2015/12/02. doi:10.1016/j.addr.2015.11.009
Leroy, B., Anderson, M., and Soussi, T. (2014). Tp53 mutations in human cancer: Database reassessment and prospects for the next decade. Hum. Mutat. 35 (6), 672–688. Epub 2014/03/26. doi:10.1002/humu.22552
Li, J. H., Liu, S., Zhou, H., Qu, L. H., and Yang, J. H. (2014). Starbase V2.0: Decoding mirna-cerna, mirna-ncrna and protein-rna interaction networks from large-scale clip-seq data. Nucleic Acids Res. 42, D92–D97. Epub 2013/12/04. doi:10.1093/nar/gkt1248
Liu, L., Wu, Y., Li, Q., Liang, J., He, Q., Zhao, L., et al. (2020). Mettl3 promotes tumorigenesis and metastasis through bmi1 M(6)a methylation in oral squamous cell carcinoma. Mol. Ther. 28 (10), 2177–2190. Epub 2020/07/06. doi:10.1016/j.ymthe.2020.06.024
Lopez-Yrigoyen, M., Cassetta, L., and Pollard, J. W. (2021). Macrophage targeting in cancer. Ann. N. Y. Acad. Sci. 1499 (1), 18–41. Epub 2020/05/24. doi:10.1111/nyas.14377
Maio, M., Blank, C., Necchi, A., Di Giacomo, A. M., Ibrahim, R., Lahn, M., et al. (2021). Neoadjuvant immunotherapy is reshaping cancer management across multiple tumour types: The future is now. Eur. J. Cancer 152, 155–164. Epub 2021/06/10. doi:10.1016/j.ejca.2021.04.035
Malfitano, A. M., Pisanti, S., Napolitano, F., Di Somma, S., Martinelli, R., and Portella, G. (2020). Tumor-associated macrophage status in cancer treatment. Cancers (Basel) 12 (7), 1987. Epub 2020/07/28. doi:10.3390/cancers12071987
Masarwy, R., Kampel, L., Horowitz, G., Gutfeld, O., and Muhanna, N. (2021). Neoadjuvant Pd-1/Pd-L1 inhibitors for resectable head and neck cancer: A systematic review and meta-analysis. JAMA Otolaryngol. Head. Neck Surg. 147 (10), 871–878. Epub 2021/09/03. doi:10.1001/jamaoto.2021.2191
Mody, M. D., Rocco, J. W., Yom, S. S., Haddad, R. I., and Saba, N. F. (2021). Head and neck cancer. Lancet 398 (10318), 2289–2299. Epub 2021/09/26. doi:10.1016/S0140-6736(21)01550-6
Ngambenjawong, C., Gustafson, H. H., and Pun, S. H. (2017). Progress in tumor-associated macrophage (Tam)-Targeted therapeutics. Adv. Drug Deliv. Rev. 114, 206–221. Epub 2017/04/30. doi:10.1016/j.addr.2017.04.010
Poeta, M. L., Manola, J., Goldwasser, M. A., Forastiere, A., Benoit, N., Califano, J. A., et al. (2007). Tp53 mutations and survival in squamous-cell carcinoma of the head and neck. N. Engl. J. Med. 357 (25), 2552–2561. Epub 2007/12/21. doi:10.1056/NEJMoa073770
Poulose, J. V., and Kainickal, C. T. (2022). Immune checkpoint inhibitors in head and neck squamous cell carcinoma: A systematic review of phase-3 clinical trials. World J. Clin. Oncol. 13 (5), 388–411. Epub 2022/06/07. doi:10.5306/wjco.v13.i5.388
Pu, Y., and Ji, Q. (2022). Tumor-associated macrophages regulate Pd-1/Pd-L1 immunosuppression. Front. Immunol. 13, 874589. Epub 2022/05/21. doi:10.3389/fimmu.2022.874589
Schantz, S. P., Shillitoe, E. J., Brown, B., and Campbell, B. (1986). Natural killer cell activity and head and neck cancer: A clinical assessment. J. Natl. Cancer Inst. 77 (4), 869–875. Epub 1986/10/01.
Schmidl, B., Siegl, M., Boxberg, M., Stogbauer, F., Jira, D., Winter, C., et al. (2022). Notch1 intracellular domain and the tumor microenvironment as prognostic markers in hnscc. Cancers (Basel) 14 (4), 1080. Epub 2022/02/26. doi:10.3390/cancers14041080
Sloas, C., Gill, S., and Klichinsky, M. (2021). Engineered car-macrophages as adoptive immunotherapies for solid tumors. Front. Immunol. 12, 783305. Epub 2021/12/14. doi:10.3389/fimmu.2021.783305
Thorsson, V., Gibbs, D. L., Brown, S. D., Wolf, D., Bortone, D. S., Ou Yang, T. H., et al. (2018). The immune landscape of cancer. Immunity 48 (4), 812–830.e14. Epub 2018/04/10. doi:10.1016/j.immuni.2018.03.023
Watermann, C., Pasternack, H., Idel, C., Ribbat-Idel, J., Bragelmann, J., Kuppler, P., et al. (2021). Recurrent hnscc harbor an immunosuppressive tumor immune microenvironment suggesting successful tumor immune evasion. Clin. Cancer Res. 27 (2), 632–644. Epub 2020/10/29. doi:10.1158/1078-0432.CCR-20-0197
Wu, P., Fang, X., Liu, Y., Tang, Y., Wang, W., Li, X., et al. (2021). N6-Methyladenosine modification of Circcux1 confers radioresistance of hypopharyngeal squamous cell carcinoma through Caspase1 pathway. Cell Death Dis. 12 (4), 298. Epub 2021/03/21. doi:10.1038/s41419-021-03558-2
Yao, Y., Kong, X., Liu, R., Xu, F., Liu, G., and Sun, C. (2022). Development of a novel immune-related gene prognostic index for breast cancer. Front. Immunol. 13, 845093. Epub 2022/05/14. doi:10.3389/fimmu.2022.845093
Yi, L., Wu, G., Guo, L., Zou, X., and Huang, P. (2020). Comprehensive analysis of the Pd-L1 and immune infiltrates of M(6)a rna methylation regulators in head and neck squamous cell carcinoma. Mol. Ther. Nucleic Acids 21, 299–314. Epub 2020/07/06. doi:10.1016/j.omtn.2020.06.001
Yu, D., Pan, M., Li, Y., Lu, T., Wang, Z., Liu, C., et al. (2022). Rna N6-methyladenosine reader Igf2bp2 promotes lymphatic metastasis and epithelial-mesenchymal transition of head and neck squamous carcinoma cells via stabilizing slug mrna in an m6a-dependent manner. J. Exp. Clin. Cancer Res. 41 (1), 6. Epub 2022/01/05. doi:10.1186/s13046-021-02212-1
Keywords: head and neck squamous cell carcinoma, m6A, immunotherapy, tumor immune microenvironment, prognosis
Citation: Yang Y, Cai Z, Huang K, Li M, Wang X, Lin Y, Chen S, Yang Z and Lin Z (2023) A composite score based on immune-related gene prognostic index and m6A risk score of head and neck squamous cell carcinoma. Front. Genet. 14:1061569. doi: 10.3389/fgene.2023.1061569
Received: 09 November 2022; Accepted: 31 January 2023;
Published: 09 February 2023.
Edited by:
Xin Hua, Shanghai Jiao Tong University, ChinaReviewed by:
Yao Xiao, Wuhan University, ChinaHongwei Chen, School of Public Health, Sun Yat-sen University, China
Copyright © 2023 Yang, Cai, Huang, Li, Wang, Lin, Chen, Yang and Lin. 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: Zhixiong Lin, zxlin5@qq.com