Head and neck squamous cell carcinoma gene expression and clinical follow-up data from The Cancer Genome Atlas (TCGA) database were downloaded on July 20, 2022. Gene expression data were normalized to log2(FPKM + 1). Furthermore, we screened samples based on the following criteria: 1) The site of Head and Neck Squamous Cell Carcinoma was the larynx, 2) Patients had well-preserved prognostic information, and 3) Patients had a non-zero survival time. As a result, a total of 81 primary laryngeal cancer samples with prognostic information were obtained.

The laryngeal cancer dataset GSE2702011 from the NCBI Gene Expression Omnibus (GEO) database was downloaded, and GSE27020, which included 109 laryngeal cancer samples with prognostic information, was used for validation analysis. Furthermore, 19 cuproptosis-related genes were collected from published literature.6

The expression levels of cuproptosis-related genes were compared between laryngeal cancer patients and healthy participants, whilst Pearson correlations between cuproptosis-related genes were also calculated. Subsequently, based on the differential expression of cuproptosis-related genes between laryngeal cancer patients and healthy participants, we performed tumor subtype analysis of the samples using unsupervised hierarchical clustering R3.6.1 ConsensusClusterPlus (version 1.54.0).12 Finally, the most suitable cuproptosis subtype was obtained when the K-value was set between 2 and 6.

In the R package (version 1.36.2), the Gene Set Variation Analysis (GSVA) algorithm13 was used to calculate enrichment scores between samples from different subtype groups and to then represent the cuproptosis score for each sample. The enrichment scores between different subtypes were investigated using the Wilcoxon test to justify the cuproptosis subtypes.

Subsequently, the Kaplan-Meier curve14 was used to assess the contrast in prognosis between the subtype samples. Furthermore, we compared differences in the clinical information of patients with different subtypes.

CIBERSORT15 was used to calculate the proportions of 22 immune cell compositions in both subtype groups. The Wilcoxon test was used to compare differences in infiltration between the groups and subsequently, the stromal and immune scores of the tumor samples were determined using the ESTIMATE algorithm.16

To observe the possible different molecular mechanisms between various cuproptosis subtypes, we used linear regression and classical Bayesian methods in the limma package (version 3.10.3)17 to perform differential gene expression analysis for each subtype. Subsequently, Benjamini and Hochberg were used to correct the results, and the adj.P.Value was obtained. Next, Differentially Expressed Genes (DEGs) with an adj.P.Value < 0.05 and |log2FC| > 0.5 were defined as inter-subtype specific expressed genes.

Furthermore, we screened genes significantly related to prognosis from specific expressed genes using univariate Cox regression analysis of the survival package in R3.6.1 (version 2.41-1).14 The significance threshold was set at p < 0.05.

A risk model was constructed with stepwise Cox regression analysis of the Survminer package in R3.6.1 (version 0.4.9) based on genes significantly associated with prognosis. The risk score was calculated using the following formula: Risk score = h0(t) × exp ( β1X1 + β2X2 + … +βnXn ).

In this formula, β, h0(t), and h(t,X) refer to the regression coefficients, baseline risk function, and risk function associated with X (covariate), respectively, at time, t.

For the samples’ risk scores in TCGA and GEO datasets, samples were divided into high- (>median risk scores) and low-risk groups (≤median risk score). Prognostic variability between different risk groups was assessed using the Kaplan-Meier curve method.14

Clinical information was compared among the different risk groups. Independent prognostic factors were then screened using univariate and multivariate Cox regression using p < 0.05 as a threshold, and forest plots were drawn. We then collated independent prognostic factors and constructed a nomogram using the rms package (version 5.1-2).18

Based on the Wilcoxon test, the expression of immune checkpoint genes and HLA family genes was compared between patients in different risk groups.

Enrichment of significant hallmark gene sets (h.all.v7.4. symbols) between different risk groups were investigated using GSEA. The threshold for screening significantly different genes set was set at p < 0.05 and ¿NES¿ > 1.

Using the ggalluvial package (version 0.12.3) in R3.6.1, we compared the proportions of the two cuproptosis subtypes between the different risk groups and investigated the association between these subtypes and subgroups.

The expression patterns of 19 cuproptosis-related genes were analyzed across laryngeal cancer patients and healthy participants, with 12 of these genes showing significantly different expression patterns (Fig. 1A) (p < 0.05). Furthermore, we investigated the correlations between these 12 cuproptosis-related genes and found negative correlations between the majority of them (Fig. 1B).

Figure 1.

Expression patterns of cuproptosis-related genes. (A) Expression of cuproptosis-related genes in laryngeal cancer and normal samples. (B) Correlation analysis of differentially expressed cuproptosis-related genes.

(0.62MB).

Based on the 12 cuproptosis-related genes differentially expressed in healthy and cancerous samples, unsupervised clustering analysis of tumor samples was performed. Subsequently, two of these subtypes were chosen as the best K-values (Fig. 2 A–C). In addition, to further validate the plausibility of the subtypes, we compared cuproptosis scores between subtypes, with these results restricting the cuproptosis scores of patients with the C1 subtype (27 samples) to be significantly higher than the C2 subtype (54 samples) (Fig. 2D) (p < 0.05).

