- Research article
- Open Access
- Open Peer Review
A 13-gene expression-based radioresistance score highlights the heterogeneity in the response to radiation therapy across HPV-negative HNSCC molecular subtypes
BMC Medicine volume 15, Article number: 165 (2017)
Radiotherapy for head and neck squamous cell carcinomas (HNSCC) is associated with a substantial morbidity and inconsistent efficacy. Human papillomavirus (HPV)-positive status is recognized as a marker of increased radiosensitivity. Our goal was to identify molecular markers associated with benefit to radiotherapy in patients with HPV-negative disease.
Gene expression profiles from public repositories were downloaded for data mining. Training sets included 421 HPV-negative HNSCC tumors from The Cancer Genome Atlas (TCGA) and 32 HNSCC cell lines with available radiosensitivity data (GSE79368). A radioresistance (RadR) score was computed using the single sample Gene Set Enrichment Analysis tool. The validation sets included two panels of cell lines (NCI-60 and GSE21644) and HPV-negative HNSCC tumor datasets, including 44 (GSE6631), 82 (GSE39366), and 179 (GSE65858) patients, respectively. We finally performed an integrated analysis of the RadR score with known recurrent genomic alterations in HNSCC, patterns of protein expression, biological hallmarks, and patterns of drug sensitivity using TCGA and the E-MTAB-3610 dataset (659 pancancer cell lines, 140 drugs).
We identified 13 genes differentially expressed between tumor and normal head and neck mucosa that were associated with radioresistance in vitro and in patients. The 13-gene expression-based RadR score was associated with recurrence in patients treated with surgery and adjuvant radiotherapy but not with surgery alone. It was significantly different among different molecular subtypes of HPV-negative HNSCC and was significantly lower in the “atypical” molecular subtype. An integrated analysis of RadR score with genomic alterations, protein expression, biological hallmarks and patterns of drug sensitivity showed a significant association with CCND1 amplification, fibronectin expression, seven hallmarks (including epithelial-to-mesenchymal transition and unfolded protein response), and increased sensitivity to elesclomol, an HSP90 inhibitor.
Our study highlights the clinical relevance of the molecular classification of HNSCC and the RadR score to refine radiation strategies in HPV-negative disease.
Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer worldwide and the second smoking-related cancer after lung carcinoma . Management of HNSCC is challenging due to the high morbidity and mortality of surgery and radiotherapy . Although intensity-modulated radiotherapy has greatly improved the quality of life of patients , especially by decreasing the incidence of early acute and late grade xerostomia , the identification of biomarkers that may help tailor radiotherapy is needed.
Several approaches have been used to predict the tumor response to radiotherapy, including functional assays (clonogenic cell survival  or DNA repair ) or omics-based biomarkers [3, 7,8,9,10]. The latter have been identified in cancer cell lines and their expression in normal cells has not been evaluated. Tumor samples from patients harbor a diverse range of stromal/normal cells that alter the tumor purity, which is known to influence the biological interpretation of results generated from genomic approaches [11,12,13]. Therefore, the potential impact of normal tissue contaminating the tumor sample on intrinsic radioresistance has to be considered.
Oropharyngeal human papillomavirus (HPV)-positive HNSCC are associated with an increased radiosensitivity and a favorable outcome compared to HPV-negative HNSCC [14,15,16]. De-escalation strategies are currently being evaluated in clinical trials in this setting [17, 18]. The higher radiosensitivity of HPV-positive HNSCC may be due to a more effective cytotoxic T-cell-related antitumor immune response , as well as specific molecular alterations, such as p16 overexpression, that decrease DNA repair by inhibiting the recruitment of RAD51 to the site of DNA damage . In HPV-negative HNSCC, the epithelial-to-mesenchymal transition (EMT) via overexpression of fibronectin 1 (FN1), as well as activation of the epidermal growth factor receptor (EGFR) have been implicated in HNSCC radioresistance [9, 21,22,23]. Genome-wide expression profiling of a large number of HNSCC led to the identification of four robust molecular classes of HNSCC [24,25,26]. In this classification, the “classical”, “basal” and “mesenchymal” subtypes exhibit canonical genomic alterations such as focal EGFR amplification, high frequency of HRAS mutations, and upregulation of EMT-related genes, respectively [25, 26]. Therefore, it is tempting to hypothesize that these HNSCC classes are more radioresistant compared to the “atypical” class characterized by the lack of EGFR amplification and enriched in HPV-positive tumors. Of note, a significant proportion of HNSCC in the “atypical” subtype are HPV-negative and the radiosensitivity of these tumors is unknown.
In this work, we made use of gene expression profiles available in public repositories to define a 13-gene expression-based radioresistance (RadR) score in HPV-negative HNSCC and evaluate its association with the current molecular classification, known recurrent genomic alterations in HNSCC, protein expression, biological hallmarks, and patterns of drug sensitivity in vitro. We found the RadR score to be (1) associated with poor disease free-survival (DFS) in HPV-negative HNSCC patients treated by surgery and radiotherapy (with or without chemotherapy) but not in patients treated with surgery alone; (2) lower in the HPV-negative atypical molecular subtype and higher in the mesenchymal subtype; and (3) significantly associated with CCND1 amplification, fibronectin expression and seven biological hallmarks, including the ENT and the actionable “unfolded protein response” hallmarks. Current molecular classification of HNSCC and the RadR score may help to tailor radiotherapy in HPV-negative HNSCC.
No ethical approval was needed for this study.
Cancer cell line datasets
We identified genome-wide expression profiles from four sets of established cancer cell lines with available radiation or drug sensitivity data. The first set included 59 cell lines established from 9 different cancer tissue types (breast, central nervous system, colon, leukemia, melanoma, non-small cell lung, ovarian, prostate, and renal) and included in the NCI-60 panel. Raw data was downloaded from the CellMiner database [27, 28] and information on the surviving cell fraction at 2 Gy (SF2) was retrieved . Radiosensitive and radioresistant cancer cell lines were defined by a SF2 < 0.2 and > 0.8, respectively, as previously described . Two or three replicates were available for each cell line. The second and third sets included 32 and 5 HNSCC cell lines. Raw data were downloaded from GSE79368  and GSE21644 , respectively. In these two sets, intrinsic radioresistance was assessed by a 96-well plate clonogenic assay and the cutoff was set at a median area under the curve of 2.0 in order to identify radioresistant versus radiosensitive cell lines.
Finally, a set of 659 established cancer cell lines from 13 different cancer tissue types (nervous system, soft tissue, digestive system, bone, lung, urogenital system, aero-digestive tract, skin, kidney, breast, pancreas, blood and thyroid) was included in the analysis. Raw data was downloaded from ArrayExpress (E-MTAB-3610) . The corresponding half maximal inhibitory concentration (IC50) for 140 drugs was retrieved from the Genomics of Drug Sensitivity in Cancer data portal .
Detailed information (type and number of samples, platform, and normalization) is provided in Additional file 1: Table S1.
Primary HNSCC datasets
For the ‘The Cancer Genome Atlas’ (TCGA) dataset, RNAseqV2 normalized data (level 3, reads per kilobase per million mapped reads) of 518 primary HNSCC (421 HPV-negative and 97 HPV-positive) and 43 paired normal mucosa surrounding tumors was downloaded using the ‘TCGA2STAT’ R package . From a previous publication of the comprehensive characterization of HNSCC , as well as from the Broad Institute TCGA GDAC Firehose using a q-value < 0.05 [34, 35], 73 genes recurrently altered at the genomic level in HNSCC were selected. Mutation status and copy number alterations from GISTIC  of these genes were downloaded using cBioPortal [37, 38]. Expression of 237 proteins and phosphoproteins generated by Reverse-Phase Protein Array, corresponding to 295 of the 421 primary HPV-negative HNSCC , was downloaded from The Cancer Proteome Atlas. For these 518 primary HNSCC, we queried the TCGA data portal as well as the cBioportal [37, 38] to retrieve clinical data, HPV status, the percentage of copy number alteration and the mutation count per sample. The molecular subtype of the 518 HNSCC was kindly provided by Vonn Andrew Walter (Lineberger Comprehensive Cancer Center, University of North Carolina at Chapel Hill, Chapel Hill, NC, USA). The percentage of tumor cells (tumor purity) was given for each tumor sample by the consensus purity estimation method and was retrieved from a previous publication .
Other datasets included 44 paired HNSCC and normal mucosa (GSE6631) , and two cohorts of 252 and 138 primary HNSCC retrieved from GSE65858  and GSE39366 , respectively. Molecular subtypes, HPV status, and clinical information were available for tumors included in GSE39366 and GSE65858.
A detailed description of these datasets is available in Additional file 1: Table S1.
Gene selection and computation of a RadR score
The following sequential approach was used to select genes associated with radioresistance in HPV-negative HNSCC. (1) In order to overcome the confounding effect of tumor purity [11,12,13], genes had to be differentially expressed in tumor versus normal mucosa; (2) among them, genes had to be differentially expressed in radioresistant versus radiosensitive HNSCC cell lines; and (3) finally, genes had to be associated with DFS in patients with HPV-negative HNSCC and treated by surgery and adjuvant radiation. Data from TCGA and GSE79368 were used as training sets for gene selection.
The single sample gene set enrichment analysis tool (ssGSEA) [42, 43] was used to compute the enrichment score (ES) using the selected genes, in the training sets (TCGA and GSE79368) and in the validation sets (GSE6631, GSE21644, NCI-60, and GSE39366). The ES of these selected genes was referred to as the RadR score. The RadR score was compared in normal versus tumor samples (GSE6631) and in radiosensitive versus radioresistant cancer cell lines (GSE21644 and NCI-60 panel), and its association with DFS was tested in GSE39366, which included patients who underwent surgical resection alone or followed by adjuvant radiotherapy alone or in combination with chemotherapy.
Integrated analysis of the RadR score with genomic alterations, biological hallmarks, and patterns of protein expression and drug sensitivity
In the TCGA dataset, we performed an integrated analysis of the RadR score with known recurrent genomic alterations in HNSCC, genomic instability assessed by the percentage of copy number alterations (CNAs) and the mutation count provided by the cBioportal [37, 38]. Using the ssGSEA tool, we then calculated the Pearson’s coefficient of correlation of the RadR score with the ES of a collection of 48 biological hallmark gene sets previously reported  in three independent datasets, namely TCGA, GSE39366, and GSE65858. We then evaluated the correlation of the RadR score with the expression pattern of a panel of 237 proteins and phosphoproteins in 295 primary HPV-negative HNSCC from TCGA, as well as with expression of 188 proteins and phosphoproteins in 59 cancer cell lines from NCI-60. Finally, the RadR score was computed in 659 established cancer cell lines from various tumor types (E-MTAB-3610) and its correlation with the IC50 for 140 drugs was tested.
Bioinformatics and statistical analysis
Bioinformatics and statistical analysis were performed using the Array Studio software (Omicsoft Corporation), Bioconductor packages in the R language (http://www.bioconductor.org)  and GraphPad Prism version 6.00 (San Diego, SA). Detailed information on processing and normalization of raw data are available in Additional file 1: Table S1. The ssGSEA tool from the GenePattern public server was used to compute the RadR score and the hallmarks’ ESs in various datasets . The unpaired non-parametrical Mann–Whitney or Wilcoxon test was used to compare continuous variables between two groups and the Kruskal–Wallis test was used when more than two groups were compared. The Pearson correlation coefficient (r) was estimated to measure the strength of a linear association between two continuous variables. In the GSE39366 study, recurrence-free survival time was defined by the time in months from tumor biopsy to death, recurrence, or loss to follow-up. In the TCGA cohort, DFS was defined by the time in months from tumor biopsy to a new tumor event. Survival distributions were estimated using the Kaplan–Meier method and compared with the log‐rank test between groups of patients defined by the median RadR score (low versus high RadR score). Multivariate cox proportional hazard model, including age, node, and tumor stage, were built to select genes associated with DFS in patients treated by surgery and radiation. All statistical tests were two-sided, and P values of more than 0.05 were considered to be statistically significant. For multiple testing, a false-discovery rate was computed using the Bonferroni–Holm method in order to adjust P value.
Identification of radioresistance-associated genes in HPV-negative HNSCC
A sequential approach was used to select genes associated with radioresistance, by selecting (1) genes differentially expressed in tumor versus normal head and neck samples; (2) genes associated with in vitro radioresistance, and (3) genes associated with a poor DFS in patients treated by surgery and radiotherapy (Fig. 1).
In order to overcome the confounding effect of tumor purity on intrinsic radioresistance , we first filtered out genes that were not differentially expressed in cancer versus normal tissue using gene expression profiles of 86 paired HNSCC and normal mucosa from TCGA (FDR < 0.05, |Fold Change (FC)| > 1, paired Wilcoxon test) (Fig. 1). Among the genes differentially expressed in tumor versus normal tissue, we selected those associated with radioresistance in vitro in 32 HNSCC cell lines (18 and 14 radioresistant and radiosensitive cell lines, respectively) analyzed at baseline, 2 and 6 h after 4 Gy radiation (GSE79368) . Using a two way-ANOVA independently of the time point (with a FDR < 0.05 and a |FC| > 1) (Fig. 1), all samples were included to increase the statistical power, allowing to identify 75 genes over-expressed and 146 genes under-expressed in radioresistant cell lines compared to radiosensitive ones (Fig. 1). These genes were also over- and under-expressed, respectively, in tumor compared to normal mucosa of the head and neck (Fig. 1). Among these 221 genes, we identified a set of clinically relevant genes by their association with DFS in 128 of the 421 patients with HPV-negative HNSCC from TCGA who were treated by surgery and adjuvant radiotherapy with or without concurrent chemotherapy, using a multicovariate Cox proportional hazard model, including age, node (pN0-N1 vs. pN2-N3) and tumor stage (pT1-T2 vs. pT3-T4) (Fig. 1). Clinical and pathological parameters of these 128 patients from TCGA are detailed in Additional file 2: Table S2. Median follow-up was 20.5 months. This led to the identification of 3 and 10 genes associated with increased (hazard ratio (HR) > 1, P < 0.05) and decreased risk (HR < 1, P < 0.05) of recurrence, respectively. The gene list is detailed in Table 1. Detailed information of the FC at each time points and P values using a Wilcoxon test for the 13 selected genes in radiosensitive and radioresistant HNSCC cell lines from GSE79368 are provided in Additional file 3: Table S3.
Validation of the RadR score in HPV-negative HNSCC
Using the ssGSEA tool, the 13-gene expression-based RadR score was computed in multiple independent datasets available from public repositories, namely GSE6631, GSE21644, GSE3966, and the NCI-60 panel. The RadR score was increased in 19/22 HNSCC tumors when compared to paired normal mucosa in GSE6631. Overall, the score was significantly increased when tumor and paired normal mucosa were compared (Wilcoxon paired test P < 0.0001) (Fig. 2a). Moreover, the RadR score was higher in radioresistant cell lines (UTSCC-77, UTSCC-24A) compared to the radiosensitive UTSCC-9 cell line (GSE21644) (Fig. 2b). Of note, among the five HNSCC cell lines that were profiled in duplicate by independent investigators (GSE21644), four were overlapping with the cell lines used in the training set (GSE79368). Thus, the RadR score was also computed in a larger panel of 59 cancer cell lines from various tumor types (NCI-60), which were profiled in two or three replicates (174 samples). Interestingly, the RadR score computed in each replicate was highly correlated to other replicates for the same cell line (r ≥ 0.9, P < 0.0001, see Additional file 4: Table S4). Using the first replicate, the RadR score was significantly different between the three groups of NCI-60 cell lines as defined by different ranges of survival fraction at 2 Gy (SF2) as shown in Fig. 2c (P = 0.0236, Kruskal–Wallis test). Radioresistant lines (SF2 > 0.8) had higher RadR scores compared to radiosensitive lines (SF2 < 0.2) (P = 0.0162, Mann–Whitney test).
The clinical relevance of the RadR score was assessed in a set of 63 primary HNSCC using GSE39366. In this study, HPV-positive tumors were excluded and patients were divided into two groups according to treatment. A high RadR score (i.e., greater than the median score) was significantly associated with a shorter recurrence-free survival (log-rank P = 0.0384) in 34 patients treated with surgery and adjuvant radiotherapy with or without concurrent chemotherapy (Fig. 2d). Interestingly, no difference was observed in 29 patients treated with surgery alone or with chemotherapy (P = 0.5218). This finding underlines the specificity of the RadR score and validates our strategy to include the selection of genes associated with in vitro radioresistance and poor DFS in patients surgically resected and treated with postoperative radiation therapy. When validated in independent cohorts, the RadR score may therefore help the clinician to tailor radiation dose according to the radioresistance level of a tumor.
In summary, these results validate the association of the RadR score with in vitro radioresistance and with the benefit from radiotherapy in patients treated with surgery and adjuvant radiotherapy but not in patients treated with surgery alone.
The RadR score is associated with the molecular classification of HNSCC
In large genomic profiling studies of HNSCC, four distinct molecular subtypes have been consistently reported, namely atypical, basal, classical, and mesenchymal [25, 26]. HPV-positive HNSCC, which arise from the oropharynx and often fall into the “atypical” class, are recognized to be more radiosensitive compared to HPV-negative tumors [14, 15]. In line with these reports, HPV-positive HNSCC showed lower RadR scores compared to HPV-negative tumors (P < 0.01) in three independent sets of HNSCC (TCGA, GSE65858 and GSE39366) (Additional file 5: Figure S1).
We then focused on the RadR score distribution according to the molecular subtypes and the anatomical subsites in HPV-negative HNSCCs from these three independent datasets (TCGA, n = 421; GSE65858, n = 179; GSE39366, n = 82). The RadR score was found to be statistically different across molecular subtypes, and consistently lower in atypical compared to classical, basal, and mesenchymal HPV-negative HNSCC (P < 0.0001 in all three datasets) (Fig. 3a). Furthermore, the RadR score was different according to the anatomical subsite with a trend to be higher in the oral cavity in the TCGA (P < 0.0001) and GSE39366 (P = 0.0002) datasets, although no difference was observed in GSE65858 (P = 0.1243) (Fig. 3b).
Overall, our results provide evidence that, among HPV-negative HNSCC, “atypical” HNSCCs may be more sensitive to radiotherapy when compared to other subtypes, whereas “mesenchymal” HNSCCs may be more radioresistant. In contrast, the association of the anatomical subsite with the RadR score was not consistent across the three sets of data analyzed.
Molecular subtypes of HPV-negative HNSCC, recurrent genomic alterations and RadR score
Since gene expression-based molecular subtypes have been shown to be enriched with some genomic alterations, such as EGFR amplification and HRAS-CASP8 co-mutations in the classical and basal subtypes, respectively [25, 26], we then tested the association of the RadR score with genomic instability and known recurrent genomics alterations in HNSCC.
No clear association was observed between the RadR score and genomic instability assessed by the percent of CNAs (r = 0.012; P = 0.803) and the mutation count (r = –0.204, P = 0.002). We then further investigated the association between the RadR score with somatic mutations and CNAs of 73 genes recurrently altered in HNSCC as reported previously , and defined by the Broad Institute TCGA GDAC Firehose [34, 35]. Somatic mutations as well as putative CNAs (GISTIC) of these genes were retrieved from cBioPortal and used to perform an integrated analysis of RadR score with overall genomic instability, mutation status and CNAs in 220 HPV-negative HNSCC (Fig. 4a). A decreased RadR score was significantly associated with NSD1 (P < 0.0001) and PIK3CA (P = 0.0042) mutations, and to a lesser degree with TP63 amplification (P = 0.0108); while HNSCC harboring CCND1 amplification had a higher RadR score (P = 0.0354) (Fig. 4b). Noticeably, percentages of NSD1 and PIK3CA mutations as well as TP63 amplification were significantly higher in the atypical and classical subtypes between HPV-negative HNSCC (P = 0.0013, P = 0.0003, and P = 0.0002, respectively), whereas the percentage of CCND1 amplification was not different between the four subtypes (P = 0.7130) (Fig. 4c).
The influence of tumor purity on genomic analysis has been well established [11,12,13]. We thus evaluated the association between the RadR score and the tumor purity given by the consensus purity estimation method . No significant association was found in the 421 HPV-negative HNSCC from TCGA (r = –0.088, P = 0.0707) (Additional file 6: Figure S2), suggesting that the RadR score is independent of the tumor purity.
In summary, the RadR score was associated with some specific known recurrent genomic alterations that are known to be present at different rates in the molecular subtypes of HPV-negative HNSCC, suggesting that genomic-driven radioresistance may be different among different molecular subtypes. The lack of correlation between the RadR score and the tumor purity greatly enhances our confidence in this result.
Molecular subtypes of HPV-negative HNSCC, biological hallmark, patterns of protein expression and the RadR score
Since few actionable associations were observed with known recurrent genomic alterations in HNSCC, we seek to gain more insights into the biological context associated with radioresistance, using an in silico functional approach based on the integrative analysis of biological pathways and proteomics data provided by The Cancer Proteome Atlas .
We computed the Pearson’s correlation of the RadR score with the ESs of 48 gene sets characterizing biological hallmarks in three independent datasets of primary HPV-negative HNSCC (TCGA, n = 421; GSE65858, n = 179; GSE39366, n = 82). In this analysis, seven hallmarks were significantly correlated with the RadR score in the three datasets (P < 0.01, Fig. 5a): “MYC_targets v2” (average r = 0.30), “E2F targets” (average r = 0.29), “unfolded protein response” (UPR gene set; average r = 0.30), “TGFb signaling” (average r = 0.38), “angiogenesis” (average r = 0.31), “DNA repair” (average r = 0.31), and EMT (average r = 0.38) (Fig. 5a). The ESs of these seven radioresistance-associated hallmarks were significantly different between the molecular subtypes (Kruskal–Wallis test, FDR < 0.0001). Interestingly, the mesenchymal subtype had higher ES for the hallmarks “EMT” and “angiogenesis”, whereas the classical subtype was associated with a higher ES for the hallmark “E2F targets” and “MYC targets”. Notably, the ES of the “UPR” hallmark were lower in the atypical tumors compared to others.
We then analyzed the expression of 237 proteins and phosphoproteins in 295 HPV-negative HNSCC from TCGA and tested their correlation with the RadR score. Interestingly, the most positively correlated proteins were fibronectin and PAI1 (r = 0.3486, P < 0.0001; r = 0.3507, P < 0.0001, respectively), while E-cadherin was negatively correlated (r = –0.3016, P < 0.0001), consistent with the association of the RadR score with EMT. Using a coefficient threshold of 0.3, we found no relevant correlation of phosphoproteins with the RadR score.
In order to assess the association of proteins and phospho-proteins with in vitro radioresistance, we looked for proteins that were correlated to the RadR score as well as to SF2 in 59 cancer cell lines from NCI60 for consistency, using a Pearson’s correlation. A significant correlation was found for Akt (vs. RadR score: r = 0.5351, P < 0.0001; vs. SF2: r = 0.2836, P = 0.0295), fibronectin (vs. RadR score: r = 0.3687, P = 0.0041; vs. SF2: r = 0.259, P = 0.0476), cyclin D1 (vs. RadR score: r = 0.3587, P = 0.0053; vs. SF2: r = 0.4027, P = 0.0016), Irs1 (vs. RadR score: r = 0.2921, P = 0.0248; vs. SF2: r = 0.2885, P = 0.0267), Rictor pT1135 (vs. RadR score: r = 0.2742, P = 0.0356; vs. SF2: r = 0.3047, P = 0.0189), and Paxillin (vs. RadR score: r = 0.2641, P = 0.0433; vs. SF2: r = 0.452, P = 0.0003). Of note, cyclin D1 and fibronectin protein levels were increased according to the range of SF2 in cancer cell lines from NCI-60 (P = 0.0015 and P = 0.0325, respectively) (Fig. 6a). Moreover, fibronectin was significantly higher in the mesenchymal subtype of HPV-negative HNSCC, whereas no difference was observed for cyclin D1 protein expression between the four subtypes.
The association of specific biological hallmarks and proteomic alterations with the RadR score and molecular subtypes of HPV-negative HNSCC suggests that drugs targeting these hallmarks may be relevant to treat tumors with a high RadR score, according to their molecular classification. We computed the RadR score in a large panel of 659 cancer cell lines with available IC50 (half maximal inhibitory concentration) for 140 drugs . By correlating the RadR score with the IC50, four drugs were positively correlated to the RadR score (r > 0.2; P < 0.01) (Additional file 4: Figure S3). They included cisplatin, a potent radiosensitizer routinely used in patients concurrently with radiation therapy, the pro-apoptotic Bcl-2 inhibitor TW37, the BRAF inhibitor PLX4720, and elesclomol, an inhibitor of HSP90 that modulates the UPR pathway [47, 48]. Of note, the RadR score was significantly higher in cell lines sensitive to elesclomol compared to resistant cell lines as previously defined  (Additional file 7: Figure S3), suggesting that targeting the radioresistance-associated UPR hallmark may be relevant in patients with a high RadR score.
Our results evidence an important role of seven specific hallmarks into the radioresistance of HPV-negative HNSCC, particularly the EMT hallmark, which is in line with an increased radioresistance observed in the mesenchymal subtype, and the UPR hallmark, which may be modulated by elesclomol in HPV-negative HNSCC with a high RadR score. The association of these hallmarks with molecular subtypes of HPV-negative HNSCC may lead to the generation of hypotheses to refine combination therapies of radiation with modulators of EMT-driven cancer cell plasticity or of the UPR pathway.
Radiotherapy is indicated and may benefit approximately 75% of patients with HNSCC , demonstrating its key role in the management of these patients. However, radiotherapy is challenging in part due to the complex anatomy of the upper aerodigestive tract, which makes it difficult to safely deliver an efficient dose to the tumor and lymph nodes when indicated. It is now well established that HPV-positive tumors are more radiosensitive [14, 15]. In this study, we focus on HPV-negative HNSCC and identify a 13-gene expression-based RadR score predictive of the benefit of radiotherapy in vitro and in patients. We show that molecular subtypes of HPV-negative HNSCC were characterized by different levels of radioresistance, where the atypical subtype was the most radiosensitive, whereas the mesenchymal was the most radioresistant, as defined by the RadR score in this population of patients. Our results suggest that the molecular classification of HPV-negative HNSCC and the RadR score may help refine radiation therapy strategies.
Gene-wide expression profiles have been previously used to identify various gene expression signatures or biomarkers associated with radiosensitivity and/or locoregional recurrence in head and neck cancers [7,8,9,10, 50]. However, the score of these signatures and biomarkers are poorly described in normal mucosa surrounding the tumor or tumor-associated normal cells, introducing bias for their interpretation in tumors with heterogeneous purity. Indeed, the influence of tumor purity on genomic analysis is well established [11,12,13]. In order to overcome this potential bias, we filtered out genes that were not differentially expressed between normal and tumor samples and showed that the RadR score was significantly higher in tumors compared to normal samples and was not correlated to tumor purity in the HNSCC TCGA dataset. Then, we selected genes that were differentially expressed between radioresistant and radiosensitive head and neck cancer cell lines and those which were associated with DFS in HPV-negative HNSCC from TCGA. However, because clinical information in TCGA has not been collected prospectively and monitored as in a clinical trial, the association of the RadR score with DFS needs to be validated in prospective cohorts.
Previous studies have shown a higher radiosensitivity in HPV-positive oropharyngeal SC, leading to de-escalation treatment protocols to decrease morbidity [14, 15, 17]. Consistently, the RadR scores were found to be lower in HPV-positive versus HPV-negative tumors. Moreover, the vast majority of HPV-positive tumors fall into the “atypical” molecular class in the most recent and extensive classification of HNSCC by the TCGA . Intriguingly, HPV-negative tumors falling into the “atypical” molecular class had lower RadR scores in three independent datasets and may therefore be more radiosensitive and benefit from similar de-escalation strategies.
Since molecular subtypes of HNSCC exhibit different canonical genomic alterations, such as EGFR amplification or HRAS mutation [25, 26], it was tempting to hypothesize that the association of radioresistance, as defined by the RadR score, with the molecular classification may be related to specific genomic alterations. Interestingly, the PIK3CA mutations, most commonly found in radiosensitive HPV-positive HNSCC , were also more commonly found in the atypical subtype, a rather radiosensitive subtype based on low RadR scores. These observations suggest that similar molecular alterations, such as PIK3CA mutations, may enhance radiosensitivity in atypical tumors irrespective of their HPV status. Similarly, NSD1 mutation and amplification of TP63 were associated with a low RadR score and more frequent in the atypical and classical subtypes.
In order to gain more insight into the biological characterization of radioresistance in HPV-negative HNSCC and because few targetable genomic alterations were positively associated with the RadR score, we correlated the RadR score with the ES of 48 distinct biological hallmarks as well as with the expression of 237 phosphoproteins and proteins. We found seven biological hallmarks that were consistently correlated to the RadR score in the three independent datasets, namely “TGFβ signaling”, “DNA repair”, “angiogenesis”, “UPR”, “E2F targets”, “MYC targets”, and “EMT”. In line with our results, angiogenesis and DNA repair pathways are known to play an important role in radioresistance in HNSCC . Moreover, the “EMT” hallmark was significantly and positively correlated to the RadR score, supporting the previously reported association of EMT with radioresistance in HNSCC [9, 23]. As expected, the mesenchymal subtype had higher ESs for the hallmarks “EMT”, “angiogenesis” and “TGFβ signaling”, as compared to other subtypes; whereas the classical subtype was characterized by higher ESs for the hallmarks “E2F targets” and “MYC targets”, thus suggesting that radioresistance may be driven by different hallmarks depending on the molecular subtype of a specific HPV-negative HNSCC.
Furthermore, we found that the RadR score was associated with patterns of drug sensitivity among the 140 drugs available for analysis in a large panel of established cancer cell lines [31, 32]. Computing the correlation between IC50 of these drugs and the RadR score, we found that four drugs may be active in cancer cell lines harboring high RadR scores, including the well-known potent radiosensitizer cisplatin. Interestingly, the IC50 of elesclomol, an inhibitor of the HSP90 chaperone protein which modulates the UPR pathway [47, 48, 51], was also negatively correlated to the RadR score, consistent with the likely important role of this pathway in HPV-negative HNSCC radioresistance. This will need to be validated in preclinical models.
We report a 13-gene expression-based RadR score that is associated with specific molecular features of HPV-negative HNSCC. To the best of our knowledge, our study is the first report suggesting the clinical relevance of the molecular classification of HNSCC in order to refine radiation strategies. The predictive value of the RadR score needs to be validated in larger retrospective and prospective cohorts. Our work is hypothesis-generating for the refinement of combination therapies of radiation with modulators of EMT-driven cancer cell plasticity or of the UPR pathway.
copy number alterations
enrichment score (computed using single sample Gene Set Enrichment Analysis)
false discovery rate
head and neck squamous cell carcinomas
- IC50 :
half maximal inhibitory concentration (a measure of drug sensitivity)
survival fraction at 2 Gy (a measure of radiosensitivity)
single sample gene set enrichment analysis (for computing an enrichment score of a gene set in a unique sample)
The Cancer Genome Atlas (Public repository)
unfolded protein response
Parkin DM, Bray F, Ferlay J, Pisani P. Global cancer statistics, 2002. CA Cancer J Clin. 2005;55(2):74–108.
Baxi SS, Sher DJ, Pfister DG. Value considerations in the treatment of head and neck cancer: radiation, chemotherapy, and supportive care. Am Soc Clin Oncol Educ Book. 2014:e296–303. doi:10.14694/EdBook_AM.2014.34.e296.
Singer S, Langendijk J, Yarom N. Assessing and improving quality of life in patients with head and neck cancer. Am Soc Clin Oncol Educ Book. 2013. doi:10.1200/EdBook_AM.2013.33.e230.
Ghosh-Laskar S, Yathiraj PH, Dutta D, Rangarajan V, Purandare N, Gupta T, Budrukkar A, Murthy V, Kannan S, Agarwal JP. Prospective randomized controlled trial to compare 3-dimensional conformal radiotherapy to intensity-modulated radiotherapy in head and neck squamous cell carcinoma: Long-term results. Head Neck. 2016;38 Suppl 1:E1481–1487.
Fertil B, Malaise EP. Inherent cellular radiosensitivity as a basic concept for human tumor radiotherapy. Int J Radiat Oncol Biol Phys. 1981;7(5):621–9.
Chavaudra N, Bourhis J, Foray N. Quantified relationship between cellular radiosensitivity, DNA repair defects and chromatin relaxation: a study of 19 human tumour cell lines from different origin. Radiother Oncol. 2004;73(3):373–82.
Skinner HD, Giri U, Yang L, Woo SH, Story MD, Pickering CR, Byers LA, Williams MD, El-Naggar A, Wang J, et al. Proteomic profiling identifies PTK2/FAK as a driver of radioresistance in HPV negative head and neck cancer. Clin Cancer Res. 2016;22(18):4643–50.
Hall JS, Iype R, Senra J, Taylor J, Armenoult L, Oguejiofor K, Li Y, Stratford I, Stern PL, O'Connor MJ, et al. Investigation of radiosensitivity gene signatures in cancer cell lines. PLoS One. 2014;9(1):e86329.
Jerhammar F, Ceder R, Garvin S, Grenman R, Grafstrom RC, Roberg K. Fibronectin 1 is a potential biomarker for radioresistance in head and neck squamous cell carcinoma. Cancer Biol Ther. 2010;10(12):1244–51.
Pramana J, Van den Brekel MW, van Velthuysen ML, Wessels LF, Nuyten DS, Hofland I, Atsma D, Pimentel N, Hoebers FJ, Rasch CR, et al. Gene expression profiling to predict outcome after chemoradiation in head and neck cancer. Int J Radiat Oncol Biol Phys. 2007;69(5):1544–52.
de Ridder D, van der Linden CE, Schonewille T, Dik WA, Reinders MJ, van Dongen JJ, Staal FJ. Purity for clarity: the need for purification of tumor cells in DNA microarray studies. Leukemia. 2005;19(4):618–27.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, Trevino V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Aran D, Sirota M, Butte AJ. Systematic pan-cancer analysis of tumour purity. Nat Commun. 2015;6:8971.
Lindel K, Beer KT, Laissue J, Greiner RH, Aebersold DM. Human papillomavirus positive squamous cell carcinoma of the oropharynx: a radiosensitive subgroup of head and neck carcinoma. Cancer. 2001;92(4):805–13.
Mirghani H, Amen F, Tao Y, Deutsch E, Levy A. Increased radiosensitivity of HPV-positive head and neck cancers: Molecular basis and therapeutic perspectives. Cancer Treat Rev. 2015;41(10):844–52.
Fakhry C, Westra WH, Li S, Cmelak A, Ridge JA, Pinto H, Forastiere A, Gillison ML. Improved survival of patients with human papillomavirus-positive head and neck squamous cell carcinoma in a prospective clinical trial. J Natl Cancer Inst. 2008;100(4):261–9.
Masterson L, Moualed D, Liu ZW, Howard JE, Dwivedi RC, Tysome JR, Benson R, Sterling JC, Sudhoff H, Jani P, et al. De-escalation treatment protocols for human papillomavirus-associated oropharyngeal squamous cell carcinoma: a systematic review and meta-analysis of current clinical trials. Eur J Cancer. 2014;50(15):2636–48.
Mirghani H, Amen F, Blanchard P, Moreau F, Guigay J, Hartl DM, Lacau St Guily J. Treatment de-escalation in HPV-positive oropharyngeal carcinoma: ongoing trials, critical issues and perspectives. Int J Cancer. 2015;136(7):1494–503.
Jung AC, Guihard S, Krugell S, Ledrappier S, Brochot A, Dalstein V, Job S, de Reynies A, Noel G, Wasylyk B, et al. CD8-alpha T-cell infiltration in human papillomavirus-related oropharyngeal carcinoma correlates with improved patient prognosis. Int J Cancer. 2013;132(2):E26–36.
Pawlik TM, Keyomarsi K. Role of cell cycle in mediating sensitivity to radiotherapy. Int J Radiat Oncol Biol Phys. 2004;59(4):928–42.
Jedlinski A, Ansell A, Johansson AC, Roberg K. EGFR status and EGFR ligand expression influence the treatment response of head and neck cancer cell lines. J Oral Pathol Med. 2013;42(1):26–36.
Perri F, Pacelli R, Della Vittoria Scarpati G, Cella L, Giuliano M, Caponigro F, Pepe S. Radioresistance in head and neck squamous cell carcinoma: Biological bases and therapeutic implications. Head Neck. 2015;37(5):763–70.
de Jong MC, Ten Hoeve JJ, Grenman R, Wessels LF, Kerkhoven R, Te Riele H, van den Brekel MW, Verheij M, Begg AC. Pretreatment microRNA expression impacting on epithelial-to-mesenchymal transition predicts intrinsic radiosensitivity in head and neck cancer cell lines and patients. Clin Cancer Res. 2015;21(24):5630–8.
Chung CH, Parker JS, Karaca G, Wu J, Funkhouser WK, Moore D, Butterfoss D, Xiang D, Zanation A, Yin X, et al. Molecular classification of head and neck squamous cell carcinomas using patterns of gene expression. Cancer Cell. 2004;5(5):489–500.
Walter V, Yin X, Wilkerson MD, Cabanski CR, Zhao N, Du Y, Ang MK, Hayward MC, Salazar AH, Hoadley KA, et al. Molecular subtypes in head and neck cancer exhibit distinct patterns of chromosomal gain and loss of canonical cancer genes. PLoS One. 2013;8(2):e56823.
Cancer Genome Atlas Network. Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature. 2015;517(7536):576–82.
Shankavaram UT, Varma S, Kane D, Sunshine M, Chary KK, Reinhold WC, Pommier Y, Weinstein JN. Cell Miner: a relational database and query tool for the NCI-60 cancer cell lines. BMC Genomics. 2009;10:277.
Reinhold WC, Sunshine M, Liu H, Varma S, Kohn KW, Morris J, Doroshow J, Pommier Y. Cell Miner: a web-based suite of genomic and pharmacologic tools to explore transcript and drug patterns in the NCI-60 cell line set. Cancer Res. 2012;72(14):3499–511.
Amundson SA, Do KT, Vinikoor LC, Lee RA, Koch-Paiz CA, Ahn J, Reimers M, Chen Y, Scudiero DA, Weinstein JN, et al. Integrating global gene expression and radiation survival parameters across the 60 cell lines of the National Cancer Institute Anticancer Drug Screen. Cancer Res. 2008;68(2):415–24.
Kim HS, Kim SC, Kim SJ, Park CH, Jeung HC, Kim YB, Ahn JB, Chung HC, Rha SY. Identification of a radiosensitivity signature using integrative metaanalysis of published microarray data for NCI-60 cancer cells. BMC Genomics. 2012;13:348.
Iorio F, Knijnenburg TA, Vis DJ, Bignell GR, Menden MP, Schubert M, Aben N, Goncalves E, Barthorpe S, Lightfoot H, et al. A landscape of pharmacogenomic interactions in cancer. Cell. 2016;166(3):740–54.
Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, Bindal N, Beare D, Smith JA, Thompson IR, et al. Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 2013;41(Database issue):D955–961.
Wan YW, Allen GI, Liu Z. TCGA2STAT: simple TCGA data access for integrated statistical analysis in R. Bioinformatics. 2016;32(6):952–4.
Broad Institute TCGA Genome Data Analysis Center: Mutation Analysis (MutSig v2.0), Broad Institute of MIT and Harvard. 2016. doi:10.7908/C1D21X1G.
Broad Institute TCGA Genome Data Analysis Center: SNP6 Copy number analysis (GISTIC2), Broad Institute of MIT and Harvard. 2016. doi:10.7908/C1V987FP.
Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12(4):R41.
Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012;2(5):401–4.
Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, Sun Y, Jacobsen A, Sinha R, Larsson E, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013;6(269):l1.
Li J, Lu Y, Akbani R, Ju Z, Roebuck PL, Liu W, Yang JY, Broom BM, Verhaak RG, Kane DW, et al. TCPA: a resource for cancer functional proteomics data. Nat Methods. 2013;10(11):1046–7.
Kuriakose MA, Chen WT, He ZM, Sikora AG, Zhang P, Zhang ZY, Qiu WL, Hsu DF, McMunn-Coffran C, Brown SM, et al. Selection and validation of differentially expressed genes in head and neck cancer. Cell Mol Life Sci. 2004;61(11):1372–83.
Wichmann G, Rosolowski M, Krohn K, Kreuz M, Boehm A, Reiche A, Scharrer U, Halama D, Bertolini J, Bauer U, et al. The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer. Int J Cancer. 2015;137(12):2846–57.
Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009;462(7269):108–12.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.
Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25.
Dudoit S, Gentleman RC, Quackenbush J. Open source software for the analysis of microarray data. Biotechniques. 2003;Suppl:45–51.
Reich M, Liefeld T, Gould J, Lerner J, Tamayo P, Mesirov JP. GenePattern 2.0. Nat Genet. 2006;38(5):500–1.
Marcu MG, Doyle M, Bertolotti A, Ron D, Hendershot L, Neckers L. Heat shock protein 90 modulates the unfolded protein response by stabilizing IRE1alpha. Mol Cell Biol. 2002;22(24):8506–13.
Born EJ, Hartman SV, Holstein SA. Targeting HSP90 and monoclonal protein trafficking modulates the unfolded protein response, chaperone regulation and apoptosis in myeloma cells. Blood Cancer J. 2013;3:e167.
Wong K, Delaney GP, Barton MB. Evidence-based optimal number of radiotherapy fractions for cancer: A useful tool to estimate radiotherapy demand. Radiother Oncol. 2016;119(1):145–9.
Eschrich SA, Pramana J, Zhang H, Zhao H, Boulware D, Lee JH, Bloom G, Rocha-Lima C, Kelley S, Calvin DP, et al. A gene expression model of intrinsic tumor radiosensitivity: prediction of response and prognosis after chemoradiation. Int J Radiat Oncol Biol Phys. 2009;75(2):489–96.
Hetz C, Chevet E, Harding HP. Targeting the unfolded protein response in disease. Nat Rev Drug Discov. 2013;12(9):703–19.
We acknowledge Vonn Andrew Walter (Lineberger Comprehensive Cancer Center, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina), who kindly provided us information on the molecular subtype for the 518 tumors from TCGA.
PS: LYric Grant INCa-DGOS-4664; J-PF was supported by fellowship grants “Année recherché” from Assistance Publique Hôpitaux de Paris (2015–2016) and “Soutien pour la formation à la recherche translationnelle en cancérologie” from INCa and AVIESAN (2016–2017).
Availability of data and materials
The datasets analyzed during the current study are available from public repositories: CellMinerTM database version 2.1 (NCI-60), Gene Expression Omnibus (GSE79368; GSE21644; GSE6631; GSE39366; GSE69858), Array Express (E-MTAB-3610), and The Cancer Genome Atlas.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Description of the datasets used in the study. (DOCX 14 kb)
Clinical and pathological parameters of the cohort from TCGA for selection of genes associated with disease-free survival. This cohort includes 128 patients suffering from human papillomavirus (HPV)-negative head and neck squamous cell carcinoma (HNSCC) and who were treated by surgery and adjuvant radiotherapy. (DOCX 16 kb)
Wilcoxon test between radiosensitive and radioresistant head and neck squamous cell carcinoma (HNSCC) cell lines at each time point (baseline: 0 Gy, 2 h after 4 Gy, and 6 h after 4 Gy). FC, fold change; P, P value. (DOCX 15 kb)
Correlation of the RadR score between the different replicates of each cell lines from NCI-60. (DOCX 36 kb)
The radioresistance score in human papillomavirus (HPV)-positive compared to HPV-negative head and neck squamous cell carcinoma (HNSCC) from TCGA, GSE39366, and GSE65858. P value is shown for each dataset. (PPTX 189 kb)
Correlation between tumor purity and RadR score in the 421 primaryHPV-negative HNSCC from TCGA. Pearson’s coefficient of correlation as well as P-value are shown. (PPTX 114 kb)
RadR score and in vitro sensitivity to Elesclomol, TW37, Cisplatine andPLX4720 in >700 cancer cell lines from the Genomics of Drug Sensitivity in Cancer database. TheRadR score was compared between sensitive and resistant cell lines to Cisplatin, TW47, PLX4720 andelesclomol. Number of tested samples is shown for each drug. (PPTX 494 kb)
About this article
Cite this article
Foy, J., Bazire, L., Ortiz-Cuaran, S. et al. A 13-gene expression-based radioresistance score highlights the heterogeneity in the response to radiation therapy across HPV-negative HNSCC molecular subtypes. BMC Med 15, 165 (2017) doi:10.1186/s12916-017-0929-y
- Head neck squamous cell carcinomas
- Molecular subtypes
- Predictive biomarker
- Radiation therapy