A seven-gene prognostic signature for rapid determination of head and neck squamous cell carcinoma survival.

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer and displays divergent clinical outcomes. Prognostic biomarkers might improve risk stratification and survival prediction. We aimed to investigate the prognostic genes associated with overall survival. A two-step gene selection method was used to develop a seven-gene-based prognostic model based on the training set collected from The Cancer Genome Atlas (TCGA). In addition, the prognostic model was validated in an independent testing set from Gene Expression Omnibus (GEO). The score based on the model successfully distinguished HNSCC survival into high-risk and low-risk groups in the training set (HR, 2.79; 95% CI, 1.98-3.92; P=4.05x10-9) and the testing set (HR, 2.05; 95% CI, 1.35-3.11; P=7.98x10-4). In addition, the score could significantly predict 5-year survival by ROC curves (AUCs for training set, 0.73; testing set, 0.66). Combining risk scores with clinical characteristics improved the AUCs beyond using clinical characteristics alone (training set, from 0.57 to 0.75; testing set, from 0.63 to 0.72). A subgroup sensitivity analysis with HPV status and tumor sites revealed that the risk score was significant in all subgroups except oral cavity tumors of the testing set. Furthermore, HPV-positive status improves survival in oropharyngeal HNSCC but not non-oropharyngeal HNSCC. In conclusion, the seven-gene prognostic signature is a reliable and practical prognostic tool for HNSCC. This approach can add prognostic value to clinical characteristics and provides a new possibility for individualized treatment.

Introduction

Head and neck squamous cell carcinoma (HNSCC) constitutes approximately 4% of all new cancer diagnoses in the United States, with approximately 62,000 new cases in 2016 (1). Each year approximately 600,000 patients are affected worldwide (2). Importantly, the 5-year survival rate of HNSCC patients is only 40–50% (3). The high mortality rate is attributable to a high rate of late diagnosis, and the survival rate for cases in late stages is only 34.9% (4). These outcomes demonstrate the need for prognostic biomarkers to help predict patient outcome and outline individualized treatment plans. Age, clinical stage and smoking status are characteristics emerging as important contributors to clinical outcome that may also help us improve survival prediction (5–7). However, the traditional clinical information has limited prediction ability due to the complex molecular regulation mechanism in cancer.

Recently, the clinical importance of messenger RNA (mRNA) expression has been reported in various types of cancer including HNSCC. They play important roles in a variety of physiological and pathological processes, such as development, differentiation, cell proliferation, apoptosis and stress responses (8). Therefore, characterization of key genes in different tumors is essential not only for an urgent requirement of precision medicine (9) but also for preclinical and pharmaceutical research (10).

Prognostic models in HNSCC have been described using different biomarkers such as somatic mutations (11), microRNAs (12,13) and proteins (4) but limited studies focus on mRNA expression according to Cancer Genetics Web (14). De Cecco et al (15) reported a gene expression survival predictor using HNSCC microarray data based on a semi-supervised survival method involving principal component method (16). However, the model comprised 172 genes and was complicated for further interpretation. Now that transcriptome sequencing technologies (RNA-Seq) are being applied widely, there is a more ideal platform for cancer genetic studies (17). In addition, The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) repositories provide abundant HNSCC case resources, which may be useful to explore reliable biomarkers.

In this study, we investigated the prognostic value of seven gene expression biomarkers (AATF, APP, GNPDA1, HPRT1, LASP1, P4HA1 and ILF3) for HNSCC overall survival. The cases mainly included oropharyngeal, laryngeal and oral squamous cell carcinoma. TCGA cohort was used as the training set to generate the prognostic model by a two-stage variable selection. We also used an independent external testing set to validate the robustness and reproducibility of the model.

Materials and methods

Study population

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.

