The latest RNA sequencing (RNA-seq) data, clinicopathologic information, and survival data of 546 HNSCC samples, including 502 tumor samples and 44 normal samples, were obtained from The Cancer Genome Atlas (TCGA) database through the Genomic Data Commons Data Portal (GDC; https://portal.gdc.cancer.gov/). RNA-seq data of 270 HNSCC tumor samples (GSE65858), corresponding survival information, and the Sequencing platform annotation information (GPL570) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The list of metabolism-related genes was collected and integrated from the Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolism-related pathways (https://www.kegg.jp/kegg/kegg1.html).

To obtain the metabolism-related differentially expressed genes (DEGs), we first obtained the DEGs (adj. p < 0.05, |log2FC| >1.0) between the 502 tumor samples and 44 normal samples of HNSCC from TCGA-HNSCC project using the limma package of R software. Then, by intersecting with the metabolism-related gene lists, the differentially expressed metabolism-related genes were extracted from the abovementioned DEGs. Volcano plots of DEGs and differentially expressed metabolism-related genes were plotted using the ggplot2 package of R. The clusterProfiler package of R was used to perform Gene Ontology (GO) and KEGG analysis of the obtained metabolism-related DEGs to explore the biological functions and processes involved in these genes.

The WGCNA was used to recognize the metabolism-related hub genes (35). The first step was to calculate the correlation coefficient (Pearson coefficient) between any two genes and construct the similarity matrix. To measure whether two genes have similar expression patterns, correlation coefficient weighted values were used during WGCNA. Next, the best “soft threshold” was determined graphically as β = 4 to ensure a scale-free network distribution. The similarity matrix and β value were adopted to construct an adjacency matrix and then transformed into a topological matrix with the topological overlap measure (TOM) describing the degree of association between genes. Here, 1-TOM was used as the distance to cluster the genes, and then the dynamic pruning tree was constructed to generate the coexpression gene modules. Finally, the similar gene modules were determined by merging modules whose distance is less than 0.25. The correlation between each module and HNSCC was calculated by the eigengene function in the “WGCNA” R package. The most significantly related module was selected for the follow-up study. In this module, the edges between two genes with weight >0.2 were selected to plot the network. The genes in the network with a degree ranking in the top 50 were the nominated hub genes. The best cutoff value of the impact of each hub gene on overall survival (OS) was obtained through the Survminer package of R. Only genes significantly related to survival (p < 0.05, log-rank test) were considered as the metabolism-related hub genes.

Among the obtained metabolism-related hub genes, the genes significantly affecting OS were identified by univariate Cox regression analysis and then utilized to construct the MRGPI by multivariate Cox regression analysis. The MRGPI of each sample was calculated as follows: MRGPI = ,where represents the regression coefficient between the gene (i) and HNSCC prognosis, and represents the expression level of the gene (i).

The HNSCC patients from TCGA-HNSCC project (training cohort) were divided into the high-MRGPI group and low-MRGPI group according to the median MRGPI. The ability of the MRGPI to discriminate the prognosis of patients was respectively evaluated by the Kaplan–Meier survival curve and the log-rank test of TCGA and GEO cohorts (validation cohort).

In addition, to verify the independent prognostic role of the MRGPI, the MRGPI and other clinicopathological factors were included in univariate and multivariate Cox regression analyses.

The gene set variation analysis (GSVA) package of R was used to calculate the common metabolic activity scores of each patient (36), including beta-alanine metabolism; fatty acid metabolism; glutathione metabolism; glycerolipid metabolism; nitrogen metabolism; purine metabolism; pyruvate metabolism; starch and sucrose metabolism; glutamate and glutamine metabolism; glycogen metabolism; glucose metabolism; and alanine, aspartate, and glutamate metabolism. Then, the difference in metabolic activities was compared between the high- and low-MRGPI groups.

Specifically, we firstly obtained the gene lists reflecting these metabolic activities ( Supplementary Table S1 ) from the KEGG (https://www.kegg.jp/) and Reactome (https://reactome.org/) pathway databases, which are open-source, open-access, manually curated, and peer-reviewed pathway databases and gained wide acceptance (37, 38). Then, we used the method of Hänzelmann et al. (36) to assess the GSVA enrichment scores of metabolic activities in each patient: 1) Input RNA-seq counts and a list of metabolic gene sets; 2) Gene expression level statistic: Kernel estimation of the cumulative density function (kcdf); 3) Rank ordered for each sample by the expression-level statistic; 4) Calculated the Kolmogorov–Smirnov-like rank statistic for every gene set; 5) Different score distributions; 6) Output a matrix containing pathway enrichment scores for each gene set and sample. The detailed method is available as a Bioconductor package for R under the name GSVA at http://www.bioconductor.org.

To discover the molecular characteristics of the different MRGPI subgroups, we firstly used the limma package of R to find out the DEGs between the high- and low-MRGPI groups. Then, the gene set enrichment analysis (GSEA) was used to enrich the signal pathways involved in different MRGPI subgroups based on the KEGG gene sets (c2. cp. kegg. v7.4) by using the clusterProfiler package of R [p < 0.05, false discovery rate (FDR) <0.25].

The data of genetic alterations (simple nucleotide variation) of HNSCC was obtained from TCGA GDC data portal, and R’s Maftools package was used to analyze the number and categories of gene mutations in the two MRGPI subgroups. The tumor mutation burden (TMB) of each sample was also calculated by the Maftools package of R.

There are several anticancer immune steps in the antitumor immune cycle, including Release of cancer antigens, Cancer antigen presentation, Priming and activation, Trafficking of immune cells to tumors, Infiltration of immune cells into tumors, Recognition of cancer cells by T cells, and Killing of cancer cells. These 7 antitumor immune processes constitute the antitumor immune cycle. The tracking tumor immunophenotype (TIP) website collected and curated gene lists that represent these seven antitumor immune steps, which are presented in Supplementary Table S2 ( , 39). Single-sample gene set enrichment analysis (ssGSEA) is a method that calculates the absolute enrichment degree of a given gene set based on these gene expression values of the sample sequencing data and then indirectly reflects the activity or enrichment degree of the given gene set (40). Therefore, we used the ssGSEA algorithm to score the seven antitumor immune processes during the antitumor immune cycle of 499 HNSCC tumor samples and to evaluate the tumor immunophenotype in high-MRGPI and low-MRGPI subgroups. Scores for each antitumor immune step for each HNSCC patient in this study were attached in Supplementary Table S3 . To further evaluate the difference in the level of immune cell infiltration between different MRGPI score subgroups, we used the CIBERSORT algorithm to quantitatively analyze the relative abundance of 22 immune cells in 499 HNSCC tumor samples (41). The specific information on the proportion of each immune cell infiltration in each patient was presented in Supplementary Table S4 .

To compare the immune state within the TME of the two MRGPI subgroups. The ESTIMATE score, immune score, and stromal score were calculated using the ESTIMATE algorithm. The ESTIMATE algorithm is a tool that uses gene expression data to predict tumor purity and the degree of infiltrating stromal/immune cells in tumor tissue. The ESTIMATE algorithm is based on ssGSEA and generates three scores: 1) stromal score (capturing the presence of stroma in tumor tissue), 2) immune score (representing the infiltration of immune cells in tumor tissue), and 3) ESTIMATE score (inferred tumor purity).

