Cancer associated fibroblast derived gene signature determines cancer subtypes and prognostic model construction in head and neck squamous cell carcinomas

Abstract

Background: Head and neck squamous cell carcinomas (HNSCC) are the most common type of head and neck cancer with an unimproved prognosis over the past decades. Although the role of cancer-associated-fibroblast (CAF) has been demonstrated in HNSCC, the correlation between CAF-derived gene expression and patient prognosis remains unknown.

Methods: A total of 528 patients from TCGA database and 270 patients from GSE65858 database were contained in this study. After extracting 66 CAF-related gene expression data from TCGA database, consensus clustering was performed to identify different HNSCC subtypes. Limma package was used to distinguish the differentially expression genes (DEGs) between these subtypes, followed by Lasso regression analysis to construct a prognostic model. The model was validated by performing Kaplan-Meier survival, ROC and risk curve, univariate and multivariate COX regression analysis. GO, KEGG, GSEA, ESTIMATE and ssGSEA analyses was performed to explort the potential mechanism leading to different prognosis.

Results: Based on the 66 CAF-related gene expression pattern we stratitied HNSCC patients into two previously unreported subtypes with different clinical outcomes. A prognostic model composed of 15 DEGs was constructed and validated. In addition, bioinformatics analysis showed that the prognostic risk of HNSCC patients was also negatively correlated to immune infiltration, implying the role of tumor immune escape in HNSCC prognosis and treatment option.

Conclusions: The study develops a reliable prognostic prediction tool and provides a theoretical treatment guidance for HNSCC patients.

Keywords: Cancer subtypes; Cancer-Associated-Fibroblast; Gene signature; Head and Neck Squamous Cell Carcinoma; Prognostic Model.

FIGURE 1

Analysis process of this study.

FIGURE 2

Identification of two HNSCC subtypes based on CAF‐related gene signature expression. (A) Consensus matrix of HNSCC patients (n = 499). (B) Cumulative distribution function. (C) The relative change of the area under CDF curve. (D) Kaplan–Meier curve indicated that patients of cluster 1 have a worse prognosis, p = 0.048 by log‐rank test (Cluster1, n = 231, Cluster2, n = 268). (E) StromalScore and TMEScore were significantly higher in cluster 1 than cluster 2, while ImmuneScore was equivalent between two clusters. **p < 0.001. ns, no significant difference. Student's t‐test.

FIGURE 3

Construction of a 15‐gene signature prognostic model for HNSCC. (A) Heatmap of DEGs between two clusters. (B) Forest plot showed DEGs related with patients' survival status. (C) Lasso regression model. (D) Cross validation showed that 15 was the optimum λ value. (E) Table of 15 genes and their corresponding coefficients included in the prognostic model. (F) Heatmap of 15 genes' expression in two risk groups and the relationship between risk score and clinical features was shown repsectively. *p < 0.05. Chi‐square test.

FIGURE 4

Combined Kaplan–Meier curve, ROC curve, risk curve and independent prognostic analysis demonstrates a good performance of the HSNCC prognostic model. (A, B) Kaplan–Meier curve showed that patients in the high risk group was associated with poor prognosis in both training (n = 499 patients) (A) and validation cohorts (n = 270 patients) (B). (C, E) ROC curves of the training cohort (C) and validation cohort (E). (D, F) Risk curve showed that the number of deaths (bottom, red dots) increased as the risk score increased (top) in both training (D) and validation cohorts (F). (G, H) Univariate COX regression (top) and multivariate COX regression (bottom) indicated risk score as an independent prognostic factor in both training (G) and validation cohorts (H). AUC, Area Under Curve.

FIGURE 5

Construction of a Nomogram. (A) A nomogram containing risk score to predict 1‐, 3‐, 5‐year survival of HNSCC patients. (B–D) 1‐year (B), 3‐year (C), 5‐year (D) surivval calibration plots showed that the prediction calibration curve (red) was close to the standard curve (gray), verified a good predictive ability of the nomogram. OS, Overall Survival.

FIGURE 6

Correlation Between Prognostic Risk and Tumor Immune Microenvironment. (A, B) The correlation between risk score and StromalScore/ImmuneScore/TMEScore in training cohort (A) or validation cohort (B) is given as indicated. (C, D) The number of immune cells (top) and the activity of immune‐related processes (bottom) in training cohort (C) and validation cohort are shown (D).