Weighted t-test (WTT) method as the first step for gene selection

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, the Cox proportional hazards model is given by λ(t|Z) = λ0(t)exp(β) and the survival function is S(t|Z) = exp{-λ0(t)exp(β)}, where λ0(t) is the basic hazard function, β is the regression coefficient and λ0(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. 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 = w 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 for the kth gene. Afterwards, 1,000 total permutations were performed and d could be got for the ith permutation. Then, we calculated the averaged order statistics, , across all 1,000 permutations. A gene was labeled as significant when was at the top 5%.

Sure independence screening (SIS) as the second step for gene selection

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.

Statistical analysis

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.

Results

Demographic and clinical characteristics

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.

Table I.

Demographic and clinical characteristics of HNSCC patients.

Characteristics	Training set (n=512)	Testing set (n=270)
Median follow-up time (years)	4.35	4.95
Censor rate (%)	70.8	67.4
Age, mean ± SD (years)	60.8±11.9	60.1±10.3
Sex, n (%)
Male	376 (73.4)	223 (82.6)
Female	136 (26.6)	47 (17.4)
Smoking status, n (%)
Never	115 (22.5)	48 (17.8)
Current/former	385 (75.2)	222 (82.2)
NA[a]	12 (2.3)	0 (0)
Tumor site, n (%)
Oropharynx[b]	80 (15.6)	102 (37.8)
Larynx	114 (22.3)	48 (17.8)
Oral cavity[b]	308 (60.2)	83 (30.7)
Others	10 (2)	37 (13.7)
HPV status, n (%)
Positive	35 (6.8)	73 (27.0)
Negative	241 (47.1)	196 (72.6)
NA[a]	236 (46.1)	1 (0.4)
T classification, n (%)
T1	48 (9.4)	35 (13.0)
T2	130 (25.4)	80 (29.6)
T3	99 (19.3)	58 (21.5)
T4	172 (33.6)	97 (35.9)
TX or NA[a]	63 (12.3)	0 (0)
N classification, n (%)
N0	174 (34)	94 (34.8)
N1	66 (12.9)	32 (11.9)
N2	165 (32.2)	132 (48.9)
N3	8 (1.6)	12 (4.4)
NX or NA[a]	99 (19.3)	0 (0)
M classification, n (%)
M0	484 (94.5)	263 (97.4)
M1	4 (0.8)	7 (2.6)
MX or NA[a]	24 (4.7)	0 (0)
TNM stage, n (%)
I	20 (3.9)	18 (6.7)
II	97 (18.9)	37 (13.7)
III	104 (20.3)	37 (13.7)
IV	278 (54.3)	178 (65.9)
NA[a]	13 (2.5)	0 (0)
Grade, n (%)
1	61 (11.9)	–
2	300 (58.6)	–
3	122 (23.8)	–
4	7 (1.4)	–
NA[a]	22 (4.3)	270 (100)
Neoadjuvant treatment, n (%)
Yes	10 (1.9)	–
No	502 (98.1)	–
NA[a]	0 (0)	270 (100)

NA, not available

oropharynx also includes tonsil and base of tongue; oral cavity also includes oral tongue, buccal mucosa, lip, alveolar ridge, hard palate and floor of mouth.

Development of biomarker signature model

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−5) 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 χ2 value defined by Kaplan-Meier survival analysis and log-rank test in the training test (26).

Figure 1.

(A) Work flow of gene selection steps, with the number of genes remaining. (B) Heatmap of the seven genes using 43 pairs of tumor and matched adjacent normal HNSCC tissue data in the training set. The upper half is the expression of normal tissues, and the lower half is the expression of tumor tissues. (C) Volcano plot showing the prognostic significance of the 886 genes selected by the WTT method. This plot depicts the hazard ratio on the x-axis and statistical significance on the y-axis, as evaluated by the Cox regression model. The blue dashed line indicates a P-value of 0.05. (D) Hazard ratio with 95% CI of the seven genes in univariable Cox regression analysis of the training set.

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−9) (Fig. 2A). In total, 16.4% in the low-risk group vs. 41.8% in the high-risk group were followed until death (χ2=38.78, P=4.76×10−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).

Figure 2.