Then, the TIDE score, dysfunctional cytotoxic T cells, exclusion cytotoxic T cells, interferon gamma (IFNG) score and M2-TAMs, cancer-associated fibroblasts (CAFs), and myeloid-derived suppressor cell (MDSC) infiltrations were assessed by the TIDE platform to compare the antitumor and tumor immune escape abilities of the two MRGPI subgroups. TIDE is a computational framework that can use its built-in signatures and algorithm to calculate the degree of infiltration of dysfunctional T cells, M2-TAMs, CAFs, and MDSCs, as well as the IFNG score (reflecting cytotoxic killing of T cells). Then, the degree of infiltrating M2-TAMs, CAFs, and MDSCs in the TME was used to estimate the strength of the TME to exclude T cells, which was called the T-cell exclusion score. Finally, the TIDE score was obtained based on the T-cell dysfunction score and T-cell exclusion score to estimate the ability of the tumors to escape from immunity. Detailed principles and algorithms can be found in the articles published by Jiang et al. (42) and Fu et al. (43).

Finally, we referred to the gene sets published by He et al. reflecting many specific antitumor immunological functions (attached in Supplementary Table S5 ), 44) and used the ssGSEA algorithm to evaluate the immune and molecular characteristics between the different MRGPI subgroups.

In order to explore the prognostic value of the MRGPI for patients after immunotherapy, we respectively performed survival analysis in the IMvigor210 cohort and a Kidney renal clear cell carcinoma (KIRC) ICI cohort (45, 46). IMvigor210 cohort is a phase II clinical trial of atezolizumab (anti-PD-L1) therapy in patients with locally advanced or metastatic urothelial carcinoma. The primary follow-up information of this clinical trial is OS. The prognostic information and RNA-sequencing data of the IMvigor210 cohort were built into and can be publicly obtained by the IMvigor210CoreBiologies R package. The KIRC ICI cohort composes three prospective clinical trials of receiving nivolumab (anti-PD-1) therapy in clear cell renal cell carcinoma (ccRCC), including CheckMate 025, CheckMate 009, and CheckMate 010. The primary follow-up data of the study are progression-free survival (PFS) and OS in patients with ccRCC who received anti-PD-1 immunotherapy. The clinicopathological data, immunotherapy information, treatment response, prognostic information, and tumor sample RNA-sequencing data of this cohort were selflessly shared by Prof. Toni K. Choueir et al. (Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School) at https://doi.org/10.1038/s41591-020-0839-y. The MRGPI scores of patients in the IMvigor210 cohort and KIRC ICI cohort were calculated according to the following formula: (MRGPI) =, where represents the regression coefficient between the gene (i) and patient prognosis and represents the expression level of the gene (i).

