LncRNA-mRNA sequencing revealing involvement of disulfide-induced cell death in HNSC
Utilizing lncRNA-mRNA RNA-Seq, we identified 1,151 upregulated lncRNAs and 2,529 upregulated mRNAs in the HNSC cell group compared to the controls. Additionally, 2,277 lncRNAs and 2,603 mRNAs were downregulated, exhibiting a fold change > 2.0 and P < 0.05 in the HNSC cell group (Figure 1A-B). Comprehensive details about the differentially expressed RNAs can be found in (Supplementary Table S2-S3).
Given that the functions of lncRNAs primarily involve the regulation of coding gene expression36, we performed GO enrichment and KEGG pathway analyses to elucidate the molecular functions of the DE lncRNAs. Remarkably, significant GO enrichment terms in all DE mRNAs encompassed cell-cell adhesion, microtubule cytoskeleton, response to oxidative stress, and dynein complex (Figure 1C). Results from the KEGG pathway enrichment analysis of the DE mRNAs revealed clusters of functional pathways, including axon guidance, cell adhesion molecules, proteoglycans in cancers, the Hippo signaling pathway, and the PI3K-AKT pathway (Figure 1D). Within the context of oxidative stress, cells undergo a shift to form protein disulfide bonds, leading to the collapse of the cytoskeleton network—a process crucial for promoting disulfide stress-induced cell death. Notably, we identified representative DE mRNAs associated with the microtubule cytoskeleton, response to oxidative stress, and three key genes linked to disulfidptosis pathways (SLC7A11, MYH10, and TLN1) (Figure 1B, E). These findings collectively suggest the potential involvement of disulfide-induced cell death in HNSC cell lines.
Identification of three DRLs and construction of the risk model
To delve into the role of disulfidptosis in HNSC, we compiled a set of eighteen previously reported disulfidptosis-related genes from existing literature (Supplementary Table S4) and extracted their expression data from HNSC samples (Supplementary Table S5). The visual representation of the study workflow through a flowchart was shown in Figure S1. Employing Pearson correlation analysis, we identified a total of 3,573 lncRNAs. Module identification was conducted using the dynamic tree clipping method, and modules with at least 50 lncRNAs and a similarity greater than 0.75 were merged. The resulting module LncRNAs tree diagram is presented in Figure 2A. Correlation analysis between modules and trait groupings generated by clustering was performed, and a heat map of the module and trait data was plotted (Figure 2B). LncRNAs in the yellow module were selected for further analysis (Supplementary Table S6), demonstrating correlations with one or more disulfidptosis-related genes (Figure 2C).
To ensure the reliability of our findings, we randomly assigned a total of 519 HNSC patients into a training set (n = 260) and a testing set (n = 259). The clinical characteristics of these patient sets are detailed in Table 1, revealing no significant differences in clinical traits between the training and testing sets. Utilizing the univariate Cox regression analysis method, we screened and identified 14 lncRNAs significantly associated with prognosis, labeling them as disulfidptosis-related prognostic lncRNAs (p < 0.05) (Figure 2D). To streamline variables and prevent overfitting in the training cohort, the LASSO Cox regression analysis method was applied, incorporating the lambda value (Figure 2E-F). Ultimately, three lncRNAs were identified according to our lncRNA sequencing results, and a predictive signature risk model was developed based on their expression patterns. The coefficients obtained through regression analysis determined the weighting for each lncRNA in the risk model. The correlation between the three lncRNAs in the signature and disulfidptosis-related genes is illustrated in Figure 2G. Additionally, the risk scores for all samples, calculated based on the expression levels of the three key DRLs, are provided in Supplementary Table S7.
Validation of the three DRLs by qRT-PCR
In the validation cohort, the levels of key DRL genes, namely LINC02434, AC245041.2, and LINC02762, were assessed using qRT-PCR. As illustrated in Figure 2H, all three genes exhibited differential expression in FaDu cells compared to NP69 cells (p < 0.0001). Except for LINC02762, the expression pattern of these three genes consistently aligned with both our RNA-sequencing and database analysis results.
Survival analysis of the risk model and three DRLs
In the survival analysis, we classified the training patients into high-risk and low-risk groups based on their median risk scores. This same categorization was applied to the testing and all patients. In the training set, patients classified as low-risk exhibited significantly better OS compared to those in the high-risk group (p < 0.001) (Figure 3A). This trend persisted consistently observed in both the testing and all patient sets (p = 0.003 and p < 0.001, respectively) (Figure 3B-C). When examining PFS, patients classified as low-risk in the training set demonstrated a higher survival rate compared to those in the high-risk group (p = 0.003) (Figure 3D). This finding was similarly observed in the testing and all patient sets (p = 0.031 and p < 0.001, respectively) (Figure 3E-F). Additionally, for DSS, patients in the low-risk group in the training set had a higher survival rate than those in the high-risk group (p = 0.002) (Figure 3G). This association remained consistent in both the testing and all patient sets (p = 0.013 and p < 0.001, respectively) (Figure 3H-I). These survival analyses highlight the prognostic significance of the risk model, providing valuable insights into their associations with OS, DSS, and PFS in HNSC patients.
Risk curve, survival status, heatmap and independent prognostic value of the risk model
The mortality rate of HNSC patients increased as the risk scores escalated in the training set (Figure 4A-D). This trend was consistently observed in both the testing (Figure 4B, 4E) and all patient sets (Figure 4C, 4F). Expression analysis revealed that LINC02434, AC245041.2 and LINC02762 were more highly expressed in the high-risk group compared to the low-risk group in the training, testing, and all patient sets (Figure 4G-I), suggesting that these DRLs could serve as poor prognostic predictors. By performing both univariate and multivariate Cox regression analyses, we determined that the risk score, calculated based on the risk model, served as an independent prognostic factor for HNSC patients in the training set. The univariate analysis demonstrated a high risk (HR) of 1.784 with a p-value less than 0.001 (Figure 5A), while the multivariate analysis showed an HR of 1.696 with a p-value less than 0.001 (Figure 5D). These findings were further validated in the testing set, with the univariate analysis yielding an HR of 1.549 and a p-value of 0.003 (Figure 5B), and the multivariate analysis resulting in an HR of 1.643 with a p-value of 0.001 (Figure 5E). Similar results were observed in the analysis of the all patients set, where the univariate analysis showed an HR of 1.670 with a p-value less than 0.001 (Figure 5C), and the multivariate analysis revealed an HR of 1.643 with a p-value less than 0.001 (Figure 5F). These results demonstrate that the risk score, derived from the risk model based on the identified DRLs, represents an independent and reliable prognostic factor for HNSC patients. This highlights the potential clinical utility of the risk score in predicting patient outcomes and guiding treatment decisions.
Diagnosis and prognosis of the risk model
To assess the diagnostic value of the risk model, ROC curve analysis was employed. In the training set, the Area Under Curve (AUC) for 1-, 3-, and 5-year survival was 0.636, 0.680, and 0.666, respectively (Figure 6A). While these values did not surpass 0.7, they demonstrated superior performance compared to individual clinical characteristics such as age, gender, grade, and clinical stage in predicting survival at 1, 3, and 5 years, respectively (Figure 6B-D). Similar results were observed in both the testing set (AUC for 1-, 3-, and 5-year was 0.683, 0.683, and 0.619) and the all patients set (AUC for 1-, 3-, and 5-year was 0.660, 0.680, and 0.627) (Figure 6E-L). Thus, the diagnostic value of the risk model surpassed that of individual clinical traits.
Moreover, the C-index, a measure of discrimination, demonstrated higher scores for the risk model compared to other clinical variables in both the training, testing, and all patient sets (Figure S2A-C). This indicates that the risk score derived from the model could serve as a robust prognostic factor, particularly for long-term clinical outcomes in HNSC patients. For practical application, a nomogram was developed, integrating clinical variables and risk scores to predict probabilities of 1-, 3-, and 5-year OS for individual patients (Figure S2D). The nomogram's predictions were validated through calibration plots, which demonstrated good agreement between the predicted and observed clinical outcomes (Figure S2E). Furthermore, survival probabilities and clinical features of HNSC patients across different variables, including age, gender, grade, T stage, N stage, M stage, and clinical stage, were compared. The results consistently showed that high-risk patients had shorter OS compared to low-risk patients across most clinical variables (Figure S3). These findings highlight the applicability of our risk model across different clinical scenarios. In conclusion, our risk model demonstrates both diagnostic and prognostic value in HNSC patients. It outperforms individual clinical traits in predicting survival outcomes, as evidenced by the AUC values, C-index scores, and nomogram predictions.
Samples of HNSC were divided into three tumor clusters
To employ a consensus clustering algorithm based on the expression of the three DRLs, we identified three distinct molecular subtypes within the cohort of 519 HNSC samples (Figure 7A) (Supplementary Table S8). Subsequent survival analysis revealed that patients in cluster 2 exhibited a significantly higher survival rate compared to those in clusters 1 and 3 (p < 0.001) (Figure 7B). This observation can be attributed to the grouping of patients, as most of cluster 2 belonged to the low-risk group, which is associated with a more favorable prognosis (Figure 7C).
To further assess the molecular differences between the identified tumor clusters, we employed PCA and T-distributed stochastic neighbor embedding (tSNE) visualization techniques. Both PCA and tSNE plots demonstrated clear distinctions among the three tumor clusters, indicating robust separability based on the identified molecular subtypes (Figure 7D-E). In addition, we conducted difference analysis among the clinical features of clusters 1, 2 and 3, displaying the proportion of their respective clinical features. We found a statistical difference in N stage (p = 0.021) (Figure 7F). These findings highlight the presence of distinct molecular subtypes within HNSC, contributing to our understanding of the heterogeneity of the disease and offering insights for personalized treatment approaches.
PCA, TMB and TIDE analysis
The PCA method was employed to assess the efficiency of the constructed signature in determining the risk status of patients. The analysis revealed that the risk status could be efficiently determined using the signature derived from the identified DRLs (Figure S5D). However, analyzing only the whole genome expression data (Figure S5A), the disulfidptosis-related genes (Figure S5B) or the DRLs alone (Figure S5C) did not provide an efficient discrimination of the risk status values.
To scrutinize the alterations in somatic mutations between high- and low-risk groups, we retrieved somatic mutation data from the TCGA database. The resulting waterfall plots depict the top 15 mutation genes for both high- and low-risk groups (Figure 8A and Figure 8B). According to the survival analysis, patients in the low-TMB group exhibited a prolonged survival period compared to those in the high-TMB group (p = 0.006) (Figure 8C). Moreover, subgroup survival analysis revealed that patients with both low-TMB and low-risk had the highest survival probability compared to other subgroups (p < 0.001) (Figure 8D). Simultaneously, the low-risk group demonstrated a lower TIDE score than the high-risk group (p < 0.05), suggesting a potentially more effective response to immunotherapy due to a reduced likelihood of immune escape (Figure 8E).
Screening potential drugs for HNSC based on the two risk groups and the three tumor clusters
Employing the "oncoPredict" package, we delved into identifying potentially effective chemotherapeutic agents for HNSC, considering both the two risk groups and the three distinct tumor clusters, which proved meaningful in both contexts. The results, systematically organized in alphabetical order, are presented in (Figure S4A-P). Significantly, Alpelisib, Axitinib, AZD1208, Doramapimod, GSK269962A, JQ1, Sinularin, and Tozasertib demonstrated higher sensitivity (lowest IC50) in low-risk groups (or the highest sensitivity in cluster 2). This implies that patients in low-risk groups (or cluster 2) with HNSC are more likely to experience therapeutic benefits from these specific drugs (All p < 0.001). Further details on the potential drugs can be found in the supplementary materials. These findings illuminate potential treatment options tailored to specific molecular subtypes, paving the way for personalized therapeutic interventions in HNSC patients.