Identification of a nomogram based on an 8-lncRNA signature as a novel diagnostic biomarker for head and neck squamous cell carcinoma

Preprocessing of the data sets

We downloaded the gene matrix of 546 samples from the TCGA-HNSC database, which included 502 tumour and 44 normal samples. We divided all HNSCC patients with complete information (n=499) in TCGA-HNSC into training cohort and validation cohort, in a random manner according to a ratio of 3:2.

Moreover, From May 2017 to August 2018, a total of 102 frozen, surgically resected tumor tissues were obtained from patients with pathological diagnosis of HNSCC at Chengdu Third People's Hospital. The specimens were frozen with liquid nitrogen immediately after removal and transferred to the −80°C refrigerator.

Differential analysis

We conducted a differential analysis of the 300 tumor and 44 normal samples. Eventually, we obtained a total of 19754 mRNAs and 14847 lncRNAs. After obtaining the expression data, we identified differentially expressed genes using the software package EdgeR, selecting genes that had at least 2-fold higher expression levels in HNSCC samples (Poisson model FDR < 0.05). Therefore, after screening, we obtained 4150 reliably expressed mRNAs and 253 lncRNAs (Figure 1A, 1B).

Identification of 8-lncRNAs for predicting HNSCC patient survival

A total of 253 lncRNAs with significant differences were identified to have prognostic significance in univariate Cox survival analysis, and 41 with P < 0.05 were screened out and applied in the following analysis (Figure 2A). As shown in Figure 2B, 2C, LASSO regression analysis identified 11 lncRNAs (lambda value=11), which were then used in the multivariate Cox regression. Finally, 8 lncRNAs for predicting HNSCC patient survival were identified, including MIR4435-2HG, LINC02541, MIR9-3HG, AC104083.1, AC099850.4, PTOV1-AS2, AC245041.2, and AL357033.4.

The role of the 8-lncRNA signature in HNSCC’s prognosis

The signature score of these 8 lncRNAs based on regression coefficients in multivariable Cox analysis was calculated as follows: signature score = (0.36314 × expression of MIR4435-2HG) + (0.23003 × expression of LINC02541)– (0.22031 × expression of MIR9-3HG) – (0.23426 × expression of AC104083.1) + (0.21344 × expression of AC099850.4) – (0.27806 × expression of PTOV1-AS2) + (0.25463 × expression of AC245041.2) – (0.31513 × expression of AL357033.4). Taking the median signature score as the dividing point, the patients were divided into high signature-score group and low-signature score group. (Figure 2D). Patients in the high-signature score group had a significantly worse OS than those in the low-signature score group (Figure 2E). Besides, the AUCs were assessed for 3years (AUC = 0.740) and 5years (AUC = 0.706) survival (Figure 2F), and the results suggest that the signature can effectively evaluate the prognosis of HNSCC patients.

Development of a prediction model integrating the 8-lncRNA signature and clinical characteristics

We evaluated age, sex, lymph node (N) status, metastasis (M) status, tumor stage (stage), and new events (which include locoregional disease, locoregional recurrence, new primary tumor, and distant metastasis) using KM analysis. Next, we found that age, metastasis, and new event play an important role in the prognosis of HNSCC (Figure 3).

The signature was regarded as a predictor for HNSCC patients. We identified the significant variables through univariate Cox analysis. The multivariate model includes candidate variables with a P-value < 0.1 in univariate analysis. (Figure 4A). Finally, the results (Table 1) suggested that the independent risk factors for HNSCC, including: stage, M stage, new event, and signature score. Moreover, we compared the multivariate Cox regression results of the two groups with and without the signature score. Surprisingly, the C-index of the signature score-containing group (0.72) was higher than that of the signature score-free group (0.71) (Supplementary Figure 1). The nomogram model was built by using the coefficients of the multivariable Cox regression model (Figure 4B). The AUC for 3-year survival reached 0.788 (Figure 4C). What’s more, the calibration curve shows that concerning the probabilities of 3-year OS and 5-year OS, the predicted values are consistent with the observed values (Figure 4D). Finally, we calculated the total risk score based on each predictor in the nomogram model. Kaplan-Meier analysis showed that patients in the high-risk group had a significantly worse OS than those in the low-risk group (Figure 4E).

Validate the signature in the internal and external validation cohorts

To determine the stability of this nomogram; we performed a similar analysis process in the validation cohort (n = 199). Taking the median signature score as the dividing point, the patients were divided into the high signature-score group (n = 100) and the low signature-score group (n = 99). with the median signature score as the cut-off point (Figure 5A). The Kaplan-Meier OS curves suggested that patients in the high-signature score group had a significantly worse OS than those in the low-signature score group (Figure 5B). The AUC value for 3-year survival exhibited by the 8-lncRNA signature reached 0.779 (Figure 5C). Besides, the calibration curve shows that concerning the probabilities of 3-year OS and 5-year OS, the predicted values are consistent with the observed values (Figure 5E). What’s more, using the same total risk score formula in the internal validation cohort, the Kaplan-Meier OS curves showed that the OS of patients with the high-risk score was significantly worse than that of patients with the low-risk score (Figure 5F). The AUC exhibited by the total risk score for 3-year survival reached 0.796 (Figure 5E).

