Effective prognostic risk model with cuproptosis-related genes in laryngeal cancer

Highlights • Cuproptosis genes stratified laryngeal cancer into diverse prognostic subtypes.• TMEM2, DACT1, STMN2, and GPR173 form our prognostic model.• TMEM2, DACT1, STMN2, GPR173 shape the risk model.• High-risk samples upregulated angiogenesis, EMT, and KRAS pathways.• Clinical_M, Clinical_cN, risk groups predict laryngeal cancer prognosis.


Introduction
Laryngeal cancer is the second most common malignant tumor of the head and neck. 1 The cure rate for early stage laryngeal cancer is higher than 80%, although the recurrence rate for advanced laryngeal cancer is 25%---50%. 2,3Furthermore, the 5-year survival rate for stage IV laryngeal cancer patients is even lower at 40%, 2 thus, an early diagnosis would increase the likelihood of curing the disease.Patients with an early diagnosis of laryngeal cancer can often be treated with surgery or radiation, which ultimately preserves the larynx.In contrast, patients with advanced disease often require multimodal treatment and are less often laryngeal preserving. 4Currently, there are no biomarkers for the early and prognostic detection of laryngeal cancer.Therefore, here we investigate whether clinical markers widely used in a variety of diseases could also play a role in laryngeal cancer.
Cell death caused by abnormal iron levels is called ferroptosis, 5 and recent studies have now also identified intracellular copper-induced cell death called cuproptosis. 6opper can bind to the lipid-acylated components of the tricarboxylic acid cycle, causing proteotoxic stress and eventually inducing ferroptosis in cells distinct from other apoptotic pathways. 7Genes related to Cu 2+ transport, such as recombinant copper transporter 1 (SLC31A1), affect dietary copper uptake and have a possible role in regulating cuproptosis. 8Additional studies have demonstrated the potential impact of cuproptosis-related genes on the immune microenvironment, clinicopathological features, and prognosis of triple-negative breast cancer. 9For example, Huang et al. developed a scoring system to predict prognosis and immunity in patients with squamous cell carcinoma of the head and neck based on cuproptosis-related genes, which showed good predictive results. 10However, no risk models have yet been constructed using copper deathrelated genes that can predict the prognoses of patients with laryngeal cancer.
In the present study, we investigated expression patterns of cuproptosis-related genes in patients with laryngeal cancer and normal participants.We expect the risk model constructed using cuproptosis-related genes to have good prognosis predictive performance for laryngeal cancer patients and provide new therapeutic methods for laryngeal cancer therapy.

Data collection and preprocessing
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 log 2 (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 GSE27020 11 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

Identification of cuproptosis subtypes in laryngeal cancer
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). 12inally, the most suitable cuproptosis subtype was obtained when the K-value was set between 2 and 6.

Prognostic analysis of of cuproptosis subtypes
In the R package (version 1.36.2), the Gene Set Variation Analysis (GSVA) algorithm 13 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 curve 14 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.

Comparison of immune microenvironment between cuproptosis subtypes
CIBERSORT 15 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

Identification specific and prognostic genes between cuproptosis subtypes
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) 17to 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). 14The significance threshold was set at p < 0.05.

Construction of prognostic model
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 = h 0 (t) × exp ( ␤ 1 X 1 + ␤ 2 X 2 + . . .+␤ n X n ).
In this formula, ␤, h 0 (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

Construction of nomogram model
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

Comparison of immune checkpoint genes and human leukocyte antigen (HLA) family genes between patients in different risk groups
Based on the Wilcoxon test, the expression of immune checkpoint genes and HLA family genes was compared between patients in different risk groups.

Gene set enrichment analysis (GSEA)
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.

Association of risk group with cuproptosis subtypes
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.

Expression pattern of cuproptosis-related genes in laryngeal cancer patients
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).

Two cuproptosis subtypes were identified in laryngeal cancer samples
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).

Correlation of cuproptosis subtypes with clinical features
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).

Immune analysis
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).

Screening for genes significantly associated with prognosis
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

Constructing a risk model for patients with laryngeal cancer
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).

Nomogram was constructed with independent prognostic factors
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).
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).

Differences in immune checkpoint genes and HLA family genes between patients in different groups
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).

Cuproptosis subtypes were associated with risk groups
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).

Differences in treatment sensitivity
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- 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).

