Nine-gene signature and nomogram for predicting survival in patients with head and neck squamous cell carcinoma

Abstract

Background: Head and neck squamous cell carcinomas (HNSCCs) are derived from the mucosal linings of the upper aerodigestive tract, salivary glands, thyroid, oropharynx, larynx, and hypopharynx. The present study aimed to identify the novel genes and pathways underlying HNSCC. Despite the advances in HNSCC research, diagnosis, and treatment, its incidence continues to rise, and the mortality of advanced HNSCC is expected to increase by 50%. Therefore, there is an urgent need for effective biomarkers to predict HNSCC patients' prognosis and provide guidance to the personalized treatment. Methods: Both HNSCC clinical and gene expression data were abstracted from The Cancer Genome Atlas (TCGA) database. Intersecting analysis was adopted between the gene expression matrix of HNSCC patients from TCGA database to extract TME-related genes. Differential gene expression analysis between HNSCC tissue samples and normal tissue samples was performed by R software. Then, HNSCC patients were categorized into clusters 1 and 2 via NMF. Next, TME-related prognosis genes (p < 0.05) were analyzed by univariate Cox regression analysis, LASSO Cox regression analysis, and multivariate Cox regression analysis. Finally, nine genes were selected to construct a prognostic risk model and a prognostic gene signature. We also established a nomogram using relevant clinical parameters and a risk score. The Kaplan-Meier curve, survival analysis, time-dependent receiver operating characteristic (ROC) analysis, decision curve analysis (DCA), and the concordance index (C-index) were carried out to assess the accuracy of the prognostic risk model and nomogram. Potential molecular mechanisms were revealed by gene set enrichment analysis (GSEA). Additionally, gene correlation analysis and immune cell correlation analysis were conducted for further enriching our results. Results: A novel HNSCC prognostic model was established based on the nine genes (GTSE1, LRRN4CL, CRYAB, SHOX2, ASNS, KRT23, ANGPT2, HOXA9, and CARD11). The value of area under the ROC curves (AUCs) (0.769, 0.841, and 0.816) in TCGA whole set showed that the model effectively predicted the 1-, 3-, and 5-year overall survival (OS). Results of the Cox regression assessment confirmed the nine-gene signature as a reliable independent prognostic factor in HNSCC patients. The prognostic nomogram developed using multivariate Cox regression analysis showed a superior C-index over other clinical signatures. Also, the calibration curve had a high level of concordance between estimated OS and the observed OS. This showed that its clinical net can precisely estimate the one-, three-, and five-year OS in HNSCC patients. The gene set enrichment analysis (GSEA) to some extent revealed the immune- and tumor-linked cascades. Conclusion: In conclusion, the TME-related nine-gene signature and nomogram can effectively improve the estimation of prognosis in patients with HNSCC.

Keywords: GEO; TCGA; bioinformatics analysis; gene signature; head and neck squamous cell carcinoma; prognosis.

FIGURE 1

Flow chart of the present study.

FIGURE 2

(A) Heatmap to show TME-related differentially expressed genes between HNSCC and normal tissue samples. (B) Volcano plot displays 1,335 TME-related differentially expressed genes in TCGA HNSCC cohort. (C) Non-negative matrix factorization (NMF) clustering was conducted to classify HNSCC patients into 2 to 10 different subtypes, and relevant heatmaps were generated. (D) NMF rank survey with multiple parameters (including dispersion, cophenetic, residuals, evar, rss, and silhouette coefficients) to determine the optimal value for consensus clustering.

FIGURE 3

(A) Consensus map of the NMF clustering data of TCGA HNSCC cohort. Rank = 2 means that HNSCC patients were separated into two groups. (B) Survival curve of OS (P=<0.001) in cluster 1 and 2. (C) Sankey plot to show the association between different subtypes and immune subtypes. (D) Differential analysis of tumor-infiltrating immune cells was conducted using MCP-counter.

FIGURE 4

(A) Deviance plot of partial likelihood by LASSO Cox regression analysis. Red dots: the partial likelihood of deviance values; gray lines: the standard error (SE); vertical dotted line on the left: the optimal values by minimum criteria; vertical dotted line on the right: 1−SE criteria. (B) LASSO coefficient profiles of the 121 TME-related genes in HNSCC.

FIGURE 5

(A) Time‐dependent ROC analysis and survival analysis of the nine‐gene signature in the training set. (B) Time‐dependent ROC analysis and survival analysis of the nine‐gene signature in the testing set. (C) Time‐dependent ROC analysis and survival analysis of the nine‐gene signature in the whole set.

FIGURE 6

(A) Nomogram for predicting the OS of HNSCC patients. For each patient, lines are drawn downward to assess the points received from the seven prognostic factors in the nomogram. The sum of these points is shown on the “Total points” axis. A line is drawn downward to determine the 1‐, 3‐, and 5‐year OS of HNSCC patients. (B) Calibration plot for the internal validation of the nomogram. The Y‐axis exhibits actual survival. The X‐axis exhibits nomogram‐estimated survival. The dotted line (45° diagonal line) signifies full agreement between actual and observed probabilities. (C) Comparison of the AUC values of the nomogram and risk score, age, gender, grade, and stage. (D) Clinical net benefit of the nomogram, risk score, and other clinical features (including age, gender, grade, and stage).

FIGURE 7

(A) Enrichment plots showing epidermis development (blue), epidermal cell differentiation (green), cornification (red), and keratinization (purple) gene enrichment in the high-risk group. (B) Enrichment plots showing the activation of immune response (red), adaptive immune response (orange), adaptive immune response based on somatic recombination of immune receptors built (green), B-cell activation (purple), and antigen receptor-mediated signaling pathway (blue) enrichment in the low-risk group.

FIGURE 8

(A,B) Subgroup analysis stratified by different tumor stages (stages I–II and stages III–IV). (C,D) Subgroup analysis stratified by different tumor grades (grades 1–2 and grades 3–4). (E,F) Subgroup analysis stratified by different age (over 65 and under 65 years). (G,H) Subgroup analysis stratified by different genders (male and female).

FIGURE 9

(A) Time‐dependent ROC analysis and survival analysis of the nine-gene signature (TME signature). (B–D) Time‐dependent ROC analysis and survival analysis of three-gene signatures established by other researchers.

FIGURE 10

(A–B) C-index results and percentile of scores of the four different signatures. (C) AUC values for predicting 1-, 3-, and 5-year OS for patients in the GSE16076 dataset. (D) Survival analysis for patients in the GSE16076 dataset.

FIGURE 11

(A) Heatmap showing the correlation between risk score and immune cells. (B) Spearman correlation analysis of TME-related genes in PTC. The number on the right vertical axis represents the Spearman correlation coefficient between two genes. The black asterisk (*) inside the circle indicates p < 0.05. Red indicates positive correlation. Blue indicates negative correlation. The stronger the correlation, the larger the circle and the deeper the color.