An Epithelial-Mesenchymal Transition Hallmark Gene-Based Risk Score System in Head and Neck Squamous-Cell Carcinoma.

BACKGROUND: Epithelial-to-mesenchymal transition (EMT) program plays a critical role in cancer. Thus, we attempted to generate a risk score system according to the expression pattern of different EMT hallmark genes in head and neck squamous-cell carcinoma (HNSC).
METHODS: Differentially expressed EMT hallmark genes were screened to generate a risk score (RS) on TCGA HNSC dataset. The relative prognostic value of the RS compared to clinicopathological characteristics was explored using multivariable Cox analysis. Functional enrichment analysis was performed to reveal the biological characteristics. An external dataset was applied to validate the prognostic value of the RS.
RESULTS: Nine genes constituted the EMT hallmark gene-based RS, which is significantly associated with poor prognosis and could successfully divide patients with HNSC into high- and low-risk groups. The RS was also an independent prognostic indicator compared to routine clinical factors.
CONCLUSION: We proposed and validated a nine-EMT hallmark gene-based risk score system in HNSC.

Introduction

Head and neck cancer is the seventh most common malignancy worldwide, with the greater being head and neck squamous-cell carcinoma (HNSC).1 Despite advances in treatments, the prognosis of HNSC remains poor with mortality rates of 40–50%.2 HNSC is a heterogeneous type of disease at phenotypic and genetic levels.3 The current clinical decision-making system in HNSC is mainly based on phenotypic heterogeneity, such as the American Joint Committee on Cancer TNM staging system4 and tumor grade. It is imperative to identify individuals at high-risk with the same phenotypes by revealing the genetic heterogeneity.

Epithelial-to-mesenchymal transition (EMT) is considered one of the hallmarks of cancer.5 The emerging evidence has shown that EMT program contributes to the induction of cancer stem cells, immune escape during cancer progression, and drug resistance in various types of cancers,6–8 including HNSC.9 Thus, the expression pattern of different EMT hallmark genes may be one of the critical genetic heterogeneity of cancers. We assumed that there is an EMT hallmark gene-based risk score system as a prognostic indicator in HNSC. To investigate the hypothesis, we used the datasets of HNSC from TCGA, including clinical information and gene expression profiles, to generate an EMT hallmark gene-based risk score system for predicting prognosis for patients with HNSC, and validated it on another independent dataset.

Materials and Methods

Data Processing

The RNA sequencing (RNA-seq) data (displayed as raw counts) and clinical information of HNSC in The Cancer Genome Atlas () were downloaded to generate an EMT hallmark gene-based risk score (RS) to predict prognosis. Another HNSC data set GSE6585810 based on platform of GPL10558 was downloaded from Gene Expression Omnibus () and used to validate the RS. The RNA-seq data was normalized using voom function from limma package11 in R software (version 4.0.2) (). The gene expression profiles in GSE65858 were normalized by the contributor. If multiple probes correspond to a gene, then the average value of these probes is considered as the expression value of this gene. EMT hallmark gene set included 200 genes () was obtained from the Molecular Signatures Database (version 7.2).12,13 The present study only included 192 EMT hallmark genes because their expression values are available in both TCGA and GSE65858. The workflow of the present study is displayed in Figure 1.

Figure 1

The workflow of the present study.

Screening Differentially Expressed Genes (DEGs)

The expression profiles of the 192 EMT hallmark genes were extracted from the TCGA HNSC dataset, and subsequently used to screen the DEGs in HNSC compared to healthy tissues using limma package. Genes with P adjusted by false discovery rate < 0.05 and log2 (fold change) >1 were considered significant.

Cox Regression and Least Absolute Shrinkage and Selection Operator Analysis

The expression profiles of DEGs and overall survival data were used to perform univariable Cox regression to identify the prognosis-associated EMT hallmark genes. Subsequently, the expression profiles of prognosis-associated EMT hallmark genes were performed with least absolute shrinkage and selection operator (LASSO) analysis using glmnet () R package to select the optimal prognostic EMT hallmark genes. Thus, the EMT hallmark gene-based risk score (RS) was created as:

RS = Exprgene1*Coefgene1 + Exprgene2*Coefgene2+ Exprgene2*Coefgene2+ …

The “Coef” is the regression coefficient of gene and is derived from the LASSO Cox regression, and “Expr” indicates the expression values of the gene. Each patient with HNSC got a RS, and was divided into the high- or low-risk group according to the median RS. The OS between the two different risk groups were compared using log-rank method. In addition, the routine clinical factors were included in the multivariable Cox regression analysis to assess whether the RS is an independent prognostic factor.

Gene Set Enrichment Analysis (GSEA)

To explore the biological state of high-risk group patients, GSEA12,14 was performed using the GSEA JAVA program (version 4.0.1) (). The hallmark gene set included 50 gene sets obtained from the Molecular Signatures Database (version 7.2)12,13 used as the reference gene set. Gene sets with nom p < 0.05 after performing 1000 permutations were considered to be significantly enriched.

