To construct an immune-related gene prognostic index (IRGPI) for head and neck squamous cell carcinoma (HNSCC) and clarify the molecular and immune characteristics and the benefit of immune checkpoint inhibitor (ICI) therapy in IRGPI-defined subgroups of HNSCC.
On the basis of The Cancer Genome Atlas HNSCC immune dataset (n = 546), 22 immune-related hub genes were identified by weighted gene coexpression network analysis. Three genes were identified to construct an IRGPI by using the Cox regression method and validated with the Gene Expression Omnibus (GEO) dataset (n = 270). Afterward, the molecular and immune characteristics and the benefit of ICI therapy in IRGPI-defined subgroups were analyzed.
The IRGPI was constructed on the basis of SFRP4, CPXM1, and COL5A1 genes. IRGPI-high patients had a better overall survival than IRGPI-low patients, consistent with the results in the GEO cohort. The comprehensive results showed that a high IRGPI score was correlated with DNA repair–related pathways; low TP53 mutation rate; high infiltration of CD8 T cells, CD4 T cells, and M1 macrophages; active immunity and less aggressive phenotypes; and more benefit from ICI therapy. In contrast, a low IRGPI score was associated with cancer and metastasis-related pathways; high TP53 and PIK3CA mutation rate; high infiltration of B cells, M0 macrophages, and M2 macrophages; suppressive immunity and more aggressive phenotypes; and less benefit from ICI therapy.
IRGPI is a promising biomarker to distinguish the prognosis, the molecular and immune characteristics, and the immune benefit from ICI therapy in HNSCC.
Immune checkpoint inhibitor (ICI) therapy has proven to be a promising treatment for head and neck squamous cell carcinoma (HNSCC). However, the major limitation of ICI therapy is the low response rate of patients to ICI therapy and high out-of-pocket medical expenses. There is an urgent need to comprehensively understand the tumor microenvironment of HNSCC and identify a valuable biomarker for predicting the benefit of immunotherapy for patients with HNSCC. In this study, we developed an immune-related genetic prognostic index (IRGPI) and analyzed its role in discriminating different molecular and immune characteristics and outcomes of HNSCC. Integrated analysis of multi-omics data demonstrated that IRGPI-high patients were characteristic of active immune response and less aggressive tumor phenotype, had a longer overall survival time, and could have more benefit from ICI therapy, while the opposite was reported for IRGPI-low patients. Therefore, IRGPI as an immune-related prognostic biomarker might bring some potential implications for the immunotherapy strategy in HNSCC.
Compared with traditional therapies, immune checkpoint inhibitor (ICI) treatments, such as those targeting programmed death-ligand 1 (PD-L1), programmed death 1 (PD1), and CTL–associated protein 4 (CTLA4), have shown significant benefits in survival (1–5). In head and neck squamous cell carcinoma (HNSCC), anti-PD1 therapy has been proven to be a promising treatment for recurrent/metastatic patients (6–9). However, a major limitation is the low response rate of patients to ICI therapy (9, 10). Multiple factors including the immune microenvironment (TME) can affect ICI effectiveness, and few biomarkers can predict patient prognosis well (11). Identification of potential prognostic markers associated with treatment benefit can allow individualization of immunotherapy for patients with HNSCC. Unfortunately, we know little about TME of HNSCC, and effective prognostic and therapeutic indicators are urgently needed.
In this study, we sought to develop a prognostic marker for HNSCC, which could predict the prognosis of both conventional therapy and immunotherapy. We focused on all immune-related genes in the transcriptome data of HNSCC and screened immune-related hub genes related to patient prognosis by weighted gene coexpression network analysis (WGCNA) to construct an immune-related gene prognostic index (IRGPI). We then characterized the molecular and immune profile of IRGPI, examined its prognostic ability in immunotherapy patients, and compared it with other biomarkers, tumor immune dysfunction and exclusion (TIDE), and tumor inflammation signature (TIS). The results showed that IRGPI was a promising prognostic biomarker for patients receiving traditional treatment and immunotherapy.
Materials and Methods
Patients and datasets
RNA sequencing (RNA-seq) data of 546 HNSCC samples, including 502 cancer samples and 44 para-cancer samples, and their clinicopathologic information were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/projects/TCGA-HNSC). RNA-seq data of 270 HNSCC samples (GSE65858) and the survival information were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The lists of immune-related genes were downloaded from the ImmPort (https://www.immport.org/shared/home) and InnateDB (https://www.innatedb.com/) databases.
Information about the relationship between mRNA and transcriptional factors (TF), miRNA, and lncRNA was downloaded from the RAID 2.0 (http://www.rna-society.org/), TRRUST v2 (https://www.grnpedia.org/trrust/), and STRING (https://www.string-db.org/) databases. Gene mutation information was downloaded from the cBioPortal database (http://www.cbioportal.org/).
The molecular subtypes of the 502 HNSCC samples were kindly provided by Vonn Andrew Walter (Lineberger Comprehensive Cancer Center, University of North Carolina, Chapel Hill, NC) and Pierre Saintigny (Department of Translational Research and Innovation, Univ Lyon, Université Claude Bernard Lyon 1, France). The immune subtypes of the 502 HNSCC samples were kindly provided by Yu-Pei Chen (Sun Yat-sen University Cancer Center, Guangzhou, China).
Identification of immune-related hub genes
On the basis of the RNA-seq data of HNSCC samples (502 tumors vs. 44 normal samples) obtained from TCGA, lists of differentially expressed genes (P < 0.05, |log2FC| > 0.585) were identified by using the limma package of R. After consideration in the context of the immune-related gene lists obtained from ImmPort and InnateDB, immune-related differentially expressed genes were obtained and analyzed by using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses with the clusterProfiler package of R.
Then, WGCNA was performed to identify hub genes. First, the similarity matrix was constructed by using the expression data by calculating the Pearson correlation coefficient between two genes. Next, the similarity matrix was transformed into an adjacency matrix with a network type of signed and a soft threshold of β = 5 and then transformed into a topological matrix with the topological overlap measure (TOM) describing the degree of association between genes. 1-TOM was used as the distance to cluster the genes, and then the dynamic pruning tree was built to identify the modules. Finally, we identified seven modules by setting the merging threshold function at 0.25. On the basis of the genes of significantly related modules (the blue and turquoise modules), the edges between two genes with weight > 0.2 were used to construct the network. The genes in the network with a degree ranking in the top 50 were the hub genes. The best cut-off value for each hub gene for overall survival (OS) was obtained by the maxstat package of R, and 22 significantly survival-associated immune-related hub genes were selected for further analysis (P < 0.05, log-rank test).
To reveal relevant genetic alterations, somatic mutations of 22 immune-related hub genes were analyzed by using the ComplexHeatmap package of R. To determine the potential regulatory mechanisms, the roles of these genes were analyzed in TF, miRNA, and lncRNA regulatory networks. KEGG pathway analysis of the genes in the network was performed using Cytoscape software (P < 0.05).
Construction and validation of the IRGPI
Among 22 immune-related hub genes, the genes significantly affecting OS were identified and used to construct an IRGPI by multivariate Cox regression analysis. The IRGPI of each sample was calculated by multiplying the expression values of certain genes by their weight in the Cox model and then adding them together. The prognostic power of the IRGPI was evaluated by Kaplan–Meier (K–M) survival curves with log-rank tests with both TCGA and GEO cohorts. To validate the independent prognostic value of IRGPI, univariate and multivariate Cox regression analyses were performed.
Comprehensive analysis of molecular and immune characteristics and ICI therapy in different IRGPI subgroups
In the signaling pathway analysis, differential expression analysis was first performed on all genes to analyze the samples with high (n = 251) and low (n = 251) IRGPI score using the limma package of R. Enrichment analysis to determine the signaling pathways in which the differentially expressed genes are involved was then carried out by using the gene set enrichment analysis (GSEA) method based on the KEGG and HALLMARK gene sets with the clusterProfiler package of R (P < 0.05 and FDR < 0.25). Single sample GSEA (ssGSEA) analysis was then performed on several representative gene sets with the GSVA package of R, and K–M survival curves were used to explore differences in survival. In the gene mutation analysis, information on genetic alterations was obtained from the cBioPortal database, and the quantity and quality of gene mutations were analyzed in two IRGPI subgroups by using the Maftools package of R. Correlation analysis were performed between IRGPI score and PD-L1 expression and total mutation burden (TMB).
To identify immune characteristics of 502 HNSCC samples, their expression data were imported into CIBERSORT (https://cibersort.stanford.edu/) and iterated 1,000 times to estimate the relative proportion of 22 types of immune cells. Then, we compared the relative proportions of 22 types of immune cells and clinicopathologic factors between the two IRGPI subgroups, and the results are presented in a landscape map.
To further define the immune and molecular function between the IRGPI subgroups, we performed ssGSEA on certain gene signatures and compared the score between two IRGPI subgroups (12–15).
To explore the prognostic value of IRGPI in patients after immunotherapy, we performed survival analyses in two urothelial cancer cohorts treated with anti-PD-L1 (16, 17). Moreover, we performed time‐dependent ROC curve analyses to obtain the AUC and compared the prognostic value among IRGPI, TIDE, and TIS with the timeROC package of R. TIDE score was calculated online (http://tide.dfci.harvard.edu/) and TIS score was calculated as an average value of log2-scale normalized expression of the 18 signature genes (18).
An independent t test was performed to compare continuous variables between two groups. Categorical data were tested using the χ2 test. TIDE score between groups was compared by the Wilcoxon test. Univariate survival analysis was performed by K–M survival analysis with the log-rank test. Multivariate survival analysis was performed using the Cox regression model. A two-sided P < 0.05 was considered significant.
Immune-related hub genes
In differential expression analysis (502 tumors vs. 44 normal samples), a total of 3,877 differentially expressed genes were obtained, of which 2,681 genes were upregulated and 1,196 genes were downregulated in the tumor samples compared with normal samples (Supplementary Fig. S1A). By intersecting these genes with the lists of immune-related genes obtained from ImmPort and InnateDB, 1,131 differentially expressed immune-related genes were obtained, of which 860 genes were upregulated and 271 were downregulated in tumor samples compared with normal samples (Supplementary Fig. S1B). The functional enrichment analysis showed that 1,131 differentially expressed genes were significantly associated with 127 GO terms and 96 KEGG pathways (details in Supplementary Table S1), and the top 10 GO terms and KEGG pathways are shown in Supplementary Fig. S1C and S1D.
To obtain the immune-related hub genes, WGCNA analysis was carried out on the candidate genes (n = 1,131). The logarithm log(k) of the node with connectivity K was negatively correlated with the logarithm log(P(k)) of the probability of the node, and the correlation coefficient was greater than 0.85. The optimal soft-thresholding power was 5 based on the scale-free network (Supplementary Fig. S2). Seven modules were then identified on the basis of the average linkage hierarchical clustering and the optimal soft-thresholding power (Supplementary Fig. S3A and S3B). A total of 1,131 genes were allocated to seven modules. According to the Pierson correlation coefficient between a module and sample feature for each module, blue and turquoise modules closely correlated with HNSCC tumors, the genes in these modules were selected for further analysis. The top 10 significantly enriched GO terms and KEGG pathways for the genes of the blue and turquoise modules are shown separately in Supplementary Fig. S3C and S3D (details in Supplementary Table S2). There were 35 genes and 248 edges for the blue module and 83 genes and 1191 edges of the turquoise module of the networks with a threshold weight > 0.2 (Supplementary Fig. S3E and S3F). We obtained the top 50 immune-related hub genes with a threshold degree of > 20. The expression of 22 immune-related hub genes was closely correlated with HNSCC patient OS as determined by K–M analysis, shown in Supplementary Fig. S4 (P < 0.05, log-rank test).
We then explored the characteristics of the 22 immune-related hub genes. As showed in Supplementary Fig. S5A, most of the 22 immune-related hub genes had amplification, deep deletion, and missense mutations, and the mutation rates of FN1, COL11A1, VCAN, and COL5A1 were greater than 5%. In the TF and ncRNA analysis, there were 26 interaction pairs between hub genes and miRNAs, 59 interacting pairs between miRNAs and lncRNAs, and 37 interaction pairs between hub genes and TFs. After removing redundancy, there were 118 interaction pairs and 77 genes in the network (Supplementary Fig. S5B). By using the ClueGO plugin of Cytoscape, KEGG pathway enrichment analysis of the genes in the regulatory network was carried out, and it was found that they were significantly enriched in the AGE–RAGE signaling pathway associated with diabetic complications, transcriptional misregulation in cancer, and longevity-regulating pathways (Supplementary Fig. S5C).
Survival outcomes in different IRGPI groups
To determine the independent prognostic genes, multivariate Cox regression analysis for OS was performed among the 22 immune-related hub genes. As showed in Fig. 1A and B, only three genes (SFRP4, CPXM1, and COL5A1) significantly affected the OS of patients with HNSCC. Then, we constructed a prognostic index for all cancer samples calculated by the formula IRGPI = expression level of SFRP4*(−0.18) + expression level of CPXM1*0.29 + expression level of COL5A1*(−0.23).
Detailed clinicopathologic characteristics of 502 patients with HNSCC in TCGA cohort were shown in Supplementary Table S3 and most of them (except HPV status) were equally distributed between the two IRGPI subgroups. Univariate Cox regression analysis showed that age, IRGPI, HPV status, and stage were significantly associated with the prognosis of HNSCC. Multivariate Cox regression analysis confirmed that IRGPI was an independent prognostic factor after adjusted for other clinicopathologic factors (Fig. 1C; Supplementary Table S4).
Taking the median IRGPI as the cut-off value, IRGPI-high patients had a better OS than IRGPI-low patients (P = 0.022, log-rank test; Fig. 1D). Then, the role of IRGPI was validated by using the GSE65858 (n = 270) HNSCC dataset. As showed in Fig. 1E, the patients in the IRGPI-high subgroup had a significantly better prognosis than those in the IRGPI-low subgroup, consistent with the result of TCGA dataset (P = 0.013, log-rank test).
Molecular characteristics of different IRGPI subgroups
GSEA was performed to determine the gene sets enriched in different IRGPI subgroups. The gene sets of the IRGPI-high samples were enriched in DNA repair and immune response-related pathways (Fig. 2A), while the gene sets of the IRGPI-low samples were enriched in cancer and tumor metastasis-related pathways (Fig. 2B; P < 0.05, FDR < 0.25). Detailed results of GSEA are listed in Supplementary Table S5.
Next, we analyzed gene mutations to gain further biological insight into the immunologic nature of the IRGPI subgroups. We found significantly higher mutation counts in the IRGPI-high subgroup than in the IRGPI-low subgroup (P = 0.002, t test). Missense variations were the most common mutation type, followed by nonsense and frameshift deletions. We then identified the top 10 genes with the highest mutation rates in the IRGPI subgroups (Fig. 2C). The mutation rates of TP53, TTN, CDKN2A, CSMD3, FRG1B, MUC16, NOTCH1, and FAT1 were higher than 15% in both groups. The mutation of the SYNE1 and NSD1 genes was more common in the IRGPI-high subgroup, while the mutation of PIK3CA and LRP1B genes was more common in the IRGPI-low subgroup.
Next, we explore the relationship between IRGPI score and PD-L1 expression and TMB. As a result, the IRGPI score was slightly correlated with TMB (r = 0.122, P = 0.006), while the IRGPI score was not significantly correlated with PD-L1, shown in Supplementary Fig. S6.
Immune characteristics of different IRGPI subgroups
To analyze the composition of immune cells in different IRGPI subgroups, we used the Wilcoxon test to compare the distribution of immune cells in different IRGPI subgroups. We found that CD8 T cells, naïve CD4 T cells, activated memory CD4 T cells, resting natural killer (NK) cells, and M1 macrophages were more abundant in the IRGPI-high subgroup, while naïve B cells, resting memory CD4 T cells and M2 macrophages were more abundant in the IRGPI-low subgroup (Fig. 3A). Characteristics related to the immune landscape, including the clinicopathologic characteristics of different IRGPI subgroups, are displayed in Fig. 3B.
We then applied certain gene signatures to define the immune and molecular function between different IRGPI subgroups. As a result, there were more CD8 T cells, MHC class I, damage repair in the IRGPI-high subgroup, while there were more immunosuppressive cells and signals and tumor and metastasis-related signals in IRGPI-low subgroup.
We further investigated whether the prognostic value of IRGPI resulted from better immune control or less aggressive cancer growth. As shown in Supplementary Fig. S7, we could find that patients with a higher score of epithelial–mesenchymal transition (EMT), TGFβ, and WNT-related signals had a worse outcome, while patients with more P53 mutation, more CD8 T cells, and more M1 macrophages had a better prognosis. Therefore, we suggested that the prognostic value of IRGPI might result from both better immune control and less aggressive cancer growth.
Relationship between IRGPI grouping and other immune and molecular subtypes
A HNSCC immune subtype classification has described the immune landscape of HNSCC according to the tumor and stromal compartments and has summarized three immune subtypes: nonimmune subtype, immune-exhausted subtype, and immune-active subtype (19). We could find from Fig. 4A that the proportion of nonimmune samples was almost equally distributed between the two groups, but there were more immune-active samples and fewer immune-exhausted samples in the IRGPI-high subgroup than in the IRGPI-low subgroup (P < 0.001, χ2 test). Then 203 immune samples were further classified according to a pan-SCC immune subtype (20). As showed in Fig. 4B, there were more IS1 and IS5 subtypes in the IRGPI-low subgroup, while more IS4 subtype in IRGPI-high subgroup (P < 0.001, χ2 test).
In large genomic profiling studies of HNSCC, four distinct molecular subtypes have been consistently reported, namely, atypical, basal, classical, and mesenchymal. We then focused on the distribution of the molecular subtypes in the IRGPI groups (21–23). In our study, the IRGPI-low subgroup comprised 16% atypical samples, 22% basal samples, 19% classical samples, and 43% mesenchymal subtype samples, while the IRGPI-high subgroup comprised 36% atypical samples, 37% basal samples, 15% classical samples, and 13% mesenchymal samples (Fig. 4C). The classical subtype was almost evenly distributed between the two groups, but there were more atypical and basal samples and fewer mesenchymal samples in the IRGPI-high subgroup than in the IRGPI-low subgroup (P < 0.001, χ2 test).
The benefit of ICI therapy in different IRGPI subgroups
We then used TIDE to assess the potential clinical efficacy of immunotherapy in different IRGPI subgroups. Higher TIDE prediction score represented a higher potential for immune evasion, which suggested that the patients were less likely to benefit from ICI therapy. In our results, the IRGPI-high subgroup had a lower TIDE score than the IRGPI-low subgroup, implying that IRGPI-high patients could benefit more from ICI therapy than IRGPI-low patients. (Fig. 5A). Besides, a higher TIDE prediction score was associated with a worse outcome. So the IRGPI-high subgroup with a low TIDE score might have a better prognosis than the IRGPI-low group with a high TIDE score. Also, we found that the IRGPI-high subgroup had a higher microsatellite instability (MSI) score, while the IRGPI-low subgroup had a higher T-cell exclusion score, but there was no difference in T-cell dysfunction between the two subgroups. Moreover, we assessed the prognostic value of IRGPI in two urothelial cancer cohorts with anti-PD-L1 therapy (16, 17). As shown in Fig. 5B and C, we could find that IRGPI-high patients had better OS than IRGPI-low patients. We could also find that the performance of TIDE was inconsistent in Alexandra and colleagues urothelial cancer cohort (ref. 17; Fig. 5D) and Mariathasan and colleagues urothelial cancer cohort (ref. 16; Fig. 5E), the AUCs for TIS were better at 12 and 18 months follow-up and that of IRGPI was better at 18 and 24 months follow-up. So we suggested that the predictive value of IRGPI was comparable with 18-gene T cell-inflamed signature (TIS) and TIDE for OS in both cohorts (16, 17).
ICI therapy has proven to be an effective treatment for patients with recurrent or refractory HNSCC (6–9). Given that the overall response rate to ICI therapy is still very low (9, 10), the identification of patients who can benefit most from these therapies is crucial. After years of evaluating different prognostic markers in HNSCC, we still have not found a validated biomarker for predicting response to immunotherapy and OS. This highlights the need to identify a prognostic biomarker for immunotherapy in HNSCC.
Multiple genes affect the TME of tumors, and WGCNA is a virtual approach that can help identify candidate immune-related biomarkers or therapeutic targets. In our study, based on HNSCC immune gene datasets, we used WGCNA to identify 22 immune-related hub genes affecting patient OS and constructed the IRGPI based on three genes (SFRP4, CPXM1, and COL5A1), which were independent prognostic factors for OS. The IRGPI proved to be a valid prognostic immune-related biomarker for HNSCC, with better survival in IRGPI-high patients and worse survival in IRGPI-low patients in both TCGA and GEO cohorts.
In the GEO cohort, although the two curves intersect at a 20% risk of death, it did not follow that 20% of patients in the IRGPI-low group have better survival than those in the IRGPI-high group. First, both patients in the IRGPI-high group and the IRGPI-low group met the proportional hazard assumption (P = 0.25). So their effect of covariates on survival would not change with time. This meant that although the curves appeared to intersect, the risk of death for patients in the IRGPI-low group was greater than that of patients in the IRGPI-high group every time point. In addition, at the intersection of the curves, the last patient died in the IRGPI-high group and there were only 4 patients left in the IRGPI-low group after the intersection. So 20% risk of death did not mean that 20% of patients died. Moreover, the patients in the IRGPI-high group had a shorter follow-up time than the 4 patients left in the IRGPI-low group after the intersection point, so when the last IRGPI-high patient died and the curve of the IRGPI-high group dropped to 0, the curve inevitably intersected the curve of the IRGPI-low group. If the follow-up time was prolonged, and there were still patients surviving in the IRGPI-high group, then the two curves might not have intersected.
IRGPI was made up of three genes, SFRP4, CPXM1, and COL5A1. Collagen type V alpha 1 chain (COL5A1), is a minor connective tissue component of nearly ubiquitous distribution, distributed in the collagen-containing extracellular matrix. A study of malignant melanoma found that mutations in COL5A1 were associated with infiltration of CD8 T cells and activated NK cells and better immunotherapy survival (24). As for secreted frizzled-related protein 4 (SFRP4), mainly distributed in the cytoplasm of various cells, it was predicted to interact with SIRT1, which has been found to accumulate at the promoters of TNFα and IL1β in response to TLR4 signaling, promote termination of NFκB-dependent transcription and suppress the immune response (25). Besides, SFRP4 has been reported to be positively correlated with Treg cells infiltration while its downregulation in tumor cells impaired the production of cytokines and the recruitment of T cells (26). And it was found to be coexpressed with the EMT-linked markers, patients overexpressing SFRP4 presented with poor OS (27). As for carboxypeptidase X, M14 family member 1 (CPXM1), it was distributed in the extracellular space and involved in protein processing and cell–cell interactions. There was experimental evidence indicating that expression of CPXM1 was epigenetically regulated in breast cancer and it might act as a tumor suppressor gene (28). In the calculation formula of IRGPI, the coefficients of COL5A1 and SFRP4 were negative numbers, while the coefficient of CPXM1 was positive number. Therefore, there was negative relationship between IRGPI and COL5A1, and SFRP4, while positive relationship between IRGPI and CPXM1. In summary, IRGPI was a biomarker associated with active immunity and tumor suppression.
To gain further biological insight into the immunologic nature of the IRGPI subgroups, we then studied gene mutations of different IRGPI subgroups. Missense variations were the most common, followed by nonsense and frameshift deletions, as reported previously (29). The largest difference in mutations between groups was in TP53 mutations, which were more common in IRGPI-low samples than IRGPI-high samples (77% vs. 64%). TP53 mutation is not only the single most common genetic event in cancer but is also linked with more aggressive disease and poorer patient outcomes in many cancers (30, 31), particularly HNSCC (32). TP53 can affect the cancer cell cycle through the p53/TGF β signaling pathway. In addition, there was a higher rate of PIK3CA mutation in the IRGPI-low subgroup than the IRGPI-high subgroup, which could mean that IRGPI-low HNSCCs promote proliferation through the PI3K–AKT signaling pathway (33). Therefore, IRGPI-low patients with high TP53 and PIK3CA mutations have a worse outcome than IRGPI-high patients with low TP53 and PIK3CA mutations, in agreement with our survival results.
Next, we explore the relationship between IRGPI and known predictive biomarkers for immunotherapy, such as PD-L1 and TMB. PD-L1+ tumors, in general, tend to better respond to anti-PD-1/PD-L1 therapies than PD-L1− tumors (34, 35). However, we could find inconsistent results in HNSCC that CHECKMATE-141 failed to show a significant correlation between PD-L1 expression and tumor response or survival when evaluating anti-PD-1 therapy in the platinum-refractory recurrent/metastatic setting (36). The lack of uniformity in the assays and the variability in the thresholds used to define PD-L1 positivity might be the most important reasons. Moreover, we thought that the intensity and location of PD-L1 expression detected by IHC were more valuable than the PD-L1 expression value measured by transcriptome data. Therefore, further researches were needed to clarify the relationship between PD-L1 and IRGPI. Besides, TMB has been recently evaluated as a potential biomarker predicting response to ICI therapy in prospective clinical trials and across many tumor types including HNSCC (37). Here, we found that the IRGPI score had a slight correlation with TMB, which suggested that TMB could help explain why IRGPI affects immunotherapy prognosis to a certain extent, but there were more other possible mechanisms involved in it.
Understanding the landscape of the TME could help in finding new ways to treat HNSCC or to alter the TME to improve the effectiveness of immunotherapies. The composition of some immune cells was different between two IRGPI subgroups. While cytotoxic CD8 T cells, CD4 T cells, and M1 macrophages were more enriched in the IRGPI-high subgroup, B cells and M0 and M2 macrophages were more common in the IRGPI-low subgroup. A substantial body of research has revealed that dense infiltration of T cells, especially cytotoxic CD8 T cells, indicates a favorable prognosis (38–40). In most tumors, M2 macrophages, a predominant subtype of macrophages, have been proven to be related to chronic inflammation and favor tumor growth and development of an invasive phenotype, and these cells have been associated with a poor outcome in breast, bladder, ovarian, gastric, and prostate cancers (39, 41, 42). Conversely, a high density of M1 macrophages might correlate with acute inflammation and imply a favorable prognosis among patients with NSCLC, HCC, or ovarian or gastric cancers (39, 41, 42). The results of our study support these conclusions. Moreover, we found that the IRGPI-high samples had the stronger ability of damage repair, while the IRGPI-low samples had more immunosuppressive cells and signals and tumor and metastasis-related signals, which implied that IRGPI-low subgroup was characteristics of immunosuppression and active tumor progression.
Integrated with other molecular and immune subtypes, IRGPI grouping could distinguish different molecular and immune subtypes of HNSCC. In terms of a HNSCC immune subtype classification, there were fewer patients with immune-exhausted subtype and more patients with immune-active subtype in the IRGPI-high subgroup than in the IRGPI-low subgroup. According to Chen study, the immune-active and immune-exhausted subtypes have significant differences in M1/M2 macrophages, B cells, and cytolytic activity, but no significant differences in T cells, CD8 T cells, and cytotoxic cells. And the immune-active subtype shows a close relationship with immune-active pathways and gene sets, and the immune-exhausted subtype is characterized by tumor-promoting signals, such as activation of the WNT/TGFβ signaling pathway suppressing the host immune response (19). As such, IRGPI-high patients may have a more vigorous immune response to the development of tumorigenesis and tumor progression and therefore benefit more from ICI therapy than IRGPI-low patients. Furthermore, according to a pan-SCC immune subtype, there were more IS1 and IS5 in the IRGPI-low subgroup, while more IS4 in the IRGPI-high subgroup. IS1 demonstrated an intermediate immune infiltration toward an immune-suppressive phenotype with elevated expression of reactive stroma and TGFβ modules and IS5 had the highest expression for most gene modules such as inflammation, reactive stroma, and TGFβ (except IFNγ), implying an immune-hot but suppressive microenvironment (20). In contrast, IS4 had a more favorable antitumor immune response with the highest T cell and IFNγ gene expression and low reactive stroma and TGFβ, suggesting a favorable immune-activated phenotype (20). The results suggested that IRGPI-high patients were characteristics of active immunity, while IRGPI-low patients were characteristics of immunosuppression, which was in agreement with our previous results. In terms of molecular subtypes, which were classified mainly based on genome-wide profiling (21–23), our study revealed that the less aggressive atypical and basal subtypes were more common in the IRGPI-high subgroup, while the more aggressive mesenchymal subtype was more common in the IRGPI-low subgroup. The results of the three research were consistent with ours, that the IRGPI-high subgroup was characterized by an active immune response and less tumor aggressiveness, while the IRGPI-low subgroup was characterized by an immune-suppressive response and more tumor aggressiveness.
Interestingly, IRGPI‐based differences in the TME might reflect different immune benefit from ICI therapy (anti-PD1 and anti-CTLA4) identified with TIDE. The TIDE prediction score correlated with T-cell dysfunction in CTL-high tumors and T-cell exclusion in CTL-low tumors and thus represents two different mechanisms of immune escape. In our study, IRGPI-low patients had less CTL infiltration and higher TIDE and T-cell exclusion score (but not T-cell dysfunction score) than IRGPI-high patients, so their lower ICI response might be due to immune evasion via T-cell exclusion (43). In contrast, the IRGPI-high subgroup had a higher MSI score and lower TIDE score than the IRGPI-low subgroup, which suggested that these patients had lower levels of immune escape and more MSI. Previous studies have proven that MSI is ubiquitous in HNSCC (44), and the high mutational burden that results from MSI renders tumors immunogenic and sensitive to anti-PD1 therapy (45, 46). To further validate the prognostic value of IRGPI, we performed survival analysis in two urothelial cancer cohorts receiving anti-PD-L1 therapy. We found that IRGPI could differentiate different outcomes in patients treated with anti-PD-L1 therapy.
It has been reported that some biomarkers, such as TIDE and TIS, could predict patient response to immunotherapy. TIDE is a creative computational method to identify factors that underlie two mechanisms of tumor immune escape: the induction of T-cell dysfunction in tumors with high infiltration of CTLs and the prevention of T-cell infiltration in tumors with low CTL levels (43). The TIDE score has been proven to predict the outcome of melanoma patients treated with first-line anti-PD1 or anti-CTLA4 antibodies more accurately than other biomarkers, such as PD-L1 level and mutation load (43). Besides, TIS developed by NanoString Technologies as a clinical-grade assay provides both quantitative and qualitative information about the TME is an 18-gene signature including genes that reflect an ongoing adaptive Th1 and cytotoxic CD8 T-cell response has shown promising results in predicting response to anti-PD-1/PD-L1 agents and has been validated in two HNSCC cohorts from prospective clinical trials (KEYNOTE-012 and KEYNOTE-055) treated with single-agent pembrolizumab showing a positive correlation with response and survival (47, 48). However, both of TIDE and TIS focused on the function and status of T cells, which could not fully reflect the complexity of the TME involved in the response to immunotherapy. Also, both of the two markers focused on patient response to immunotherapy rather than patient survival time, and life expectancy was also important in making treatment decisions. In our study, the predictive value of IRGPI was comparable with TIDE and TIS, and IRGPI might be a better predictor of OS at longer follow-up. Moreover, IRGPI was composed of only three genes and so it was easier to detect than TIDE and TIS.
In conclusion, IRGPI is a promising immune-related prognostic biomarker. IRGPI grouping may help in distinguishing immune and molecular characteristics, predicting patient outcomes. IRGPI might be a potential prognostic indicator of immunotherapy, but further studies are needed to clarify this point (Fig. 6).
No disclosures were reported.
Y. Chen: Conceptualization, software, visualization, writing-original draft, writing-review and editing, collection and assembly of data. Z.-Y. Li: Writing-original draft, writing-review and editing, collection and assembly of data, data analysis and interpretation. G.-Q. Zhou: Conceptualization, writing-review and editing. Y. Sun: Conceptualization, supervision, funding acquisition, project administration.
This work was supported by grants from the Health & Medical Collaborative Innovation Project of Guangzhou City, China (no. 201604020003, 201803040003), the Special Support Program of Sun Yat-sen University Cancer Center (16zxtzlc06), the Natural Science Foundation of Guangdong Province (no. 2017A030312003), the Innovation Team Development Plan of the Ministry of Education (no. IRT_17R110), and the National Key R&D Program of China (2016YFC0902000).
The authors are grateful to Vonn Andrew Walter (Lineberger Comprehensive Cancer Center, University of North Carolina, Chapel Hill, NC) and Pierre Saintigny (Department of Translational Research and Innovation, Univ Lyon, Université Claude Bernard Lyon 1, France) for providing the molecular subtypes of the 502 HNSCC samples. The authors thank Yu-Pei Chen (Sun Yat-sen University Cancer Center, Guangzhou, China) for providing the immune subtypes of the 502 HNSCC samples. Finally, the authors thank Ze-Xian Liu (Sun Yat-sen University Cancer Center, Guangzhou, China) for helping them to check the literature and analyze the data during revision, and thank Li Lin (Sun Yat-sen University Cancer Center, Guangzhou, China) for helping them to check the literature and revise the article during revision.