(A and B) Training set. (C and D) Testing set. Upper left panel: risk score distribution of the seven-gene model classifier and patient survival status. All scores are standardized (‘score’ = score - 0.36) to make the low-risk group negative and high-risk group positive. Lower half panel: heatmap showing expression of the seven genes among tumor patients. Right panel: Kaplan-Meier analysis for the patients. The patients are divided into low-risk (red) and high-risk (blue) groups.

Table II.

Multivariable Cox regression analysis of clinical characteristics and risk score.

	Training set
	All cases (n=512)		Cases with HPV status (n=276)		Testing set (n=270)
Characteristics	HR (95% CI)	P-value	HR (95% CI)	P-value	HR (95% CI)	P-value
High risk score	2.86 (1.99–4.12)	<0.0001	3.17 (1.90–5.30)	<0.0001	1.94 (1.27–2.96)	0.002
Age (per year)	1.02 (1.00–1.04)	0.039	1.01 (0.99–1.04)	0.141	1.03 (1.00–1.06)	0.020
Gender (Female)	1.11 (0.76–1.62)	0.578	1.07 (0.64–1.79)	0.797	1.02 (0.57–1.79)	0.952
Smoking status
(Current/former smoker)	1.26 (0.80–1.98)	0.324	1.56 (0.87–2.79)	0.133	0.92 (0.49–1.75)	0.815
Clinical stage (per stage)	1.09 (0.91–1.31)	0.331	1.02(0.82–1.25)	0.878	1.77 (1.21–2.59)	0.003
HPV status (positive)	–	–	0.65 (0.29–1.47)	0.301	0.43 (0.24–0.79)	0.006

Estimation of confidence intervals and P-values are based on a bootstrap estimate (500 resamples) of the variance-covariance matrix.

Validation of the prognostic signature

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−4) (Fig. 2C). In total, 25.1% in the low-risk group vs. 46.3% in the high-risk group were followed until death (χ2 = 11.62, P=6.53×10−4) (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.

Figure 3.

(A and B) Training set. (C and D) Testing set. Upper half panel: time-dependent ROC curves are used to evaluate patient survival with the risk score on different time, obtained by the nearest neighbor method. Lower half panel: ROC curves for clinical model and improved model. Clinical model contains age, sex, smoking status and clinical stage. Improved model contains the risk score and characteristics of clinical model.

Subgroup sensitivity analysis with HPV status and tumor site

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−6; HPV−, P=1.24×10−7; 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).

Figure 4.

Subgroup sensitivity analysis with HPV status. (A) Training set (HPV+). (B) Training set (HPV−). (C) Testing set (HPV+). (D) Testing set (HPV−). Cross tables of tumor sites and risk group proportions in each figure were summarized and tested by the Fishers exact test.

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−5) (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).

Figure 5.

Subgroup sensitivity analysis with different tumor sites. (A and B) Laryngeal tumors. (C and D) Oral cavity tumors. (E and F) Oropharyngeal tumors. (A, C and E) Training set. (B, D and F) Testing set. Cross tables of HPV status and risk group proportions in each figure were summarized and tested by Fishers exact test.

Discussion

In the present study, we developed an HNSCC prognostic risk model that includes seven mRNAs and validated it using an independent external data set. Integrating multiple biomarkers into an aggregated model could improve prognostic value compared with single biomarker (29). Results showed that the risk score was significantly associated with patient overall survival. HNSCC patients with higher risk scores tended to have a poorer clinical outcome. In addition, this score could improve model performance combined with clinical characteristics based on 5-year overall survival.

To screen out the survival-related genes from over 20,000 total genes, we used a two-stage screening method. The WTT method was used as the first step to identify a subset of genes that were not only differentially expressed in the tumor and matched normal tissues but also had an impact on patient survival by weighting the clinical covariates. Using permutation procedures that have widely been used in biomedical data analysis (30), error rates can be controlled (22). The results show that all the selected genes were differentially expressed and a number (22.5%) were significantly associated with prognosis. The WTT method performs better than a traditional t-test, which will flag a large number of genes unrelated to disease. Afterwards, SIS was used to reduce the number of genes included in the final model as the second step. Compared with a traditional penalized regression like lasso or elastic net models, SIS is ungraded on the basis of penalized regression to reduce dimensionality from high to a moderate scale that is below the sample size (24). It has improved both speed and accuracy, and has a stronger association with disease (31). Of the seven genes in the training set, all of them are significantly differentially expressed in tumor and normal tissue. In addition, six of them are significant in a univariate Cox model and the last one shows suggestive significance (P=0.056). Our results show that the two-step gene selection method is amenable to deal with a high dimension problem.