Figure 2.

Identification of two cuproptosis subtypes in laryngeal cancer samples. (A) Cumulative Distribution Function (CDF) for consensus clustering. (B) Consensus clustering matrix when K is 2. (C) Proportion of ambiguous clustering analysis. (D) Cuproptosis scores of C1 and C2 cuproptosis subtypes.

(0.38MB).

Kaplan-Meier curves showed that patients with laryngeal cancer in the C2 subtype had a significantly better prognosis than those in the C1 subtype (Fig. 3A) (p < 0.05). We also compared the expression patterns of 12 copper death genes in patients of both subtypes (Fig. 3B). Subsequently, the proportions of patients with the C1 and C2 subtypes were also compared under different clinical characteristics. The proportion of C2 subtypes was higher for patients with a history of smoking (Fig. 3C and Supplementary Fig. 1).

Figure 3.

Relationship between cuproptosis subtypes and survival of laryngeal cancer patients. (A) Comparison of survival differences between patients with C1 and C2 subtypes using Kaplan-Meier curve. (B) Heat map of the distribution of expression levels of 12 differential cuproptosis-related genes in subtype samples. (C) Comparison of the proportion of patients with C1 and C2 subtypes of laryngeal cancer among patients with different smoking histories.

(1.17MB).

The fraction of immune cells in patients with subtypes C1 and C2 was compared. Between the two subtypes, 5 of the 22 immune cells were significantly different (Fig. 4A). Specifically, the proportion of plasma cells, memory B cells, and CD8 T cells was higher in the C2 subtype, whereas the proportion of M0 macrophages and CD4 memory resting T-cells was higher in the C1 subtype (p < 0.05). Further investigation also indicated that stromal, immune, and ESTIMATE scores were always higher in the C1 subgroup than in the C2 subgroup (Fig. 4B).

Figure 4.

Immune analysis. (A) Fraction of 22 immune cells in C1 and C2 samples. (B) Comparison of the stromal, immune, and ESTIMATE scores between patients in the C1 and C2 subtypes.

(0.53MB).

By comparing the gene expression between subtypes C1 and C2, we obtained 218 DEGs (211 upregulated and 7 downregulated genes) (Fig. 5A, Supplementary Table 1). We screened the prognostic genes from 218 DEGs using a univariate Cox analysis (Fig. 5B). Finally, 14 genes were defined as significantly related to patient prognosis, of which 13 genes: Transmembrane 2 (TMEM2), Stanniocalcin 1 (STC1), Dishevelled binding Antagonist of β-Catenin 1 (DACT1), Parvalbumin (PVALB), Myosin Light chain 4 (MYL4), Stathmin 2 (STMN2), Secretogranin II (SCG2), Potassium voltage-gated Channel subfamily D member 2 (KCND2), Gastrin Releasing Peptide Receptor (GRPR), Myelin-associated Oligodendrocyte Basic Protein (MOBP), Interleukin 1 Receptor I (IL1R1), G Protein-coupled Receptor 173 (GPR173), and Matrilin-3 (MATN3); were risk factors, whilst 1 gene (MAGIX) was a protective factor.

Figure 5.

Identification of the prognostic genes of patients with laryngeal cancer. (A) Volcano map was used to show differentially expressed genes between C1 and C2 subtypes. (B) Forest plot was used to show prognosis-related genes screened by univariate Cox regression analysis.

(0.39MB).

We used the stepwise Cox regression algorithm to screen for the best combination of the 14 genes significantly associated with prognosis for risk model construction. TMEM2, DACT1, STMN2, and GPR173 were used to construct the risk model (Supplementary Fig. 2). After classifying the samples into different groups, we found that there were more deaths among the patients in the high-risk group (Fig. 6A) and moreover, the high-risk group had a worse prognosis (Fig. 6B). The Area Under the Curve (AUC) of the model predicting patients at 1, 3, and 5 years were 0.901, 0.798, and 0.771, respectively (Fig. 6C). Based on the median risk score, patients in the GEO dataset were grouped to validate the analysis results in the TCGA dataset. GEO results also indicated that the high-risk group had more deaths and a worse overall prognosis (Fig. 6 D,E). The Receiver Operating Characteristic (ROC) curve in the GEO dataset also proved to be a good prediction performance model (Fig. 6F).

Figure 6.

Predictive performance of risk model for prognosis of laryngeal cancer patients. (A) Survival status of patients with laryngeal cancer in high- and low-risk groups in TCGA database. (B) Kaplan-Meier curve of laryngeal cancer patients in high- and low-risk groups in TCGA database. (C) ROC curves were used to present the predictive performance of the risk model for 1-, 3-, and 5-year survival of patients with laryngeal cancer in TCGA database. (D) Survival status of patients with laryngeal cancer in high- and low-risk groups in GEO database. (E) Kaplan-Meier curve of laryngeal cancer patients in high- and low-risk groups in GEO database. (F) ROC curves were used to present the predictive performance of the risk model for 1-, 3-, and 5-year survival of patients with laryngeal cancer in GEO database.