We also validated the robustness of the signature in GSE65858 (n = 270), which had an AUC of 0.785 for 3-year OS (Figure 6A, 6C). Moreover, the OS of patients with high-signature score was worse than those of patients with the low-signature score (Figure 6B). The Kaplan-Meier OS curves manifested that patients in the high total risk score group had a significantly worse OS than patients in the low total risk score group (Figure 6F). Similarly, the calibration curve showed good agreement between the predicted and observed values (Figure 6E), and the AUC exhibited by the total risk score for 3-year survival reached 0.811 (Figure 6D).

Furthermore, we measured the expression of these eight lncRNAs in 102 HNSCC samples by qRT-PCR (Figure 7A). The Kaplan-Meier curve showed that the OS of the patients with a high-signature score was significantly worse than that of the patients with a low-signature score (Figure 7B). The AUC for 3-year survival reached 0.942 (Figure 7C). The Kaplan-Meier curve showed that the OS of the patients with a high-risk score was significantly worse than that of the patients with a low-risk score (Figure 7E). The calibration curve performs well (Figure 7D), and the AUC exhibited by the total risk score for 3-year survival reached 0.896 (Figure 7F).

Moreover, the above verification process was also performed for the entire TCGA-HNSC set (n=499) and revealed good results (Supplementary Figure 2).

WGCNA

The gene co-expression system was established by WGCNA to screen the biologically significant gene modules related to the lncRNAs in the signature. To create a scale-free system, we set the soft threshold beta to 3(Figure 8A). Besides, genes with similar patterns were clustered in different modules (Figure 8B). The minimum cluster size was determined to be 30 per module. The gene modulus was determined by the dynamic shearing method. The module eigengene (ME) was calculated to explore the similarity of all modules (Figure 8C). Eigengenes were calculated to be correlated with clinical factors. Finally, a robust correlation between the gene significance and grade and signature score was identified (Figure 8D). The ten modules were clustered into two groups (Figure 8E). In order to evaluate the correlation between gene expression and survival time, we calculated the gene significance (Figure 9A). Then, we found that there was a strong correlation between the module members of the brown module and the genetic significance of OS. (cor-value = -0.47, P = 5.3e − 12). The red module, whose hub gene contains MIR4435-2HG, was also negatively correlated with the OS (cor-value = -0.2, P = 0.032) (Figure 9B). Finally, we explore the GO term and KEGG pathway through functional enrichment analysis. (Figure 9C–9F). The results indicated that the biological processes (BP) of these genes mainly involved cell chemotaxis, leukocyte migration, immune response, cell-cell signaling, and so on. The results suggested that the molecular functions (MF) of these genes were related to actin binding, chemokine activity, chemokine receptor binding, ATPase binding, and so on. The results showed that the cellular components (CC) included collagen-containing extracellular matrix, plasma lipoprotein particle, growth cone and site of polarized growth. KEGG pathway functional enrichment showed that leukocyte transendothelial migration, cytokine-cytokine receptor interaction, cell adhesion molecules (CAMs), and the chemokine signaling pathway were mainly related to the genes in these modules.

We conducted a similar analysis process to estimate the correlation between gene expression and grade (Figure 10A). A strong correlation was found between the gene significance for grade and module membership in the turquoise module (which contains MIR9-3HG, AC099850.4 and PTOV1-AS2) (cor-value = 0.41, P = 7.2e − 23); the black module (which contains LINC02541) (cor-value = 0.35, P = 0.00047) and the red module (cor = 0.28, P = 0.0024) were both positively correlated with grade (Figure 10B). We constructed the lncRNA-mRNA network (weight>0.1) diagram of the hub lncRNAs in the turquoise module (Figure 10C). We also carried out functional enrichment analysis to explore the GO term and KEGG pathway (Figure 10D–10G). The results indicated that BP mainly involved cell proliferation, cell division, positive regulation of cell migration, and regulation of the cell cycle. The results showed that MF was related to catalytic activity, acting on DNA, protein binding, and DNA replication origin binding. The results showed that CC included proteinaceous extracellular matrix, chromosome, centromeric region, and extracellular matrix. Moreover, KEGG pathway functional enrichment showed that the cell cycle, the p53 signaling pathway, Cellular senescence, Mismatch repair, and DNA replication were mainly involved.

Data acquisition

The RNA-sequencing data of HNSCC patients were acquired from The Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov/) and The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) [21]. GSE65858 from GEO was conducted on the GPL10558 platform. Besides, we also followed 102 HNSCC patients in the Pathology Department and the Otolaryngology Department of Chengdu Third People's Hospital. The clinical features of patients with HNSCC are presented in Table 2.

Differential analysis

The edgeR package in R software [22] were used to analyze the differentially expressed RNAs in HNSCC and adjacent normal tissues of the TCGA. Significantly expressed RNAs were identified by setting adjusted P values < 0.05 and |log2FC (fold change) | > 1 (|log2FC > 1| and the adjusted FDR < .05) [23, 24].

