{'fraction_non_alphanumeric': 0.0364434856306341, 'fraction_numerical': 0.019801217351105634, 'mean_word_length': 4.89464123524069, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 7, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 system 4 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][7][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.

Data Processing
The RNA sequencing (RNA-seq) data (displayed as raw counts) and clinical information of HNSC in The Cancer Genome Atlas (https://www.cancer.gov/) were downloaded to generate an EMT hallmark gene-based risk score (RS) to predict prognosis. Another HNSC data set GSE65858 10 based on platform of GPL10558 was downloaded from Gene Expression Omnibus (https://www.ncbi. nlm.nih.gov/geo/) and used to validate the RS. The RNAseq data was normalized using voom function from limma package 11 in R software (version 4.0.2) (https://www.r-pro ject.org/). 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 (Table S1) 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.

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 (https://CRAN. R-project.org/package=glmnet) R package to select the optimal prognostic EMT hallmark genes. Thus, the EMT hallmark gene-based risk score (RS) was created as: RS = Expr gene1 *Coef gene1 + Expr gene2 *Coef gene2 + Expr gene2 *Coef gene2 + … 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

4220
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, GSEA 12,14 was performed using the GSEA JAVA program (version 4.0.1) (http://www.broadinstitute.org/gsea). 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 (https://CRAN. R-project.org/package=timeROC). All tests were twosided and p < 0.05, unless otherwise stated, was considered to indicate statistical significance.

Multiple EMT Hallmark Genes Upregulated in HNSC
A total of 93 EMT hallmark genes (Table S2) 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).

4221
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 prognosisassociated 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).

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).

Discussion
Previous studies constructed a prognostic stratification system from multiple perspectives, such as immunityrelated 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 highrisk 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 complementdependent 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

4223
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.

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.

Data Sharing Statement
The data for this study can be obtained from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/ geo/) and in The Cancer Genome Atlas (https://www. cancer.gov/).

Funding
There is no funding to report.

Disclosure
The authors declare that they have no competing interests.