(0.52MB).

Primary Tumor (T), clinical regional lymph Nodes (cN), distant Metastasis (M), age, sex, clinical stage, grade, and risk groups were all included in the univariate Cox regression analysis. Significant results were further included in this analysis (Fig. 7A). Finally, distant metastasis M, clinical regional lymph nodes cN, and risk groups were identified as independent prognostic factors for patients with laryngeal cancer and were used to construct a subsequent nomogram (Fig. 7B). The C-index of this nomogram was 0.809 (Fig. 7C), and the calibration curve showed a good agreement between the predicted and observed values of the column line graph model for the prediction of overall survival at 1, 3, and 5 years in patients with laryngeal cancer (Fig. 7D).

Figure 7.

Predictive performance of nomogram model for prognosis of laryngeal cancer patients. (A and B) Univariate and multivariate Cox regression analysis to screen prognostic independent factors in patients with laryngeal cancer. (C) A nomogram constructed by the prognostic independent factors. (D) Calibration curves were used to demonstrate the difference between the 1-, 3-, and 5-year survival predicted by nomogram for patients with laryngeal cancer and the actual survival. (E) Kaplan-Meier curve of laryngeal cancer patients in high- and low-nomogram groups. (F) ROC curves were used to present the predictive performance of the nomogram model for 1-, 3-, and 5-year survival of patients with laryngeal cancer.

(0.42MB).

Patients were then grouped based on this nomogram. The Kaplan-Meier curve indicated that the prognoses of laryngeal cancer patients in the high nomogram group were significantly poorer than those in the low nomogram group (Fig. 7E). In addition, the ROC curve further demonstrated that the prognosis of laryngeal cancer patients at 1, 3, and 5 years was well predicted using the nomogram, with respective AUCs of 0.814, 0.859, and 0.782 (Fig. 7F).

The expression of Signal regulatory protein alpha (MYD1) and Receptor induced by lymphocyte activation (4-1BB) was significantly greater in the high-risk group than in the low-risk group (Supplementary Fig. 3A) (p < 0.05), although further studies showed no significant differences in HLA gene family expression between groups (Supplementary Fig. 3A) (p < 0.05).

Hallmark gene sets were screened with p < 0.05 and ¿NES¿ > 1. Finally, we obtained 11 hallmark gene sets that were significantly enriched in the high-risk group, such as angiogenesis (NES = 1.8842, NP = 0.0000), epithelial-mesenchymal transition (NES = 1.9842, NP = 0.0000), TGF beta signaling (NES = 1.9992, NP = 0.0020), and Kirsten Rat Sarcoma viral oncogene homolog (KRAS) signaling (NES = 1.7271, NP = 0.0020) (Fig. 8A). Oxidative phosphorylation (NES = −1.8737, NP = 0.0159) was the only hallmark gene set that was significantly enriched in the low-risk group (Fig. 8A).

Figure 8.

(A) Gene set enrichment analysis. (B) Proportion of cuproptosis subtypes C1 and C2 in patients in high- and low-risk groups. (CF) Evaluation of treatment response in patients with laryngeal cancer. (C, D, and E) IC50 of BX.759, pazopanib, and PLX4720 for patients with laryngeal cancer in high- and low-risk groups. (F) TIDE score of patients with laryngeal cancer in high- and low-risk groups.

(0.39MB).

For the different risk groups, we used the chi-square test to compare the proportions of subtypes C1 and C2. The results indicated that the proportion of patients with the C2 subtype was significantly higher for the low-risk group than in the high-risk group (Fig. 8B, p < 0.05).

We investigated the sensitivity of patients in the different risk groups to 138 chemotherapeutic agents. The IC50 results revealed a significant difference in the sensitivity of patients in the different groups to 37 chemotherapeutic agents (Supplementary Table 2). The top three most significant drugs were N-[3-[[5-Iodo-4-[[3-[(2-Thienylcarbonyl)amino]propyl]amino]-2-pyrimidinyl]amino]phenyl]-1-pyrrolidinecarboxamide (BX.795), pazopanib, and N-[3-[(5-Chloro-1H-pyrrolo[2,3-b]pyridin-3-yl)carbonyl]-2,4-difluorophenyl]-1-propanesulfonamide (PLX4720) (Fig. 8D–F). The Tumor Immune Dysfunction and Exclusion (TIDE) study revealed that patients in the high-risk group had significantly higher TIDE scores than those in the other groups (Fig. 8G) (p < 0.05).

We assessed patient sensitivity to chemotherapeutic agents using the Greedy Double Subspaces Coordinate Descent (GDSCD) (https://www.cancerrxgene.org/). To quantify IC50 values, the pRRophetic package41 in R was used, whilst the Wilcoxon test was used to investigate the differences in IC50 values among 138 chemotherapeutic agents across risk groups.

The TIDE database was used to assess patient responses to immune checkpoint therapy. For patients in different risk groups, we used the Wilcoxon test to compare their respective TIDE scores.