Discussion
Laryngeal cancer accounts for 25% of squamous cell carcinomas of the neck and head, yet there is still a lack of effective prognostic markers. 19In recent years, cuproptosis, a specific cell death mechanism, has been associated with the development of several solid tumors. 20Here, we used cuproptosis-related genes to construct a cuproptosis subtype that differs from previous clinical staging and could classify patients into different subgroups.There were sig- nificant differences found in the prognoses of patients with laryngeal carcinoma.In addition, a risk model constructed using cuproptosis-related genes that were significantly associated with laryngeal cancer patients performed well in the prediction of patient prognoses.This risk model was also a prognostic independent factor for patients with laryngeal cancer.Independent prognostic factors, including the risk model, Clinical M, and Clinical N, were used to construct a nomogram.
The risk model we constructed using cuproptosisrelated genes included four genes.TMEM2 is a protein with a single-channel transmembrane domain protein. 21ther investigators have previously demonstrated that TMEM2 promoted developmental angiogenic Vegf signaling by regulating hyaluronic acid turnover. 22TMEM2 also promoted breast cancer cell metastasis, 23 although no correlation between TMEM2 expression and laryngeal cancer has been reported as yet.In the present study, TMEM2 expression was a risk factor for laryngeal cancer.High expression of the DACT1 may have been invasive and metastatic in nasopharyngeal carcinoma cells. 24DACT1 is methylated in all oral squamous cell carcinomas. 25TMN2 encodes a stathmin family member of phospho-proteins and is upregulated in squamous cell carcinoma and lung adenocarcinoma. 26GPR173 is part of a highly conserved receptor family. 27In addition, no correlation between GPR173 and laryngeal cancer has been previously reported.
MYD1 is known as Signal Regulatory Protein Alpha (SIRPA) and was upregulated in the high-risk group.Furthermore, MYD1 was found to be a receptor for CD47 and was significantly upregulated in hot tumors. 28The binding of SIRPA to CD47 triggers inhibitory pathways in some cancers, ultimately leading to the immune escape of tumor cells. 29n addition, the phagocytic activity of macrophages and other phagocytes in tumors is enhanced when CD47-SIRPA interactions are blocked. 30The expression of 4-1BB in activated T-cells can cause a signaling cascade that leads to the upregulation of anti-apoptotic molecules.Targeting 4-1BB with agonist antibodies in tumor models can induce clear and durable antitumor immunity. 31Agonistic anti-human 4-1BB antibodies and second-/third-generation 4-1BB/CD3 Chimeric Antigen Receptor (CAR) T-cells have entered clinical trials. 32The two groups responded differently to numerous chemotherapy drugs.BX-795 is an inhibitor of 3-phosphoinositide-dependent kinase 1, which plays an active role in the treatment of many solid tumors. 33ifferent sensitivities to chemotherapeutic agents may also indicate different outcomes for patients in different risk groups receiving treatment.
GSEA revealed that high-risk samples showed upregulated expression of pathways involved in angiogenesis, Epithelial-Mesenchymal Transition (EMT), and KRAS signaling.Angiogenesis is the process by which new capillaries grow out of existing vessels, and abnormal angiogenesis is an important sign of cancer development. 34Studies have found increased vascularity and invasion in cases of head and neck tumors. 35Schlüter et al. found that high expression of angiogenic markers was associated with a high risk of squamous cell carcinoma of the larynx. 36timulation with bile acids and pepsin can promote carcinogenesis in the pharynx through EMT and is a risk factor for laryngeal cancer. 37,38Aberrant KRAS expression leads to uncontrolled cell proliferation and ultimately moves the cells in a direction that favors metastasis. 39It has also been found that KRAS overexpression was associated with the progressive dedifferentiation of laryngeal squamous cell carcinoma. 40r analysis primarily relies on retrospective data from specific databases, which may not represent the global diversity of laryngeal cancer cases.Additionally, the molecular mechanisms behind the observed associations with cuproptosis-related genes require further elucidation.While our findings offer promising avenues for future research and potential therapeutic targets, they need to be validated through prospective studies and integrated with other known biomarkers of laryngeal cancer.

Conclusion
In conclusion, cuproptosis-related genes were used to subtype laryngeal cancer.The constructed cuproptosis risk model showed a good predictive performance for the prognosis of patients with laryngeal cancer.The risk model and nomogram produced here may be useful for guiding future clinical practice.This study also provides potential targets for laryngeal cancer therapy and a basis for studying the mechanisms by which cuproptosis plays a role in laryngeal cancer.

Therapeutic analysis
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 package 41 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.

Figure 1
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.

Figure 2
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.

Figure 3
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.

Figure 4
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.

Figure 5
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.

Figure 6
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.

Figure 7
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.

Figure 8 (
Figure 8 (A) Gene set enrichment analysis.(B) Proportion of cuproptosis subtypes C1 and C2 in patients in high-and low-risk groups.(C F) 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 highand low-risk groups.