Flowchart of the study protocol of TGF-β-related characteristics related to the prognosis of HNSCC is listed in Figure 1. We used open datasets from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/) and Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) databases, including the GSE65858, GSE75538 and GSE117973 chip datasets. The TGF-β mRNA expression status and correlating tumor immune microenvironment evaluation indicators in HNSCC and corresponding normal tissues were analyzed through the SangerBox database (http://sangerbox.com/Tool). We identified the p value as 0.001, with a fold change of 1.5.

RNA-seq data from TCGA HNSCC and GEO HNSCC with unreliable and incomplete clinical data were removed. The clinical statistical information of HNSCC patients after data pretreatment is shown in Table 1.

To construct the TGF-β associated prognostic model, we downloaded HNSCC transcriptional data from public database. All selected data were screened and selected with complete gene expression information and clinical information. We applied 54 TGF-β -related genes, univariate survival analysis and Lasso Cox multivariate analysis to determine the risk score of seven key genes. High and low risk groups are divided by selecting the best threshold. Therefore, we determined that λ = 0.0086 would obtain the optimal model. In addition, we selected seven genes at λ = 0.0086 as targets for the next step. The final 7-gene signature formula is as follows: risk score = 0.0003 ∗ BCAR3-0.062 ∗ ID2 + 0.112 * NOG +0.009∗ SERPINE1 + 0.177 SLC20A1 + 0.106 ∗ THBS1-0.351 ∗ TR IM33.

The HNSCC RNA-seq data from TCGA were identified as the training set. We adopted univariate Cox proportional risk regression to construct and predict a TGF-β-associated gene model (54 in total). Overall survival data were analyzed by using the R package survival COXPH function. We defined the threshold for filtering as p < 0.05. In this study, we selected five genes with hazard ratios (HRs) larger than one and two genes with HRs less than 1.

Univariate Cox regression analyses of the correlation between TGF-β-associated gene expression and HNSCC clinical prognosis information were carried out by Cox proportional hazards regression analysis. Genes with p < 0.05 were considered significant. Determination of genes in the TGF-β-specific module closely related to prognosis in HNSCC were explored.

We utilized the R software package glmnet for least absolute shrinkage and selection operator (LASSO) Cox regression analysis to determine the greatest impact on the prognosis of HNSCC. LASSO applies an L1-regularization penalty, ρ, to estimate a penalized precision matrix to illustrate indicators with the highest contribution (Love et al., 2016). In the LASSO model, the minimum criterion (λ) based on 10-fold cross validations and 1,000 iterations is chosen. The selected genes were then included in a multivariate Cox regression model, and those gene sets with the best prognostic value were identified by positive selection and reverse elimination methods.

To determine correlations between TGF-β-associated gene expression and HNSCC prognosis, we applied the Wilcoxon test and divided HNSCC patients into two subgroups based on TGF-β expression level. The optimal cutoff point for gene expression was obtained based on the R package “survminer” (cutoff = −0.007). The threshold is −0.007, in which the groups with a risk score greater than −0.007 were high risk groups, and those with a risk score less than −0.007 was identified as low risk groups We used Kaplan-Meier curves to evaluate the prognostic value of various clinical features through the R package ‘survminer’ (CRAN.R-project.org/package = survminer). This method has been described in a previous study (Alcala et al., 2019; Zhuang et al., 2020).

To evaluate somatic mutations in HNSCC, we applied the package TCGA biolinks in R and downloaded mutation annotation files. Somatic single-nucleotide polymorphisms and indels in tumor samples were called using the MuTect2 (http://www.broadinstitute.org/cancer/cga/mutect) pipeline Genome Analysis Toolkit (GATK; Broad Institute, Cambridge, MA, United States). We counted differences in the number of mutant genes in the samples. Furthermore, we screened out genes with mutation frequencies greater than three and used the chi-square test to screen genes with significantly high-frequency mutations in each subtype, with a selection threshold of p < 0.05.

GSEA was performed to examine different biological processes with GSEA software (http://software.broadinstitute.org/gsea). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used for hierarchical analysis for high- and low-risk groups of HNSCC patients. The R package clusterProfiler (https://guangchuangyu.github.io/software/clusterProfiler) (v3.14.0) was employed to process the KEGG and GO analyses. We identified the number of random sample permutations as 1,000, and enriched gene sets with a nominal p < 0.05 and 25% cutoff on false discovery rate (FDR) were defined as significant.

We applied the “ESTIMATE” R package to assess overall immune infiltration indicators based on the medium scores of ImmuneScore, StromalScore, and ESTIMATEScore (Ghatalia et al., 2019). We also calculated the proliferation score of all cells. These proliferation scores were obtained from previous study (Thorsson et al., 2018).

The TIDE algorithm was used to link individual immunotherapy responses with the TIDE web tool (Netherlands Cancer Institute, Amsterdam, Netherlands, available from http://shinyapps.datacurators.nl/tide/). TIDE is used to estimate the spectrum and frequency of small insertions and deletions (indels) generated in a pool of cells by genome editing tools such as CRISPR/Cas9, TALENs and ZFNs. In this study, we applied TIDE to estimate the likelihood of immunotherapy response.

The human nasopharyngeal carcinoma cell NPC and human immortalized nasopharyngeal epithelial cell NP69 were purchased from FuHeng (Shanghai, China). The cell line NP69 was cultured in KM medium, and the cell line NPC was cultured in Dulbecco’s modified Eagle medium supplemented with 10% fetal bovine serum (ThermoFisher Scientific, Waltham, MA, United States). Cell lines were grown at 37 °C in a humidified incubator containing 5% CO₂.

Total RNA was extracted from NP69 and NPC cells using TRIzol reagent (Invitrogen Life Technologies, Waltham MA, United States), followed by reverse transcription according to the manufacturer’s instructions (Takara, Japan). The specific quantitative primers used are listed in Supplementary Table S1. Samples were assessed by quantitative real-time reverse transcriptase–polymerase chain reaction (qRT-PCR) using an Agilent Mx3005P using SYBR qPCR Mix (MQ10201s, Monad Biotech, Wuhan, China). Human glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an endogenous control. Relative expression levels were defined according to the 2^−ΔΔCt method. Each experiment was performed in triplicate.

We calculated correlations with the Pearson correlation coefficient, and differences between subgroups were determined by the Wilcoxon test or Kruskal–Wallis test. The data are expressed as means ± SD. A two-tailed p-value less than 0.05 was defined as statistically significant.

The workflow of this study is depicted in Figure 1. To select proper protective factors and risk factors, we applied TCGA-HNSCC datasets as a training set; GSE65858, GSE75538, and GSE102995 were used as external independent validation datasets. These datasets were obtained from GEO database, which were used for verified the accuracy of the model.

We adopted LASSO Cox regression to obtain the change trajectory of each independent variable, as illustrated in Figure 2A. The expression levels of the identified seven genes in human nasopharyngeal carcinoma cell NPC and human immortalized nasopharyngeal epithelial cell NP69 are shown in Figure 2B. We applied 10-fold cross-validation techniques to avoid performance bias with all prediction methods (Figure 2C). We also calculated the risk score for each sample in the dataset TCGA-HNSCC according to the TGF-β expression level and plotted the risk score distribution of the samples, as shown in Figure 2D. The different expression of seven different signature genes with the increase of risk value was assessed. High expression of BCAR3, NOG, SERPINE1, SLC20A1, and THBS1 was identified as a risk factor associated with a high risk score. Conversely, high expression of ID2 and TRIM33 was associated with low risk, constituting a protective factor (Figure 2D). Furthermore, we used the R software package and receiver operating characteristic (ROC) curve analysis to evaluate prognostic factors. We analyzed the classification efficiency for prognosis prediction, and areas under the curve (AUCs) at 1, 3 and 5 years were 0.64, 0.65, and 0.55, respectively, as indicated in Figure 3C. An AUC = 0.64 at 1 year was the most reliable indicator for survival prediction with this 7-gene model (Figure 2E). To further investigate clinical prognosis between the high- and low-risk score subgroups, we applied the R package “survminer” with the cutoff = -0.007. The Kaplan-Meier curve indicated that the difference between the high-and low-risk score subgroups was highly significant (p < 0.0001) (Figure 2F). In total, 318 samples were assigned to the high-risk group and 199 to the low-risk group. We downloaded validation datasets such as GSE65858, GSE75538, and GSE102995 and applied the same gene models and coefficients. As the GSE65858 dataset lacked NOG gene expression, we used a six-gene model in further validation. The results of the three external validation datasets were highly consistent with those of the training sets. Overall, the high-risk score subgroup had a poorer prognosis than the low-risk score subgroup (Figures 2G–I). The forest plot showed that TGF-β signaling pathway-associated gene-based characteristics might be major risk factors for HNSCC (HR = 2.069, 95%CI = 1.582–2.698, p < 0.001) (Figure 2I).

To gain insight into the mutational mechanisms of HNSCC, we calculated the tumor mutation burden (TMB) as the number of somatic mutations for each patient. However, patients in the high- and low-risk score subgroups showed no significant difference in TMB (p = 0.07) (Figure 3A) or number of mutated genes (p = 0.067) (Figure 3B). Given a set of mutational signatures, we calculated the presence of somatic mutations to reveal signatures for these five genes (TP53, CDKN2A, SYNE1, PKHD1L1, and PEG3) and found nonsense mutations and missense mutations to be the most common mutations. In addition, the most frequent genetic lesion in HNSCC was in TP53, followed by CDKN2A, SYNE1, PKHD1L1, and PEG3 (Figure 3C).

In further investigating clinical signatures between the high- and low-risk score subgroups, we determined that there was no significant difference in T stage (p = 0.62), N stage (p = 0.84), M stage (p = 0.39), sex (p = 0.62), or age (p = 0.98) (Supplementary Figure 1A-F). In contrast, a significant difference in grade was observed (p = 0.00045) (Supplementary Figure 1G). In addition, we further analyzed overall survival as related to the clinical signature between these two subgroups. The Kaplan mire curve revealed a significant difference for male sex (p < 0.0001) ( Figure 4A), age over 65 years (p = 0.0023) (Figure 4B), younger and equal to 65 years old (p = 0.00044) (Figure 4C), T1-T2 stags (p = 0.031), T3-T4 stages (p < 0.0001) (Figure 4E), N1 stage (p = 0.0078) (Figure 4F), M0 stage (p < 0.0001) (Figure 4G), M1-M2 stages (p = 0.037) (Figure 4H), grade 1–2 (p = 0.0053) (Figure 4I), and grade 3–4 (p = 0.00083) (Figure 4J). Overall survival analysis showed no significant difference between the high- and low-risk score subgroups for female sex e and N0 stage (Figure 4L).

To assess the independence of the clinical signature as a prognostic factor, we applied multivariate analyses and the Cox proportional hazards regression model and determined that clinical characteristics such as age (p = 0.022, HR = 1.4 (1.05.1.88)), stage (p = 0.002, HR = 2.16 (1.34.3.5)) (Figure 4M and N), and risk score type were significantly associated with prognosis. Overall, the seven-gene signature model has good predictive value for HNSCC.

When analyzing signaling pathway enrichment, we found that the high-risk score subgroup was enriched in hypoxia and EMT pathways (Figure 5A). To deeply explore prognosis-associated biological pathways, we obtained 13 positively correlating genes. The heatmap in Figure 5B illustrated that the expression level of these genes in the high- and low-risk score subgroups. GO classification includes biological process (GO-BP), cell component (GO-CC), and molecular function (GO-MF) categories. We identified 10 pathways as enriched in BP annotations (Figure 5C), seven in CC annotations (Figure 5D) and three in MF annotations (Figure 5E). KEGG analysis showed six biological signaling pathways to be involved in HNSCC regulation (Figure 5F).

To investigate the relationship between risk score and immune cell infiltration and inflammatory characteristics in patients with HNSCC, we evaluated three immune microenvironment scores, namely, the stromal score, immune score, and ESTIMATE score, and found that the stromal score (Figure 6A) and ESTIMATE score (Figure 6B) were significantly higher in the high-risk group than in the low-risk group. Conversely, no significant difference was found for the immune score (Figure 6C). Moreover, immune cell infiltration evaluations indicated significant differences in the proportions of 11 types of immune cells in the subgroups (Figures 6D,E). Among them, the proportions of resting memory CD4⁺ T cells, resting NK cells, M0 macrophages and activated mast cells in the low-risk group were significantly lower than those in the high-risk group, whereas proportions of natural B cells, memory B cells, CD8⁺ T cells, follicular helper T cells, activated NK cells, regulatory T cells and resting mast cells were significantly higher in the low-risk group.

To probe inflammatory activity associated with risk scores, we examined relationships between seven metagene clusters, whereby differences represent different inflammatory and immune responses. The characteristics based on the TGF-β pathway were explored, and detailed information on these metagenes is provided in Figure 6F.

To verify the gene expression details observed, gene set variation analysis (GSVA) was utilized to calculate scores for the corresponding clusters of seven metagene clusters. Our results showed that the risk score correlated positively with IgG, LCK and MHC-II and with IgG, interferon and MHC-I (Figure 6G). At the same time, we compared differences of these seven scores in the high- and low-risk subgroups and found that IgG and LCK scores were significantly higher in the latter but that interferon and MHC-I scores were significantly higher in the former (Figure 6H).

To evaluate potential differences in clinical effects between immunotherapy and chemotherapy, we adopted TIDE software and found that the high-risk score subgroup had a markedly higher TGF-β response (p = 3.4e-23) (Figure 7A), proliferation score (p = 0.038) (Figure 7B), wound healing score (p = 0.016) (Figure 7C), and exclusion score (p = 1.2 e-21) (Figure 7D) compared with low-risk score subgroup. However, macrophage regulation (Figure 7E) and dysfunction (Figure 7F) showed no significant difference between the subgroups. The TIDE prediction therapy survival curve suggested that the FALSE group had a better prognosis (p = 0.021) (Figure 7G), and the TIDE score was much higher in the high-risk score group than in the low-risk score group (p = 5.7e-14) (Figure 7H).

We further analyzed differences in immunotherapy and chemotherapy among different subtypes of immune molecules. Subclass mapping was used to compare the similarity between the high-low risk subgroups in our TCGA-HNSC dataset and immunotherapy patients in the GSE78220 dataset: the lower the p value, the higher the similarity. For TCGA-HNSC, the low-risk score group was more sensitive to PD1 treatment (Figure 7I). This result was consistent with the TIDE results, indicating that low-risk group patients will benefit more from immunotherapy. The results also showed that the high-risk score subgroup was more sensitive to the traditional chemotherapy drugs cisplatin (p = 0.69e-10) (Figure 7J), erlotinib (p = 0.0017) (Figure 7K), paclitaxel (p = 3.4e-11) (Figure 7L) and crizotinib (p = 8.3e-08) (Figure 7M). Overall, this risk score model provides a novel basis for HNSCC patient treatment options.