Validation of the EMT Hallmark Gene-Based RS

As it was in the TCGA HNSC data set, each patient in GSE65858 got an RS according to the above formula, and was divided into high- or low-risk groups. The OS between the two different risk groups were compared and the multivariable Cox regression analysis was carried out.

Statistical Analysis

All analyses were performed using R software (version 4.0.2). The unpaired t-test from limma package was used to screen DEGs. Kaplan–Meier survival analysis and Log rank test were used to compare survival between the two groups of patients. Time-receiver operating characteristic (tROC) curve analysis was performed using the timeROC package (). All tests were two-sided and p < 0.05, unless otherwise stated, was considered to indicate statistical significance.

Results

Multiple EMT Hallmark Genes Upregulated in HNSC

A total of 93 EMT hallmark genes () were differentially expressed in HNSC compared to healthy paracancer tissue, including 16 downregulated and 77 upregulated genes (Figure 2A). This indicates that EMT plays a crucial role in promoting HNSC due to the fact that most EMT hallmark genes are upregulated. The expression heat map of the DEGs shows that it has a promising effect of distinguishing tumor from paracancer tissue (Figure 2B).

Figure 2

Differentially expressed EMT hallmark genes. (A) The volcano plot of the differentially expressed genes. Red indicates up-regulated, and blue indicates down-regulated. (B) The expression heat map of the differentially expressed EMT hallmark genes.

EMT Hallmark Gene-Based RS as an Independent Prognostic Factor

After univariable Cox regression analysis, nineteen EMT hallmark genes were identified as the prognosis-associated genes (Table 1). Unsurprisingly, most of (17 from 19) the prognosis-associated EMT hallmark genes show an association with poor prognosis in HNSC. Subsequently, nine EMT hallmark genes (SFRP1, TGFBR3, DKK1, PCOLCE2, PTX3, CAP2, PLOD2, VEGFC, and IL6) were considered as optimal features through LASSO Cox analysis (Figure 3A). Thus, all patients got RS according to the coefficients. The RS is significantly associated with poor prognosis (Hazard Rate (HR) = 3.254, 95% CI = 2.367–4.473, p < 0.001). The RS showed promising prognostic value with AUC approximately 0.7 (Figure 3B), and the AUC of 5-year tROC was 0.660 (Figure 3C). The high-risk group HNSC patients showed significantly shorter OS than the low-risk group HNSC patients (Figure 3D). Furthermore, the RS remained independent compared to some routine clinical factors, including TNM staging system, tumor grade, and tumor primary subdivision (Figure 4).

Table 1

The Results of Univariable and Least Absolute Shrinkage and Selection Operator Cox Analysis

Predictor	Univariable Cox Analysis				LASSO Analysis
	β	HR	HR95% CI	P value	Coefficient
BASP1	0.121	1.128	1.021–1.246	0.017
CAP2	0.124	1.132	1.025–1.249	0.014	0.06442
DKK1	0.132	1.141	1.080–1.204	0.000	0.08468
FAP	0.096	1.101	1.012–1.198	0.025
FSTL3	0.092	1.096	1.006–1.195	0.036
IL6	0.082	1.086	1.016–1.160	0.015	0.01665
INHBA	0.106	1.112	1.032–1.198	0.005
ITGA5	0.188	1.207	1.078–1.351	0.001
NT5E	0.117	1.125	1.039–1.217	0.004
PCOLCE2	0.106	1.112	1.049–1.178	0.000	0.07212
PLOD2	0.204	1.226	1.082–1.389	0.001	0.0227
PTX3	0.121	1.129	1.062–1.199	0.000	0.06458
SERPINE1	0.114	1.12	1.033–1.215	0.006
SERPINH1	0.146	1.158	1.006–1.332	0.041
SFRP1	−0.061	0.941	0.895–0.989	0.017	−0.06221
TGFBI	0.092	1.096	1.009–1.191	0.030
TGFBR3	−0.11	0.896	0.810–0.991	0.033	−0.06494
TNFRSF12A	0.194	1.214	1.069–1.379	0.003
VEGFC	0.092	1.097	1.019–1.181	0.014	0.00709

Abbreviations: LASSO, least absolute shrinkage and selection operator; HR, hazard ratio; CI, confidence interval.

Figure 3

The nine-EMT hallmark gene-based risk score in the head and neck squamous-cell carcinoma data set from The Cancer Genome Atlas. (A) Nine genes were considered as the optimal features in the least absolute shrinkage and selection operator Cox analysis. (B) The time-dependent receiver operating characteristic (ROC) curve analysis for the risk score. (C) 5-year The time-dependent ROC curve analysis for the risk score. (D) The Kaplan–Meier curves with the Log rank test of the high- and low-risk groups.

Figure 4

The results of the multivariable Cox analysis for the nine-EMT hallmark gene-based risk score (RS) and routine clinical factors. *P < 0.05, **P < 0.01, and ***P < 0.001.

Biological Phenotypes of High-Risk HNSC

