Information for the HNSCC training set was obtained from TCGA on November 13, 2016 (18). Gene expression data were extracted from IlluminaHiSeq_RNASeqV2 platform and normalized by RSEM method (19). In addition, we performed quality control with a total of 20530 genes. Genes with more than half of values as zero were removed, 17711 genes remained with quantile normalization. Patients with complete follow-up information and gene expression values for tumor tissues were included in the study. Information for the HNSCC testing set was collected from GSE65858 (20) in GEO. Gene expression data were extracted from Illumina HumanHT-12 V4.0 expression beadchip and normalized using the robust spline normalization (RSN) method (21). Consecutive patients with primary and metachronous secondary HNSCC of oral cavity, larynx, oro- and hypopharynx were included, while tumor cell lines and those with low quality assays were excluded. All gene expression values were log2-transformed and standardized for comparability between the training and testing sets.

To select differentially expressed genes combined with clinical information, WTT was used to select genes based on the method of Hu et al (22). For the ith subject with a covariate vector Z_i, the Cox proportional hazards model is given by λ(t|Z_i) = λ₀(t)exp(β^TZ_i) and the survival function is S(t|Z_i) = exp{-λ₀(t)exp(β^TZ_i)}, where λ₀(t) is the basic hazard function, β is the regression coefficient and λ₀(t) is the cumulative baseline hazard function. Then we constructed a Cox regression model for each subject based on clinical information (age, sex, smoking status and clinical stage) only and defined h_i = β^TZ_i. The weights for n patients totally were calculated accordingly:

which were assigned for the tumor cases but not the normal cases.

With the weighted tumor expression exp_wi = w_i x expi, Students t-test was conducted for each gene to measure the difference between tumor and matched normal expression level. We also used t-test with no weight adjustment and examined the difference between the t-test statistics after and before weight adjustment, d_k = t_adjust - t_unadjust for the kth gene. Afterwards, 1,000 total permutations were performed and d_ki could be got for the ith permutation. Then, we calculated the averaged order statistics, d¯k, across all 1,000 permutations. A gene was labeled as significant when |dk-d¯k| was at the top 5%.

After the WTT selection, there were still over 800 genes left, which were too many and not robust to build the prognostic sigature in HNSCC. The traditional univariate or multivariate Cox regression was not suitable to select the prognosis-associated genes because it easily led to overfitting and produced instable results (23). SIS was used to choose those which were truly associated with disease from the 5% genes remaining for further modeling (24). This is a two-step screening approach: it first screened all genomic features and discarded the irrelevant features whose correlation with overall survival were weak, and secondly applied LASSO penalized regression to estimate the sensitivity from the selected genomic instability data. We could significantly reduce the number of genes in the final model by the SIS method.

Continuous variables are described as mean ± SD, and categorized variables are summarized by frequency (n) and proportion (%). Chi-square test was used for rate or proportion comparison. Associations between the characteristics and the overall survival were evaluated by Cox proportional hazard models. Survival curves were drawn with the Kaplan-Meier method and were compared among subgroups using log-rank tests. To evaluate the robustness of the results, we used the bootstrap method with ‘bootcov’ function that computed a bootstrap estimate of the covariance matrix for a set of regression coefficients in rms package. The bootstrap procedure were carried out with 500 re-samplings for the multivariable Cox regression. We predicted 5-year patient survival using the nearest neighbor method for receiver operating characteristic (ROC) curves of censored survival data (25) and estimation of confidence intervals and P-values of area under the curve (AUC) was based on bootstrap resampling. In the subgroup analysis, we used the Fishers exact test to compare the proportions of different HPV status or tumor sites.

Statistical analyses were performed using R version 3.3.1 (The R Foundation). P-values are two-sided and P<0.05 indicates statistical significance.

The analysis included 512 HNSCC cases from TCGA training set and 270 cases from the GEO testing set (Table I). Cases in the training set had an average age of 60.8±11.9 years, ranging from 19 to 90 years; 149 (29.1%) individuals were followed until death. Cases in the testing set had an average age of 60.1±10.3 years, ranging from 35 to 87 years; 88 (32.6%) individuals were followed until death.