In addition, time-dependent receiver operating characteristic curve (ROC) curve analyses were used to obtain the area under curve (AUC) of various prognostic models, which can reflect the model’s ability of prediction. The AUC of the MRGPI was compared with the AUC of two immune-related scores, TIDE and TIS. TIDE is a score obtained by comprehensively assessing the ability of the TME to exclude T cells and the infiltration levels of dysfunctional T cells within the TME. It is often used to evaluate the immune evasion ability of tumors and the potential efficacy of immunotherapy and can be calculated online (http://tide.dfci.harvard.edu/) (43). The TIS is a set of 18 genes that are highly correlated with clinical response to ICIs, including IFN-gamma signaling pathway, T cells, and Natural killer (NK) cell-related genes. The TIS score is calculated as the average log2 scale-normalized expression of 18 characteristic genes (42).

We explored the relationship between gene expression and drug sensitivity using the Gene Set Cancer Analysis (GSCA) online tool (http://bioinfo.life.hust.edu.cn/GSCA/#/drug). The GSCA tools collected the IC50 of 367 small molecules in 987 cell lines from Genomics of Drug Sensitivity in Cancer (GDSC) and the IC50 of 481 small molecules in 860 cell lines from Cancer Therapeutics Response Portal (CTRP) (47–49). Gene expression sequencing data of these cell lines were simultaneously collected from the Cancer Cell Line Encyclopedia (CCLE) database. Afterward, the mRNA expression data and drug IC50 data were merged. Pearson correlation analysis was performed to get the correlation between gene mRNA expression and drug IC50. The p-values with adjusted FDR were obtained. A negative correlation means that gene expression is suppressed indicating sensitivity to that drug and vice versa.

The Connectivity map (cMAP) database (https://portals.broadinstitute.org/cmap/) is a public online tool that was used to identify candidate inhibitors for the MRGPI subgroups based on gene expression profiles. Differential expression genes between the high- and low-MRGPI groups were used to query in cMAP to get a list of potential inhibitors (50). All compounds’ mechanism of action (MoA) and drug targets were downloaded from clue (https://clue.io).

Differential expression analysis was performed between HNSCC tumor samples (n = 502) and normal samples (n = 44), and finally, a total of 5,897 DEGs were obtained, including 4,473 upregulated genes and 1,424 downregulated genes ( Supplementary Figure S1A ). A total of 188 metabolism-related DEGs were then acquired by intersecting the lists of metabolism-related genes integrated from the KEGG database with the DEGs. Finally, 93 metabolism-related genes were found upregulated in tumor samples, while 95 genes were downregulated ( Supplementary Figure S1B ). The functional enrichment analysis showed that 188 metabolism-related DEGs were significantly associated with 690 GO terms and 54 KEGG pathways (details in Supplementary Table S6) , and the top 8 GO terms and KEGG pathways are shown in Supplementary Figures S1C, D .

To obtain the metabolism-related hub genes, WGCNA was carried out on the candidate genes (n = 188). The power value when the correlation coefficient between the connectivity K and the logarithm log(P(k)) reaches 0.85 is set as the β value. The optimal soft-thresholding power was 4 based on the scale-free network ( Supplementary Figure S2A ). Two modules were identified according to the average linkage hierarchical clustering and the optimal soft-thresholding power ( Supplementary Figures S2B, C ). According to the Pearson correlation coefficient between a module and sample feature for each module, the turquoise module is closely correlated with HNSCC tumors. Thus, the genes in the turquoise module were selected for further analysis. There were 23 genes and 99 edges of the turquoise module of the networks with a threshold weight >0.2 ( Supplementary Figure S2D ). The top 9 significantly enriched hallmark pathways for the genes in the turquoise module were shown in Supplementary Figure S2E (details in Supplementary Table S7 ). Among these 23 genes of the turquoise module, the expression of 17 metabolism-related hub genes was closely correlated with HNSCC patient OS as determined by Kaplan–Meier survival analysis, as shown in Supplementary Figure S3 (p < 0.05, log-rank test).

To further identify the metabolism genes that play as independent prognostic factors, multivariate Cox regression was carried out among the 17 metabolism-related hub genes. Finally, seven genes were recognized as independent factors, including HPRT1, AGPAT4, AMY2B, ACADL, CKM, PLA2G2D, and ADA. Therefore, the MRGPI of patients’ tumor samples was calculated by the formula: MRGPI = the expression level of HPRT1 * (0.57) + the expression level of AGPAT4 * (0.24) + the expression level of AMY2B * (-0.50) + the expression level of ACADL * (0.59) + the expression level of CKM * (0.07) + the expression level of PLA2G2D * (-0.21) + the expression level of ADA * (0.21).

The 499 HNSCC patients from TCGA were separated into the high-MRGPI group and low-MRGPI group according to the median MRGPI of all patients. Details of clinicopathological characteristics of the two MRGPI groups were exhibited in Supplementary Table S8 , and the distribution of these characteristics has no difference between the two groups. Univariate Cox analysis indicated that Stage and MRGPI are significantly associated with the prognosis of HNSCC patients ( Figure 1A ). Further multivariate Cox regression analysis confirms that Stage and MRGPI are independent prognostic factors after being adjusted by other clinicopathologic factors ( Figure 1B ). Importantly, the OS of patients with low MRGPI was significantly higher than that of patients with high MRGPI (p < 0.001, log-rank test; Figure 1C ). Similarly, the disease-specific survival (DSS) of patients with low MRGPI was also significantly better than that of patients with high MRGPI (p < 0.001, log-rank test; Supplementary Figure S4A ). Moreover, according to the previous treatment of HNSCC patients in TCGA database, we divided the patients into radiation therapy subgroups, non-radiation therapy subgroups, molecular targeted therapy subgroups, and non-molecular targeted therapy subgroups. Excitingly, the MRGPI had an excellent ability to distinguish patient outcomes in all 4 subgroups, i.e., patients in the low-MRGPI group had better OS than those in the high-MRGPI group (p < 0.05, log-rank test; Supplementary Figures S4C –S4F).

In addition, we used an HNSCC dataset (GSE65858; n = 270) from the GEO database to externally validate the predictive effect of the MRGPI score on the prognosis. Coinciding with the results from TCGA cohort, low-MRGPI patients have better OS than that of high-MRGPI patients (p = 0.013, log-rank test; Figure 1D ). Additionally, the patients with a low MRGPI have a longer disease-free interval (DFI) than that of high-MRGPI patients (p < 0.001, log-rank test; Supplementary Figure S4B ). The model discriminatory accuracy was self-verified in TCGA training cohort using the AUC, resulting in values of 0.66 (95% CI: 0.59–0.72), 0.66 (95% CI: 0.60–0.71), 0.69 (95% CI: 0.63–0.75), 0.67 (95% CI: 0.59–0.75), and 0.75 (95% CI: 0.63–0.87) at 1, 2, 3, 5, and 10 years, respectively ( Figure 1E ), and verified by the GEO validation cohort, resulting in values of 0.62 (95% CI: 0.53–0.72), 0.61 (95% CI: 0.52–0.69), 0.65 (95% CI: 0.56–0.74), and 0.76 (95% CI: 0.62–0.90) at 1, 2, 3, and 5 years, respectively ( Figure 1F ), which both reflected satisfactory accuracy.

To gain further biological insight into the metabolic character of the MRGPI subgroups, we estimated the scores of common metabolic pathways. The metabolic scores indicated that the patients of the high-MRGPI group have higher levels of Purine metabolism; Glutamate and glutamine metabolism; Glycogen metabolism; Glucose metabolism; and Alanine, aspartate, and glutamate metabolism, while the patients in the low-MRGPI group have higher Fatty acid metabolism level ( Figure 2A ).

GSEA was performed to determine the specific gene sets enriched in each MRGPI subgroup. Against the background of gene sets included in the hallmark dataset, the gene set enriched in the high-MRGPI subgroup contains pathways related to tumor epithelial-mesenchymal transition (EMT) conversion, glycolysis, and hypoxia ( Figure 2B ), while the gene set of the low-MRGPI sample is enriched in interleukin-2 (IL-2) channels ( Figure 2C ; p < 0.05, FDR <0.25). When using the KEGG database as a background, the gene sets enriched in the high-MRGPI subgroup contains pathways related to PRIMARY_IMMUNODEFICIENCY, FOCAL ADHESION, and immune response ( Supplementary Figure S5A ), while the gene set of the low-MRGPI sample is rich in cancer and tumor metastasis-related pathways ( Supplementary Figure S5B ). The detailed results of GSEA are listed in Supplementary Table S9 .

To further discover the different immunological natures brought about by metabolic characteristics, we explored the gene mutation landscapes of the different MRGPI subgroups. First, our results intuitively show that there are more samples with genetic mutations in the high-MRGPI subgroup than those in the low-MRGPI group (94.47% vs. 89.39%). Missense mutation accounted for the largest proportion (87.4%) of all mutation types, followed by Nonsense mutation (7.1%) and Frameshift deletions (2.9%) in the high-MRGPI subgroup. In the low-MRGPI group, the order of the proportions of the different mutation types is the same. Missense mutation is the most common mutation type (87.2%), followed by Nonsense mutation (7.0%) and Frameshift deletions (3.4%) in the high-MRGPI subgroup. Next, we explored the top 20 genes with the most frequent mutation rates in the two MRGPI subgroups ( Figure 2D ). Two MRGPI subgroups have obviously distinct mutant gene atlases. For example, although the gene with the highest mutation frequency in the two groups is TP53, the mutation rate of TP53 is as high as 74% in the high-MRGPI subgroup, which is much higher than the 51% mutation rate in the low-MRGPI group. Moreover, an important tumor suppressor gene, CDKN2A, which controls the cell cycle has an incidence rate of 20% in the high-MRGPI group and only 16% in the low-MRGPI group. Additionally, an important mutation that can cause the activation of the notch pathway, NOTCH1-mutation, has an incidence rate of 16% in the high-MRGPI group and only 14% in the low-MRGPI group.

The anticancer immune response can be conceptualized as a series of stepwise events, including the release of cancer antigens (step 1), cancer antigen presentation (step 2), initiation and activation (step 3), transportation of immune cells to the tumor (step 4), immune cells infiltrate the tumor (step 5), T cells recognize cancer cells (step 6) and kill cancer cells (step 7). The activity of these steps reflects the ability of the body’s immune cells to mobilize, transport, and kill, which together determine the success of the antitumor immune response. To analyze the impact of tumor metabolic activity on the seven-step anticancer immune response process of the 26 immune cells, we obtained signatures representing the seven-step anticancer immune process from the TIP database (attached in Supplementary Table S2 ) and used the ssGSEA algorithm to evaluate the ability of each antitumor step in patients. Then, we used the Wilcoxon test to compare the status of the seven antitumor immune processes in the different MRGPI subgroups ( Figure 3 ). We found that the ability to release cancer antigens was stronger in the high-MRGPI group ( Figure 3A ). The ability of immune cell priming and activation of the low-MRGPI patients is significantly superior to that of the high-MRGPI patients ( Figure 3C ). The total ability of trafficking immune cells to tumors is stronger in the low-MRGPI group than in the high-MRGPI group. In detail, the trafficking abilities of T cell, CD4 T cell, CD8 T cell, Th1 cell, DC, NK cell, and Th2 cell to tumors are stronger in the low-MRGPI patients than those in the high-MRGPI patients, while the trafficking abilities of Neutrophil and Treg cell to the tumor are weaker in the low-MRGPI patients ( Figure 3D ). In the next step, the degree of immune cell infiltration into tumors is naturally higher in the low-MRGPI group ( Figure 3F ). However, the capabilities of cancer antigen presentation, recognition of cancer cells by T cells, and killing of cancer cells between the two MRGPI groups showed no difference ( Figures 3B–G, H ).

To further assess the immune characteristics of the TME in the different MRGPI subgroups, we firstly used the CIBERSORT algorithm to evaluate the proportions of 22 immune cell infiltration of each HNSCC case in the different MRGPI subgroups ( Figure 4A ). We found that naive B cells, plasma cells, CD8 T cells, memory activated CD4 T cells, follicular helper T cells, Treg cells, gamma delta T cells, resting Mast cells, and Neutrophils were significantly more abundant in the low-MRGPI subgroup, while memory resting CD4 T cells, resting NK cells, M0 macrophages, and activated Mast cells showed a significantly high infiltration level in the high-MRGPI subgroup ( Figure 4A ).

Moreover, the ESTIMATE algorithm was used to evaluate the tumor microenvironmental components of each HNSCC sample. Finally, we calculated the ESTIMATE score (that infers tumor purity), immune score, and stromal score to compare the differences in the immune microenvironment between the low-MRGPI subgroup and the high-MRGPI subgroup. ESTIMATE score is the sum of the immune score and stromal score, which reflects the proportion of immune cells and stromal cells within the TME, and is also an indicator of tumor purity. A higher ESTIMATE score represents a higher proportion of immune cells and stromal cells within the TME and thus a lower proportion of tumor cells (lower tumor purity). That is, the ESTIMATE score is inversely proportional to tumor purity. We found that the low-MRGPI subgroup has a significantly higher ESTIMATE score (lower tumor purity) and immune score, which is consistent with the immune infiltration results above (p < 0.001) ( Figures 4B, C ). There was no difference in stromal scores between the two subgroups ( Figure 4D ).

In addition, we further used certain gene signatures that can reflect specific antitumor immunological functions to evaluate the immune and molecular characteristics between the different MRGPI subgroups. The results showed that there is better antigen presentation ability [activated dendritic cells (aDCs), antigen presenting cells (APC)], more chemokine receptors (CCR), stronger anti-inflammatory and cytotoxic effects, and more IFN responses in the low-MRGPI group ( Supplementary Figure S6A ).

As the MRGPI can successfully distinguish the OS rate of HNSCC patients in different subgroups and is related to multiple steps in the antitumor immune response processes, we naturally guessed whether the MRGPI can predict the potential clinical efficacy of ICI therapy.

First, a correlation analysis was performed between the MRGPI and immune checkpoint proteins, including PD-1 (PDCD1), PD-L1 (CD274), CTLA-4, lymphocyte activation gene-3 (LAG-3), T-cell immunoglobulin, and immunoreceptor tyrosine-based inhibition motif domain (TIGIT) in the 502 HNSCC patients. There were significant negative correlations between the MRGPI and these immune checkpoint proteins (p < 0.01), except for PD-L1, which was not associated with the MRGPI (p > 0.05) ( Supplementary Figures S6B–F ). Specifically, a low MRGPI correlates with a high expression of checkpoint proteins, while a high MRGPI correlates with a low expression of checkpoint proteins. Many studies have demonstrated that a higher expression of immune checkpoint proteins such as PD-1 and CTLA-4 is associated with higher response rates or better efficacy of immunotherapy (51, 52). Therefore, these results suggest that patients with a low MRGPI may be more likely to benefit from immunotherapy.

We then used the TIDE algorithm to assess the potential clinical efficacy of immunotherapy in the different MRGPI subgroups. We found that the low-MRGPI patients with higher immune cell infiltration had higher T-cell dysfunction scores, while the high-MRGPI patients with less immune cell infiltration level had higher T-cell exclusion scores ( Figures 5A, B ). Moreover, two signatures for assessing the ability to release IFN-γ, IFNG, and Merck18 both suggest that the patients of the low-MRGPI subgroup have a stronger IFN-γ release ability, which indirectly reflects that they have a stronger tumor-killing effect ( Figures 5C, D ). In addition, three kinds of suppressor cells that contribute to tumor immune escape, M2-TAMs), CAFs, and MDSCs, respectively, suggest that the high-MRGPI tumors have a stronger immunosuppressive microenvironment and a more formidable immune evasion tendency ( Figures 5E–G ). Finally, the above results are comprehensively analyzed to get the TIDE prediction score of each HNSCC patient. A higher TIDE prediction score reflected a higher potential for immune evasion, which suggested that the patients were less likely to benefit from ICI therapy. In our results, the TIDE scores of the low-MRGPI subgroup were surprisingly slightly higher than those of the high-MRGPI subgroup, implying that the low-MRGPI patients may not benefit from ICI therapy compared with the high-MRGPI patients ( Figure 5H ).

As the result of the ICI therapy benefit predicted by TIDE is inconsistent with the result of the antitumor immune response in the MRGPI subgroups, we further estimated the prognostic value of the MRGPI in two cancer cohorts that had received anti-PD-L1 therapy. In the IMvigor210 cohort of bladder cancer, patients with a low MRGPI have better immunotherapy efficacy and better OS prognosis ( Figure 5I ). In the renal clear cell carcinoma (KIRC) cohort treated with nivolumab anti-PD-1 monoclonal antibody, we found that patients with a low MRGPI have significantly better ICI therapy efficacy and higher OS rates than patients with a high MRGPI ( Figure 5J ). To validate the performance of the MRGPI on patient prognosis, we compared the performance of the MRGPI with TIDE and TIS in TCGA HNSCC cohort and IMvigor210 cohort. The AUC of the MRGPI was better than the AUC of TIS and TIDE at 3 years’ follow-up in TCGA cohort that included patients with comprehensive therapy ( Figure 5K ). However, in the cohort receiving ICI therapy, the AUC of the MRGPI was between TIS and TIDE at 1-year follow-up ( Figure 5L ). Hence, we considered that the MRGPI is an ideal predictive index whose predictive power for OS is better than TIDE in the ICI therapy cohort.

IC50, the 50% inhibitory concentration, represents the concentration of a drug that is required for 50% inhibition of tumor cells. It is commonly used as a measure of drug sensitivity. That is, the lower the IC50 value of a drug on tumor cells, the more sensitive the tumor cells are to this drug. For drug sensitivity analysis using the GDSC database, the IC50 values for the treatment of 987 tumor cells by 367 drugs, respectively, and transcriptome sequencing data for these tumor cell lines [based on RNA microarrays and normalized using a robust multi-array analysis (RMA) algorithm] were first obtained from the GDSC database. Then, the expression values of the candidate genes in each cell line were combined with the IC50 values of each drug for each cell line using the cell line name as a reference. Finally, correlation analysis was performed on gene expression values and drug IC50 values using Pearson correlation, and then, the correlation between candidate genes and each drug was calculated. For drug sensitivity analysis using the CTRP database, the IC50 values for the treatment of 860 tumor cells by 481 drugs were obtained from the CTRP database, and transcriptome sequencing data for these tumor cell lines [based on RNA microarrays and normalized, quantified by the RNA-Seq by Expectation Maximization (RSEM) algorithm] were obtained from the CCLE database. Then, the candidate gene expression levels and drug IC50 values were combined, and their correlations were calculated according to the above method. Because the drugs and tumor cells tested in the GDSC database and the CTRP database are different, the results obtained are also different. For ACADL, only one drug in both databases has a meaningful correlation of IC50 with its expression and is therefore not exhibited in Figures 6A, B . The AMY2B gene is not embodied in the GDSC database, so the correlation between AMY2B and drug sensitivity is not shown in Figure 6A . The drugs that sensitized the rest of the MRGPI genes were ranked by the integrated level of correlation coefficient and FDR (adjust. p-value). The top 30 ranked drugs were shown in Figures 6A, B . In the figure, the purple dots represent a negative correlation between the gene and the IC50 of the candidate drug. This means that the gene is positively correlated with the sensitivity of the candidate drug. We found that the ADA expression was positively correlated with sensitivity to classical chemotherapeutic agents, such vinblastine, camptothecin, gemcitabine, and methotrexate, as well as novel small-molecule inhibitors, such as AT-7519 [Cyclin-dependent kinases (CDK) inhibitor], AZD7762 [Checkpoint Kinase (CHK) inhibitor], and AZD8055 [mammalian target of rapamycin (mTOR) inhibitor]. PLA2G2D was negatively associated with the resistance to vincristine, methotrexate, GSK461364 (PLK1 inhibitor), parbendazole (microtubule inhibitor), and so on. HPRT1 was significantly positively correlated with sensitivity to CD-437 [Retinoic Acid Receptor gamma (RARγ) agonist], COL-3 [Matrix metalloproteinase (MMP) inhibitor], and manumycin A (antibiotic). CKM was significantly positively correlated with sensitivity to ciclopirox (antifungal agent), CR-1-31B (eIF4A inhibitor), and narciclasine. AGPAT4 was significantly positively correlated with resistance to KIN001-102 and ciclopirox and negatively correlated with resistance to camptothecin, chlorambucil, and triazolothiadiazine.

Potential compounds targeting the MRGPI subgroups were identified with the DEGs between the high- and low-risk groups by querying the cMAP database. The potential compounds specific to the MRGPI and the corresponding mechanism of action are shown in Figure 6C . The negative mean value indicated that changes in expression profiles of these drug-treated cell lines in the cMAP database were reversed from those in the high-risk group and the positive value was just the opposite. Therefore, alimemazine, doxylamine, nimesulide, NU-1025, and ondansetron may serve as potential inhibitors targeting the MRGPI in the high-risk group.