Among the seven candidate genes, six have positive coefficients in the prognostic model and are associated with worse survival. AATF, also called Che-1, is a critical regulator of apoptosis driven by genes coding for PAR4 and p53 (32), and promotes tumor cell survival by sustaining mutant p53 transcription and inhibiting DNA damage response activation (33). It has activity in transcriptional regulation, cell cycle control, DNA damage responses, and in the execution of cell death programs (34). In addition, it interacts with NRAGE that has been found as a tumor marker in different cancers. AATF has been reported associated with multiple cancers, such as colon carcinoma (35), gastric cancer (36), hepatocellular carcinoma (37) and breast cancer (38).

APP was initially found to be associated with Alzheimers disease, but it also contributes to regulating cell growth, apoptosis, and motility of cancer cells (39). Several studies have confirmed APP as an invaluable marker for oral carcinogenesis that promotes the proliferation and carcinogenesis of oral squamous cell carcinoma (OSCC) (40–42). Notably, in addition to its effects on promoting oral carcinogenesis, APP expression could be negatively regulated by tea in OSCC, which has been demonstrated to be effective in preventing animal carcinogenesis in different experimental systems (41).

LASP1, a recognized cancer biomarker functioning in cell structure, physiological processes, and cell signaling, contributes to cancer aggressiveness by overexpression (43). Increased LASP1 levels occur in OSCC and more than ten other tumor types (44). It appears to involved in regulation of cancer cell metastatic propensity and perturb the architecture and dynamics of focal adhesion that triggers cell migration and invasion (45). In OSCC, LASP1 plays an essential role in tumor cellular growth by mediating G2/M transition.

This is not the first report for P4HA1 that was associated with HNSCC prognosis (46). It is involved in hydroxylation of collagen fibers and upregulated by HIF1 under hypoxic conditions directly that drive a series of different biological processes related with malignant progression. P4HA1 modulates target genes in cancer cell growth and tumor progression (47) and its expression increases during the invasion and metastasis of breast cancer and hepatoma as well (48,49).

Furthermore, GNPDA1 is an allosteric enzyme that catalyses the reversible conversion of D-glucosamine-6-phosphate into D-fructose-6-phosphate and ammonium (50). It has been reported upregulated in colorectal cancer cells with western blotting and immunofluorescence assay (51). The protein encoded by HPRT1 is a transferase, which plays a central role in the generation of purine nucleotides through the purine salvage pathway. It still needs further experiments to validate its prognostic value.

In contrast, ILF3 confers an onco-protective effect. Downregulation of ILF3 can delay cell cycle progression, inhibit cell proliferation and reduce tumorigenic capacity in vivo (52). ILF3 is also involved in HPV-induced oncogenesis and p53-mediated apoptosis. It is a positive regulator of HPV E6 expression and its depletion leads to the accumulation of active p53 (53). Since HPV is effective in HNSCC, targeting on this gene may be useful to control the cancer.

The present study includes some limitations. First, in the subgroup analysis, the results may not strongly be robust due to small sample size of some groups. Second, due to the different experimental methods (RNASeq versus microarray) used between the two data sets, some bias may exist. Third, the differences between the two populations may indicate a need for further validation in another independent cohort for the current prognostic signature. Finally, the prognostic value of the seven genes in HNSCC still warrants further biological functional experiments.

In conclusion, our results showed that the seven-gene prognostic score significantly distinguishes HNSCC patients prognosis and predicts 5-year overall survival in both training and testing sets. Thus, this score may be a novel biomarker based on gene expression levels and it warrants further investigation for establishing its relevance for clinical application.