The results of GSEA showed that the EMT hallmark gene set was significantly enriched in high-risk HNSC (Figure 5A). In addition to this, hallmark gene set of angiogenesis (Figure 5B), coagulation (Figure 5C), glycolysis (Figure 5D), hypoxia (Figure 5E), MTORC1 signaling (Figure 5F), unfold protein response (Figure 5G), and UV response up (Figure 5H) were also enriched in high-risk HNSC.

Figure 5

The results of gene set enrichment analysis for the high-risk head and neck squamous-cell carcinoma. Eight hallmark gene sets enriched in the samples of high-risk head and neck squamous-cell carcinoma, including (A) epithelial–mesenchymal transition, (B) angiogenesis, (C) coagulation, (D) glycolysis, (E) hypoxia, (F) MTORC1 signaling, (G) unfold protein response, and (H) UV response up.

The RS Was Validated in an Independent Data Set

As it was in the TCGA HNSC data set, the RS was generated for each individual in GSE65858 according to the formula. The RS remained significantly associated with poor prognosis (HR = 13.261, 95% CI = 2.136–82.352, p = 0.006). The high-risk group HNSC patients remained significantly shorter OS than the low-risk group HNSC patients in GSE65858 (Figure 6A). The RS also remained an independent prognostic factor compared to routine clinical factors (Figure 6B).

Figure 6

The nine- EMT hallmark gene-based risk score in GSE65858. (A) The Kaplan–Meier curves with the Log rank test of the high- and low-risk groups. (B) The results of the multivariable Cox analysis for the nine-EMT hallmark gene-based risk score (RS) and routine clinical factors. *P < 0.05, and **P < 0.01.

Discussion

Previous studies constructed a prognostic stratification system from multiple perspectives, such as immunity-related gene-based signature,15 microRNA-based signature,16,17 and microenvironment-based system.18 Few studies focused on EMT-related gene-based signature. In our present study, a nine-EMT hallmark gene-based RS was generated to successfully identify the relatively high-risk HNSCs. It was an independent prognostic factor compared to routine clinical factors; moreover, it was validated in an external data set. This may provide more references for clinical decision-making.

Unsurprisingly, some of the nine EMT hallmark genes reported were associated with HNSC or other types of cancer. CAP2 was considered to be related to multistage hepatocarcinogenesis.19 Elevated DKK1 was reported as an independent unfavorable prognostic indicator of survival in HNSC,20 which is consistent with our result. IL-6 plays an important role in HNSC tumor proliferation and metastasis, and IL-6/STAT3 signaling may be a potential target for treating HNSCC patients.21 PLOD2 was found to contribute to drug resistance in laryngeal cancer by promoting cancer stem cell-like characteristics.22 PTX3 is an extrinsic oncosuppressor regulating complement-dependent inflammation in cancer.23 The loss of SFRP1 expression is associated with colorectal cancer, prostate cancer, and renal cell cancer.24,25 Our analysis found that it may also play a tumor suppressor role in HNSC in association with a good prognosis. Gene TGFBR3 was found to play a dual role in bladder cancer, acting as both a tumor suppressor and as a tumor promoter.26 However, a previous study showed it can block lymph node metastasis in head and neck cancer.27 Our analysis shows that TGFBR3 is a protective gene in HNSC. VEGFC may contribute to HNSC growth and motility.28 There are few reports of functional experiments regarding PCOLCE2 in HNSC. In the present study, we found that PCOLCE2 is associated with poor prognosis in HNSC.

Furthermore, we also conducted GSEA to explore the biological characteristics of patients identified as high-risk patients by this nine-EMT hallmark gene-based RS. Based on the results of GSEA, high-risk HNSC is characterized by active EMT program, angiogenesis, MTORC1 signaling, unfold protein response (UPR) program, and high hypoxia. Hypoxia is common in HNSC cells and contributes to malignant behaviors, such as tumor progression, invasion, metastasis, and resistance to chemotherapy and radiotherapy.29 Recent studies suggest that the UPR may affect many hallmarks of cancer, including metastasis, genome stability, angiogenesis, inflammation, and drug resistance.30 These partially explain the reasons for the poor prognosis of the high-risk group. There are numerous ongoing efforts to target the mTOR signaling pathway for cancer therapy,31 however, whether active MTORC1 signaling in high-risk HNSC indicates the response to therapies for targeting mTORC1 signaling still needs further exploration.

Although the present study may provide new insight into the prognostic systems in HNSC, it has several noticed limitations. First, not all EMT hallmark genes were included in the analysis due to the data sets coming from different centers. Thus, the RS may be improved in further study. Secondly, the role of some of these nine genes in HNSC is not yet clear, therefore, it is not clear whether these genes are causal or merely markers for predicting prognosis in HNSC. Thirdly, the present study lacks molecular experimental verification of candidate molecules.

In conclusion, we proposed and validated a nine-EMT hallmark gene-based risk score system for predicting prognosis for patients with HNSC, and also preliminarily explained the biological characteristics of high-risk patients.