RNA sequencing data sets together with the corresponding clinical characteristics of patients with HNSCC were downloaded from The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/). Manual reannotation of RNA-seq data sets was used to separate expression data into mRNA and lncRNA. The expression levels were transformed as log₂(x+1) and standardized (16). Originally, RNA-seq data sets from 546 patients with HNSCC was collected. Only, 453 cases with complete follow-up data which included survival time ≥ 30 days, clinicopathological profile were included in the follow-up analysis.

The GSE15659 dataset, downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) DataSets, contains five T cell subtypes, including 2 suppressive Treg cell subpopulations, and nonsuppressive Treg cell subpopulations, named as Group 1 and Group 2 in our study, respectively. The dataset was based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array). Specifically, we reannotated the probe set by the affymetrixHgU133Plus2.0 array and ensured that the probes mapped to the genome were unique. We analyzed the expression differences of transcripts including mRNA and lncRNA between groups 1 and 2 according to the standard of | log2FC | > 1 and P value < 0.05. The differential mRNAs in GSE15659 were crossed with mRNAs obtained from the expression profiles of patients with HNSCC from the TCGA database, and the crossed mRNAs were listed as Treg-related mRNAs.

All patients with HNSCC were randomly divided into training and testing cohorts according to their survival status at the ratio of 1:1 (16). In the training cohort, the correlation score between lncRNAs obtained from TCGA data set and Treg-related mRNAs were calculated by Person correlation coefficient test. LncRNAs that correlated with Treg-related mRNAs were selected for subsequent analysis according to the following criteria of | correlation coefficient | > 0.6 and P < 0.001 (17–19). Finally, a Treg-related mRNAs-lncRNAs co-expression network was constructed according to the criteria of | correlation coefficient | > 0.3 and P < 0.001 (20).

To identify survival-related LncRNAs among Treg-related LncRNAs, we performed univariate Cox proportional hazard analysis (21). Treg-related LncRNAs in univariate analysis P<0.01 were included in the Least Absolute Shrinkage and Selection Operator (Lasso) regression. A multivariate Cox regression model was finally used to construct a prognostic signature based on the candidate Treg-related lncRNAs generated from the Lasso regression results (22–24).

A multivariate Cox regression model was used to establish an independent prognostic signature. The risk score for each patient sample was calculated as the expression value of each lncRNA multiplied by the sum of their weights in the multivariable Cox model. The median risk score was used to separate patients into high- and low-risk groups. To validate the predictive value of the model, we performed the Kaplan–Meier log-rank test and time-dependent ROC curve analysis which were used to compare survival between high and low-risk groups in the training, testing, and total cohorts.

Univariate and multivariate Cox regression analyses were performed on the clinical data to determine whether the risk score is an independent indicator of prognosis. ROC curve analysis by univariate and multivariate Cox regression was performed to analyze the correlation between prognosis and clinicopathological factors including age, gender, tumor size (T), lymph node metastasis (N), distant metastasis (M), clinical stage, and risk score. Time-dependent receiver operating characteristic (ROC) curves were plotted to evaluate the accuracy of different clinicopathological factors and risk scores in predicting survival time by using the survival ROC R package (25). Finally, we constructed a prediction model by using nomogram. We then tested the accuracy of the prediction model through 3-year and 5-year calibration curves.

The expression of the 17 Treg-related lncRNAs were correlated with clinicopathological factors using the ggpubr R package. Moreover, we analyzed the correlation of 17 Treg-related lncRNAs the overall HNSCC patient survival.