To exclude a large number of genes unrelated to disease, we assessed the TCGA training set in two steps after quality control (Fig. 1A). First, of all 17711 genes with normalization, we used the WTT method to select the top 5% significant genes (n=886) from 43 pairs of tumor and matched adjacent normal tissue data. All 886 genes were significantly differentially expressed (all P≤1.15×10⁻⁵) and 199 genes were significant in univariate Cox regression analysis (P<0.05; Fig. 1B). Second, the SIS method was used for further dimension reduction. After iterative process for different genes and LASSO penalized regression with 10-fold cross-validation to select the best parameter, seven genes remained after selection. All of them were significantly overexpressed in tumor tissue (Fig. 1C). In addition, they were significantly associated with overall survival except LASP1 (P=0.056) (Fig. 1D). A Cox regression model was used to generate model coefficients. The biomarker signature model was calculated as risk score = 0.198×AATF + 0.244×APP + 0.252×GNPDA1 + 0.314×HPRT1 + 0.136×LASP1 + 0.110×P4HA1 - 0.388×ILF3. We categorized the patients into low-risk and high-risk groups and defined the cut-off value (score=0.36). This was selected by the optimum cut point according to the highest χ² value defined by Kaplan-Meier survival analysis and log-rank test in the training test (26).

As a linear combination model of seven mRNAs, the risk score was significantly associated with the TCGA patient survival (HRunadjust = 2.79; 95% CI, 1.98–3.92; P=4.05×10⁻⁹) (Fig. 2A). In total, 16.4% in the low-risk group vs. 41.8% in the high-risk group were followed until death (χ²=38.78, P=4.76×10⁻¹⁰) (Fig. 2B). With the bootstrap adjustment for clinical characteristics, the results remained significant for all cases (HRadjust = 2.86; 95% CI, 1.99–4.12; P<0.0001) or cases with available HPV status (HRadjust = 3.17; 95% CI, 1.90–5.30; P<0.0001) (Table II).

In the GEO testing set, risk scores were calculated for each patient. Using the same cut-off value (score=0.36), the score showed a 2.05 times higher risk of death for the high-risk group compared to the low-risk group in univariate Cox regression (HRunadjust = 2.05; 95% CI, 1.35–3.11; P=7.98×10⁻⁴) (Fig. 2C). In total, 25.1% in the low-risk group vs. 46.3% in the high-risk group were followed until death (χ² = 11.62, P=6.53×10⁻⁴) (Fig. 2D). Results retained statistical significance with further adjustment for covariates, including age, sex, smoking status, HPV status and clinical stage (HRadjust = 1.94; 95% CI, 1.27–2.96; P=0.002) (Table II).

Furthermore, prognostic prediction ability for 5-year survival was evaluated. The time-dependent AUCs of risk scores for HNSCC cases were 0.73 (95% CI, 0.68–0.78; P<0.001) in the training set (Fig. 3A) and 0.66 (95% CI, 0.59–0.73; P<0.001) in the testing set (Fig. 3B). Besides, we combined the scores with clinical characteristics to see whether they could improve the predictive value. In the training set, prognostic score plus clinical characteristics had a higher AUC (AUC, 0.75; 95% CI, 0.70–0.80) than the clinical characteristics (age, sex and stage) alone (AUC, 0.57; 95% CI, 0.51–0.64) (Fig. 3C). The testing set also displayed improvement in AUC from 0.63 (95% CI, 0.57–0.70) to 0.72 (95% CI, 0.65–0.78) (Fig. 3D). In brief, the risk score could better distinguish HNSCC prognosis beyond clinical information alone.

Next, we examined whether the risk score could help improve prognostication in the two datasets separately by subgroup sensitivity analysis. HPV-positive HNSCC has been widely recognized as associated with better prognosis than HPV-negative HNSCC (27). The risk score could significantly distinguish patient prognosis among 274 cases with available HPV information in the training set (HPV⁺, P=8.65×10⁻⁶; HPV⁻, P=1.24×10⁻⁷; Fig. 4A and B) and 269 cases in the testing set regardless of HPV status (HPV⁺, P=0.004; HPV⁻, P=0.014; Fig. 4C and D).

In different tumor sites, high vs. low-risk score significantly distinguished outcomes in patients with tumor at larynx (training set, P=0.001; testing set, P=0.019), oropharynx (training set, P=0.0004; testing set, P=0.003) and oral cavity in the training set (P=1.80×10⁻⁵) (Fig. 5). However, the result was not significant in oral cavity of the testing set (P=0.178), possibly due to the failure to distinguish patients with longer survival. In addition, we found that HPV-positive status improves survival in oropharyngeal HNSCC but not in non-oropharyngeal HNSCC (training set, P<0.001; testing set, P=0.010; Fig. 5E and F), which was consistent with a previous report (28).