A key genomic signature associated with lymphovascular invasion in head and neck squamous cell carcinoma

Abstract

Background: Lymphovascular invasion (LOI), a key pathological feature of head and neck squamous cell carcinoma (HNSCC), is predictive of poor survival; however, the associated clinical characteristics and underlying molecular mechanisms remain largely unknown.

Methods: We performed weighted gene co-expression network analysis to construct gene co-expression networks and investigate the relationship between key modules and the LOI clinical phenotype. Functional enrichment and KEGG pathway analyses were performed with differentially expressed genes. A protein-protein interaction network was constructed using Cytoscape, and module analysis was performed using MCODE. Prognostic value, expression analysis, and survival analysis were conducted using hub genes; GEPIA and the Human Protein Atlas database were used to determine the mRNA and protein expression levels of hub genes, respectively. Multivariable Cox regression analysis was used to establish a prognostic risk formula and the areas under the receiver operating characteristic curve (AUCs) were used to evaluate prediction efficiency. Finally, potential small molecular agents that could target LOI were identified with DrugBank.

Results: Ten co-expression modules in two key modules (turquoise and pink) associated with LOI were identified. Functional enrichment and KEGG pathway analysis revealed that turquoise and pink modules played significant roles in HNSCC progression. Seven hub genes (CNFN, KIF18B, KIF23, PRC1, CCNA2, DEPDC1, and TTK) in the two modules were identified and validated by survival and expression analyses, and the following prognostic risk formula was established: [risk score = EXP_DEPDC1 * 0.32636 + EXP_CNFN * (- 0.07544)]. The low-risk group showed better overall survival than the high-risk group (P < 0.0001), and the AUCs for 1-, 3-, and 5-year overall survival were 0.582, 0.634, and 0.636, respectively. Eight small molecular agents, namely XL844, AT7519, AT9283, alvocidib, nelarabine, benzamidine, L-glutamine, and zinc, were identified as novel candidates for controlling LOI in HNSCC (P < 0.05).

Conclusions: The two-mRNA signature (CNFN and DEPDC1) could serve as an independent biomarker to predict LOI risk and provide new insights into the mechanisms underlying LOI in HNSCC. In addition, the small molecular agents appear promising for LOI treatment.

Keywords: Head and neck squamous cell carcinoma; Hub genes; Lymphovascular invasion; TCGA; Weighted gene co-expression network analysis.

Fig. 1

Determination of soft-threshold power in WGCNA. a Scale-free index analysis for soft-threshold power (β) in HNSCC. b Mean connectivity analysis for various soft-threshold powers. c Histogram depicting connectivity distribution when β = 7. d Checking scale-free topology when β = 7

Fig. 2

Visualization of WGCNA results. a mRNA clustering dendrogram obtained by hierarchical clustering of Topological Overlap Matrix (TOM)-based dissimilarity, with the corresponding module colors indicated by colored rows. Each colored row represents a color-coded module containing a group of highly connected mRNAs. Each color represents a module in the constructed gene co-expression network. b The heatmap depicts TOM among all genes in WGCNA. Light color represents low overlap and progressively darker color represents higher overlap

Fig. 3

Correlation analysis of module–trait associations and clinical characteristics. a The column corresponds to the LOI phenotypic trait. Heatmap of each cell contains the P value of that module and the LOI phenotypic trait. Correlations between the turquoise module with the LOI phenotypic trait (cor = 0.25; P = 5E− 07) and the pink module with the LOI phenotypic trait (cor = − 0.23; P = 4E− 06) were significant. b Bar plot of the significance level of 10 co-expression modules associated with LOI status. c and d Correlation analysis between gene significance of LOI status and module membership in the turquoise (c) and pink (d) modules

Fig. 4

GO function and KEGG pathway analyses. a GO enrichment analysis of the turquoise module in the biological process category. b GO enrichment analysis of the turquoise module in the KEGG pathway. c GO enrichment analysis of the pink module in the biological process category

Fig. 5

Hub genes identified by the PPI network. a and b PPI network interaction of DEGs in the turquoise (a) and pink (b) modules

Fig. 6

Prognostic value and expression analysis of seven hub genes in HNSCC. a Ten-year cumulative survival of HNSCC patients with or without LOI. b–h Ten-year survival analysis of CNFN (b), KIF18B (c), KIF23 (d), PRC1 (e), CCNA2 (f), DEPDC1 (g), and TTK (h). i mRNA expression levels of the seven hub genes in HNSCC samples (n = 519, red box) and normal tissue samples (n = 44, blue box) based on GEPIA. j Immunohistochemistry images of the seven hub genes based on the Human Protein Atlas database. k Protein expression levels analyzed by immunohistochemistry based on the Human Pathology Atlas database. **P < 0.01 and *P < 0.05

Fig. 7

Distribution of risk score, survival status, and time-dependent ROC analysis of the integrated two-mRNA signature. a Risk score distribution b Overall survival (OS) status of 330 patients. c Kaplan–Meier curve of OS between the low- and high-risk groups split by the median risk score. d Time-dependent ROC analysis for 1-, 3-, and 5-year OS probability