The construction of the lncRNA-based prognostic signature

The prognostic value of 253 differentially expressed lncRNAs was first calculated in the univariate Cox analysis, and 41 lncRNAs with P < 0.05 were identified as seed lncRNAs for LASSO regression analysis, which identified 11 lncRNAs (R ‘glmnet’, ‘survival’ packages). To determine the prognostic value of the lncRNAs, multivariate Cox regression was further performed using the R survival package based on each “significant” lncRNA identified in the above steps. A lncRNA with P < 0.05 was defined as significant. The corresponding hazard ratios (HRs), 95% confidence intervals (CIs), and P-values were calculated.

Prognostic evaluation using the 8-lncRNA signature

The signature score for each patient in the training group is calculated based on the formula (signature score = expGene1 ×βGene1 + expGene2 × βGene2 + expGenen × βGenen (where exp is the prognostic gene expression level and β represents the multivariate Cox regression model regression coefficient)). All samples are randomly divided into high- and low- signature score sets, with the median signature scores as the cut-off value [25]. The survival analysis of each group was evaluated through the Kaplan-Mayer curve and the log-rank test. Receiver operating characteristic (ROC) curve analysis was employed to assess the specificity and sensitivity of the survival predictions according to the lncRNA signature scores (R package “survivalROC”). A P-value <.05 was considered significant.

Development of a prediction model based on the 8-lncRNA signature and clinical characteristics

The gene signature score as a predictor for HNSCC patients was analyzed in the model. We determined the significant variables through univariate Cox regression analysis. The multivariate model includes candidate variables with a P-value < 0.1 on univariate analysis. Finally, the multivariable Cox regression model began with the clinical candidate predictors as follows: stage, M stage, new event, and signature score. The nomogram model was built with the coefficients of the multivariable Cox regression model (using the R packages “rms”, “Hmisc”, “lattice”, “Formula”, and “foreign”). Then, we calculated the total risk score based on each predictor in the nomogram model and divided the HNSCC patients in the training and internal validation sets into two groups with the median risk score as the cut-off point. Kaplan-Meier curves and the log-rank test were used to compare the survival outcomes of the two groups. Receiver operating characteristic (ROC) curve analysis was employed to assess the accuracy and precision. of the survival predictions according to the total risk scores. Calibration curves were plotted to assess the calibration of the nomogram (R package “rms”). To quantify the discrimination performance of the nomogram, Harrell’s C-index was measured. A P-value <.05 was considered significant.

Validation of the 8-lncRNA signature

The same risk formula was used to validate the internal validation set TCGA-HNSC (n = 199), the entire set TCGA-HNSC (n = 499), the external validation set GSE65858 (n = 270) and the qRT-PCR set (n=102).

Real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR)

Total RNA was reverse-transcribed into cDNA with random primers using the Transcriptor First Strand cDNA Synthesis Kit (Roche, Penzberg, Germany) following the manufacturer’s instructions. The expression levels of the 8 lncRNAs were measured by qRT-PCR using FastStart Essential DNA Green Master mix (Roche, Penzberg, Germany) on a Roche LightCycler 480 (Roche, Penzberg, Germany). Relative expression was determined using inter-experiment normalization to GAPDH. All quantitative PCRs were conducted in triplicate. Divergent primers, rather than the more commonly used convergent primers, were designed for the lncRNAs. Primer specificity was verified using BLAST, with a single peak in the melting curve indicating the generation of a specific product. Three experimental replicates were performed for each sample. Primers used in the study were presented in Supplementary Table 1.

Construction of a weighted gene coexpression network

The procedure of WGCNA [26] included identifying the gene expression similarity matrix, adjacency matrix, and co-expression network. We set the cut-off as a Person correlation coefficient > 0.9 and P < 0.001 to screen gene coexpression with lncRNAs. Then, differentially expressed gene (DEG) analysis was performed among these genes, and we used the expression matrix composed of 4150 differential genes and the above 8 lncRNAs as input files. The power value of the adjacent matrix soft threshold is determined to be 9 to meet the scale-free topology standard. Hierarchical clustering analysis based on average linkage used the dynamic tree cut method for branch cutting (deep split = 2, cut height = 0.25, minimum cluster size = 30). If the similarity of the modules is > 0.9, they are merged. Based on the level of expression of each gene in each sample, we calculated the correlation between the genes in these modules and the individual phenotypes to measure the correlation between the gene and the phenotype (gene significance). The associations between the modules and variables were assessed to select the relevant modules. The lncRNA-mRNA network visualization was performed via Cystoscope software version 3.7.2 (https://cytoscape.org/) [16].

Module function annotation

The enrichment analysis was conducted by DAVID [version 6.8] (https://david.ncifcrf.gov) [27] GO consists of three parts: biological processes (BP), molecular function (MF), and cellular composition (CC). Besides, all important GO or KEGG terms or genes are filtered into the meaning of P < .05 and at least two mRNAs associated.

Ethics statement

As the data (TCGA and GEO datasets) are publicly available, no ethical approval was required.