Principal component analysis (PCA) was used for effective dimensionality reduction, pattern recognition, and exploratory visualization of high-dimensional data of the whole-genome expression profiles retrieved from TCGA, 462 Treg-related coding genes, and 17-lncRNA signature expression profiles. Gene set enrichment analysis of (GSEA, http://www.gsea-msigdb.org/gsea/index.jsp) in high- and low-risk groups. GSEA was used for gene functional annotation and is a powerful analytical method for comparing genes with predefined gene sets obtained from whole genome expression profiles (26). In our experience, GSEA tends to have high repeatability and explanatory power in the analysis of molecular map data. Finally. gene expression matrix data were screened and analyzed by CIBERSORT (https://cibersortx.stanford.edu/) (27). Specifically, immune cell populations of infiltrating T cells in high- and low-risk groups were compared to access the relationship between 17-lncRNA signature and immune cell infiltration.

The flowchart in Figure 1 describes the order of computational steps we undertook to identify 17 lncRNAs. First, we obtained 18,392 mRNA and 9,357 lncRNA expression profiles which corresponded to clinical data from 453 patients registered in TCGA database. In the GSE15659 dataset from GEO, a total of 648 differentially expressed transcripts were obtained ( Table S1 ). The data was divided into 281 up-regulated and 367 down-regulated transcripts and visualized using Volcano plot ( Figure 2A ). The top 50 differentially expressed transcripts were graphed as heat map using the heatmap R package ( Figure 2B ). 462 mRNA present in both GSE15659 and TCGA datasets were identified using Venn Software Analysis. The resulting overlapping mRNA were named Treg-related mRNAs for HNSCC ( Figure 2C ). Data enrichment was used to gather functional information (GO and KEGG) utilizing the Database for Annotation, visualization and integrated discovery (DAVID). Biological pathways related to Treg cell expression, included insulin secretion pathway (28, 29), pancreatic secretion pathway, Hippo signaling pathway (30–32) and others ( Figure S1A and Table S2 ). In addition, the Treg-related mRNA transcripts related to microglial cell activation, cell adhesion and cell-cell communication were detected ( Figure S1B and Table S2 ).

The total cohort of 453 patients with HNSCC was randomly divided into training and testing cohorts at ratio of 1:1. In the training cohort, 1652 lncRNAs correlated with Treg-related mRNAs in patients with HNSCC ( Table S3 ) and were designated as Treg-related lncRNAs. Following that, 35 lncRNAs were identified by Univariate Cox regression analysis ( Table S4 ). Lasso regression analysis was carried out on 35 prognostic lncRNAs to improve confidence in the prediction. We identified 17 lncRNAs prognostic factors for patients with HNSCC ( Figures 3A, B ) (22). The correlation between Treg-related mRNAs and the associated lncRNAs is shown in Figure S2 and Table S5 . Of 17 assigned lncRNAs, additional multivariate Cox regression analysis of defined transcripts LINC00460 and AC092115.3 as significant independent prognostic factors ( Figure 3C ). These lncRNAs were then used to establish lncRNA prognostic signature. Six lncRNAs, namely, CAVIN2-AS1, AC007878.1, LINC00460, AC092115.3, AC068446.2, and LINC01976, were used to calculate risk factors for the prognosis of patients with HNSCC using the cutoff of HR>1. Meanwhile, 11 lncRNAs (AL157414.1, LINC01281, GLYCTK-AS1, LINC02325, AC026362.1, AL049552.1, STARD4-AS1, AC103809.1, AC104083.1, AC004461.2, and LINC02202) were protective factors for the prognosis of patients with HNSCC with the cutoff of HR<1 ( Figure 3C and Table S6 ). The risk score for prognosis was calculated as Risk score = Σ i E x p i ( l n c R N A i ) ∗ C o e f ( l n c R N A i ) where Expi is the expression value of each lncRNA, and Coef is the regression coefficient of the multivariate Cox analysis for the target lncRNA (21, 23).

Median risk score for 17-lncRNA obtained in prognosis model was used to divide patients with HNSCC into high- and low-risk groups. We first used scatter plot ( Figures 4A–C ) and risk curve ( Figures 4D–F ) to describe the risk score and survival status of each patient with HNSCC, in the training, testing and total cohorts. The risk coefficient and mortality rate in the low-risk group were lower than those in the high-risk group. The observed mortality rate correlated with the risk score ( Figures 4A–F ). The heatmap showed that CAVIN2-AS1, AC068446.2, LINC01976, AL157414.1, LINC00460, and AC092115.3 were highly expressed in the high-risk group, while AC007878.1, LINC01281, GLYCTK-AS1, LINC02325, AC026362.1, AL049552.1, STARD4-AS1, AC103809.1, AC104083.1, AC004461.2, and LINC02202 were highly expressed in the low-risk group ( Figures 4G–I ). To assess further the accuracy of 17-lncRNA prognostic signature, we constructed the K-M survival curve using the R package “survival”. In the training cohort, the overall survival (OS) of the low-risk group was better than that of the high-risk group (P < 0.001, Figure 4J ). These results in the testing cohort (P < 0.001, Figure 4K ) and total cohort (P < 0.001, Figure 4L ) were consistent with the results of the training cohort. The risk score correlated with the OS time of patients with HNSCC and hence is a good predictive value for HNSCC prognosis. To further evaluate model quality, we calculated area under the Curve (AUC) for the 3, 5 and 8-year survival curves using the time ROC R package. The AUC values were 0.829, 0.766, and 0.758 in the 3-, 5-, and 8-year follow-ups of the training cohort ( Figure 4M ) and 0.615, 0.578, and 0.58 in the testing cohort ( Figure 4N ). In the total cohort, the values were 0.721, 0.679, and 0.668, respectively ( Figure 4O ). Hence, the 17-lncRNA signature is a reliable measure of prognostic signature of HNSCC.

To further evaluated the risk score as an independent prognostic marker for patients with HNSCC, we performed univariate COX regression analysis comparing the risk score to other clinicopathological factors (age, gender, tumor size (T), lymph node metastasis (N), distant metastasis (M), clinical stage) using the survival ROC R package ( Figures 5A–C ). Age, gender, M, and risk score were independent prognostic indicators for patients with HNSCC. Through the multivariate COX regression analysis, we found that M and risk score were correlated with OS in the training cohort ( Figure 5D ). The multivariate COX regression analysis result of ‘M’ for OS showed p values <0.05 in the training, testing ( Figure 5E ) and total cohorts ( Figure 5F ), indicating that M was a significantly independent prognostic factor for HNSCC. In the same way, ‘risk score’ was also a significantly independent prognostic factor for HNSCC in the training, testing, and total cohorts. The multivariate COX regression analysis suggested that the risk score containing 17-lncRNA signature was a significant prognostic factor for HNSCC, independent of clinicopathological parameters. The sensitivity and specificity of the risk score in predicting the prognosis of patients with HNSCC was investigated by comparing the area under the ROC curve ( Figure 5G ) which measured changes in the risk score and other clinicopathological factors in predicting the overall survival of HNSCC patients. The AUC values of the risk score, age and M were 0.693, 0.536, and 0.506, respectively. The risk score showed a better AUC than other clinicopathological factors, indicating that risk score is more effective in predicting HNSCC prognosis.

In attempt to develop a clinically applicable method for predicting the survival probability of patients, we generated a nomogram using the rms R package ( Figure 5H ), plotting risk score, age, and M. The risk score had the biggest contribution to the 3- and 5-year OS of patients with HNSCC. In addition, we supplemented our model with the 3-year and 5-year calibration charts. The 3-and 5-year OS calibration curves were well predicted compared with the ideal models in all cohorts ( Figure 5I ). The results showed that the nomogram can independently evaluate survival of patients with HNSCC, which may help doctors to make better medical decisions and follow-up plans.

The 17-lncRNA signature has been revealed effective for predicting HNSCC prognosis in aforementioned algorithm. However, whether each lncRNA has correlation with clinicopathological factors still remains unclear and needs to be further assessed separately ( Figure S3 ). Eleven lncRNAs were associated with disease progression, and included AL157414.1, LINC01281, CAVIN2-AS1, GLYCTK-AS1, LINC02325, AC026362.1, AC007878.1, LINC00460, AL049552.1, AC092115.3, STARD4-AS1, AC103809.1, AC104083.1, AC004461.2, AC068446.2, and LINC01976. Furthermore, three lncRNAs, namely, AC092115, GLYCTK-AS1, and LINC02202, were associated with different TNM classification of HNSCC. The expression of six lncRNAs such as LINC02325, AC026362.1, AL049552.1, AC103809.1, LINC01976, and LINC02202 was statistically significant between different genders, while no statistically significant difference was observed for 17 lncRNAs among different age groups ( Figure S4 ). Kaplan-Meier curves reveal that high expression of LINC00460, AC092115.3 and low expression of LINC01281, LINC02325, AC026362.1, AC007878.1, AL049552.1, STARD4-AS1, AC104083.1, AC004461.2 was associated with poor survival (P<0.05). Only ten of the 17 individual lncRNAs had a survival prediction power for HNSCC.

Principal Component Analysis (PCA) was performed to evaluate expression differences among the low- and high-risk groups. Whole-genome expression profiles from the TCGA and 462 Treg-related genes was insufficient to separated high and low risk groups ( Figures 6A, B ). However, the expression differences in 17-lncRNA signature could obviously distinguish high-risk patients from low-risk patients ( Figure 6C ). The GSEA analysis confirmed that mRNAs related 17-lncRNA detected in the low-risk groups were enriched for transcripts related to immune-related biological processes- immune system development, leukocyte mediated immunity, and regulation of immune system process ( Figures 6D–F ). and in addition, the low-risk group had a better overall immune response than the high-risk group. The infiltrating immune cells in the HNSCC microenvironment were analyzed by CIBERSORT to explore the differences of immune cell infiltration between high-and low-risk group. Genes related to the naive B cells, CD8 T cells, memory activated CD4 T cells, follicular helper T cells, and regulatory T cells (Tregs) were found to be expressed at higher levels in low-risk group rather than in high-risk group ( Figure 6G ). In addition, expression of naive T cells, memory resting CD4 T cells, macrophage M0, and activated mast cells was significantly lower in low-risk group than in high-risk group ( Figure 6G ). These results suggest that 17-lncRNA signature model is capable of distinguishing tumor microenvironment in the low and high-risk groups of HNSCC patient population.