{'fraction_non_alphanumeric': 0.035274735924023644, 'fraction_numerical': 0.01506928965985076, 'mean_word_length': 5.214694369165914, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of an endoplasmic reticulum stress-associated gene signature to predict the immune status and prognosis of cutaneous melanoma

Besides protecting normal cells from various internal and external perturbations, endoplasmic reticulum (ER) stress is also directly related to the pathogenesis of cutaneous melanoma (CM). However, due to the lack of specific molecular biomarkers, ER stress has not been considered a novel treatment target for CM. Here, we identified ER stress-related genes involved in the prognosis of CM patients and constructed an effective model for the prognostic prediction of these patients. First, gene expression data of CM and normal skin tissues from the Genotype-Tissue Expression (GTEx) and The Cancer Genome Atlas (TCGA) databases were retrieved to identify differentially expressed ER stress-related genes in CM. Meanwhile, an independent cohort obtained from the Gene Expression Omnibus (GEO) database was used for validation. The ER stress genes ( ZBP1 , DIABLO , GNLY , FASLG , AURKA , TNFRSF21, and CD40LG ) that were associated with CM prognosis were incorporated into our prognostic model. The functional analyses indicated that the prognostic model was correlated with patient survival, gender, and cancer growth. Multivariate and univariate Cox regressions revealed that the constructed model could serve as an independent prognostic factor for CM patients. The pathway enrichment analysis showed that the risk model was enriched in different immunity and cancer progression-associated pathways. Moreover, the signature model was significantly connected with the immune subtypes, infiltration of immune cells, immune microenvironment, as well as tumor stem cells. The gene function analysis revealed that seven ER stress genes were differentially expressed in CM patients and were significantly associated with prognosis and several antitumor drugs. Overall, our current model presented predictive value for the prognosis of CM patients and can be further used in the development of novel therapeutic strategies for CM. genes, FC = fold change, GSEA = gene set enrichment analysis, LASSO = least absolute shrinkage and selection operator, GO = Gene Ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes, BP = biological process, CC = cellular component, MF = molecular function, ROC = receiver operating characteristic, UPR = unfolded protein response, FOXO = forkhead family transcription factor, TME = tumor microenvironment, ZBP1 = Z-DNA Binding Protein 1, DIABLO = Diablo IAP-Binding Mitochondrial Protein, GNLY = Granulysin, FASLG = Fas Ligand, AURKA = Aurora Kinase A, TNFRSF21 = TNF Receptor Superfamily Member 21, CD40LG = CD40 Ligand, CSCs = cancer stem cell-like cells.


Introduction
Cutaneous melanoma (CM) is an aggressive malignant tumor that threatens human life 1 . The pathogenesis and development of CM are negatively correlated with skin pigmentation. Most cases comprise patients with low skin pigmentation and who were exposed to ultraviolet radiation 2 . The production of melanin involves the synthesis and interaction of multiple proteins in melanosomes. It starts by transforming L-DOPA and L-tyrosine onto melanin polymers that can protect melanocytes from various physical and chemical threats [3][4][5] . Through modulating aerobic glycolysis, oxygen consumption 6,7 , intermediates of melanogenesis 8,9 , and the interaction with several metabolic pathways 10 , melanogenesis can also affect the behavior of malignant and normal melanocytes, and ultimately influence the outcome of melanoma treatments. In 2018, a total of 287,723 people were diagnosed with melanoma worldwide, and 60,709 died due to this disease 11 . The 10-year overall survival (OS) rates of stages I and II CM patients are 75 and 98%, respectively 12 . On the other hand, only 24 and 88% of CM patients in stages IIIA to IIID survived after 10 years, compared to those in stages I and II, suggesting that the early diagnosis of CM might affect its outcome. It has been suggested that skin pigmentation is also involved in tumorigenesis and progression, while the pathogenesis is affected by genetic susceptibility, family history, and acquired melanocytic nevi 13,14 . However, the precise pathogenic mechanisms behind CM remain unknown. Hence, accurate diagnosis in relatively early stages can significantly influence CM therapies. Recently, many investigators have attempted to identify novel biomarkers that can be used for prognostic prediction and personalized therapy of CM patients, but only a few biomarkers of clinical significance were identified 15 . Therefore, the identification of new biomarkers that can accurately predict the prognosis of CM patients is urgently needed.
The endoplasmic reticulum (ER) is a prominent organelle in eukaryotic cells and is involved in the folding and synthesis of transmembrane and secretory proteins, calcium homeostasis regulation, and lipid biosynthesis 16 . During internal and external perturbations, the ER can be induced by several types of stress, such as the disruption of redox homeostasis, nutrient deprivation, and inflammatory stimuli 17 . In skin tissues, the ER stress-induced unfolded protein response (UPR) is necessary for cell differentiation. However, chronic and continuous activation of UPR can ultimately develop into a cell death mechanism 18 and result in some skin diseases, including melanoma 19 . Recently, increasing evidence has shown that ER stress participates in melanoma malignancy and progression by interacting with GRP78/BiP, autophagy, and forkhead family transcription factor (FOXO) pathways [20][21][22] . Moreover, ER stress has become a potential and prevalent target in cancer treatment 23 . Although the connection between ER stress and melanoma is widely accepted, few specific genetic markers of ER stress have been identified in melanoma.
Through bioinformatic analysis, many disease-specific biomarkers have been identified. However, ER stress genes associated with melanoma progression or prognosis were not previously identified by systematic studies. Therefore, we conducted the analysis of differential gene expression and univariate Cox regression in this study, and identified the genes differentially expressed and correlated with the prognosis of CM patients. Then, hub ER stress-related genes were characterized, and a gene risk model was constructed using the least absolute shrinkage and selection operator (LASSO). The prognostic value and clinical significance of this model were further validated in CM patients. Moreover, we analyzed the connections between the ER stress-related gene signature and immune infiltrates, immune microenvironment, relationship with m6A genes, tumor stemness, and cancer chemoresistance. Currently, most models are generated according to tumor immunity and miRNAs, while a thorough analysis of ER stress genes in CM has not been performed yet. Thus, our present study demonstrated that ta risk signature with ER stress-related genes can be used to predict the prognosis of CM patients.

2.1.Datasets
The RNA-sequencing data and related clinical information of CM patients (n = 471) and normal control (n = 1) were retrieved from the TCGA database (https://portal.gdc.cancer.gov) on June 30, 2020. Transcriptome data for 812 normal skin samples were downloaded from the GTEx database. Meanwhile, clinical information and gene profiles of CM patients (n = 214) were obtained from the GEO database (ID: GSE59455 and GSE65904) and used as a cohort for external validation.
Additionally, to remove batch effects, the "sva" R package was used to perform log2transformations and normalize the results 24,25 . Furthermore, the protein domains of ER stress genes (n = 583) were retrieved from the human gene database -GeneCards

2.2.Construction of the prognostic gene signature
A differential expression analysis using the "limma" package was conducted to identify the Differentially expressed ER stress genes (DEERGs). An FDR < 0.05 and | log2 fold change (FC)| > 1 were used to identify the candidate DEERGs. The Kaplan-Meier survival package and univariate Cox regression were employed to identify the prognostic genes associated with ER stress in the TCGA-CM cohort (a p < 0.05 was set as the cutoff). Overlapping genes between the candidate DEERGs and Cox regression analysis were selected and visualized in a Venn diagram. Then, to generate a risk signature for CM and identify hub ER stress-related genes, LASSO analysis was used to integrate the selected genes. The formula risk score = Σexpgenei * βi was used to calculate the risk score for CM patients. In this formula, the relative expression of genes related to pyroptosis is represented as "expgenei" and the regression coefficients as i and β 26 . Next, we divided all CM patients into two groups according to risk scores.
The same regression coefficients and median risk scores were applied to stratify the patients in the validation cohort into two risk groups.

2.3.Evaluation of the risk signature
The "Ggplot2" and "Rtsne" R packages were applied to explore the distribution of the two risk groups in the constructed models by t-SNE and Principal Component Analysis (PCA). The "survival" R package was used to compare the outcomes between the two risk groups, according to their risk levels. The predictive accuracy and clinical characteristics of the risk signature were evaluated using the "timeROC" R package and Cox regression analyses. Moreover, the "ggpubr" R package was used to visualize the connections between clinical characteristics and the risk signature in the TCGA-CM cohort.

2.4.Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses
The Gene Set Enrichment Analysis (GSEA) software version 4.1 was used to compare the KEGG enrichment between the two risk groups. Meanwhile, the Integrated Discovery, Visualization, and Annotation database (version 6.8) 27 was used to establish the biological function of ER stress genes based on GO and KEGG enrichment analyses.
Statistical significance was considered when both FDR and p were < 0.05.

2.5.Tumor Microenvironment (TME) and immune response analyses
The infiltration of stromal and immune cells into tumors was determined using stromal and immune scores 28 . Furthermore, Spearman correlations were employed to examine the relationship between the score of risk and the scores of stromal and immune. Twoway ANOVA was used to identify the connections between risk scores and immune infiltration subtypes. Finally, the two-way ANOVA analysis was carried out to identify the connections of risk scores and immune infiltration subtypes. To determine the gene signature's stem cell-like features, Spearman correlations were conducted to evaluate the connections between tumor stemness and scores of signature.

2.6.Chemotherapy sensitivity analysis
First, 218 chemotherapy drugs (Supplementary were performed to determine the sensitivity of hub ER stress genes to chemotherapy drugs. The "ggplot2", "limma", "ggpubr", and "impute" R packages were employed to visualize these results. .

3.1.Screening of prognosis-associated DEERGs
The datasets analyzed in the present study are presented in Figure 1. A total of 583 ER stress genes were selected from 813 normal and 471 CM tissues, then analyzed to identify DEERGs. Sixty-two ER stress genes were differentially expressed in the TCGA dataset. Among them, eight were connected with the OS of CM patients ( Figure   2A). The distribution of these genes in normal and tumor samples is shown in Figure   2B. Meanwhile, their association with prognosis was verified using univariate Cox analysis ( Figure 2C). These associations were also observed using correlations ( Figure   2D). Finally, eight overlapping ER stress genes were identified. According to the median risk scores, the patients in TCGA (Figure 3A and B) and GEO ( Figure 3C and D) cohorts were divided into high-and low-risk groups. The different directions distributed in the two groups were found in both TCGA ( Figure 3E and F) and GEO (Figure 3G and H) cohorts during the PCA and t-SNE.

3.3.Correlations between the risk scores and clinical characteristics of CM patients
In the TCGA cohort, reduced OS was observed for high-risk CM patients (p < 0.001; Compared with the TCGA cohort, the AUC slightly decreased in the validation cohort, presenting an AUC of 0.619 for one, 0.624 for two, and 0.610 for three years of followup ( Figure 3L). These results confirmed that our signature model could sensitively and specifically predict the OS of CM patients.
Multivariate and univariate Cox regressions showed that, in CM patients, the risk score was an independent prognostic factor, whereas the risk signature was associated with prognosis ( Figure 4A and B). Moreover, higher risk scores were observed in male CM patients (p < 0.05, Figure 4C). Meanwhile, a significant correlation between higher T stages and higher risk scores was observed in CM patients (p < 0.05, Figure 4D).
These analyses indicated that our risk signature is connected to CM development.

3.4.Correlations with immune response, tumor stemness, and m6A-related genes
Compared to the low-risk score subgroup, the levels of nearly all related functions, pathways, and the proportion of immune cell subpopulations (except mast cells) were greatly reduced in the high-risk score group, in both TCGA (Figure 5A and B) and validation ( Figure 5C and D) cohorts. Similar results were obtained using EPIC, MCP counter, XCELL, CIBERSORT, QUANTISEQ, and TIMER ( Figure 5C). Regarding immune checkpoints, all ER stress-related genes identified were more expressed in the low-risk subgroup, except for CD276 ( Figure 5F). Moreover, for tumor promotion and suppression, the immune infiltrates, including wound healing (C1), INF-g dominant (C2), inflammatory (C3), and lymphocyte-depleted (C4), were calculated to explore the relationship between the risk signature and immune components 29 . The results showed that the risk scores were significantly higher for C1 and C4 subtypes and lower for C2 ( Figure 5G). Meanwhile, considering the connection between the ER stress genes identified and immune components, the genes ZBP1, GNLY, FASLG, and CD40LG were upregulated ( Supplementary Figure 2A-D,) and significantly connected with the C2 immune subtype. On the other hand, C1 and C4 subtypes were connected to the downregulation of these genes. Meanwhile, no significant differences were detected between immune subtypes and the risk scores of DIABLO, AURKA, and TNFRSF21 genes (Supplementary Figure 2E-G).
Other possible regulators of CM progression, such as tumor immune microenvironment (stromal and immune scores), m6A-related genes, as well as tumor stemness (DNA methylation pattern and RNA stemness score). The gene signature constructed was significantly and negatively associated with the immune microenvironment (p < 0.05; Figure 5H and I), but positively associated with tumor stemness (p < 0.05; Figure 5J and K). The correlation between hub ER stress genes and tumor immune microenvironment and stemness (Supplementary Figure 3) showed that ZBP1, GNLY, FASLG, TNFRSF21, and CD40LG were positively and FTO were more expressed and WTAP was less expressed in the high-risk group in comparison to the low-risk group ( Figure 5L).
Moreover, considering the roles of the immune checkpoint proteins PD-L1 and PD-L2 in immune progression, we analyzed the correlation between these loci and the gene signature. Both PD-L1 and PD-L2 were significantly less expressed in the high-risk group and were negatively related to the CM risk signature (Figure 6A-D).

3.5.Functional enrichment analyses
In the GO analysis, predominant enrichment of hub genes was observed in different biological processes (BP), including neuron apoptotic process, positive regulation of endothelial cell apoptotic process, and necroptotic process (Figure 7A and B). In the cellular components (CC) category, hub ER stress-related genes were enriched in the CD40 receptor complex, pronucleus, and spindle pole centrosome. Moreover, tumor necrosis factor receptor binding, cytokine activity, and cytokine receptor binding were enriched in the molecular functions (MF) category. Meanwhile, the KEGG enrichment analysis indicated that the genes associated with pyroptosis were enriched in the allograft rejection, autoimmune thyroid disease, cytokine-cytokine receptor interaction, apoptosis, and necroptosis pathways (Figure 7C and D). Additionally, the KEGG enrichment analysis carried out using the GSEA software ( Figure 7E) revealed that, in the high-risk group, 19 pathways, including pyrimidine metabolism, nucleotide excision repair, and RNA polymerase, were significantly enriched (Supplementary Table 3). Meanwhile, in the group with low-risk scores, 46 pathways, including chemokine signaling pathways, Jak stat signaling pathway, and natural killer cellmediated cytotoxicity, were enriched. Altogether, these results indicated possible underlying mechanisms of ER stress genes in CM.

3.6.Prognostic value of the hub genes selected
In the high-risk group, significantly elevated gene expressions were found for the ER stress genes DIABLO and AURKA compared to the low-risk group. Meanwhile, ZBP1, GNLY, FASLG, TNFRSF21, and CD40LG were less expressed (Figure 8A-G). Then, the correlation analyses showed that the expression of DIABLO and AURKA was positively associated, and ZBP1, GNLY, FASLG, TNFRSF21, and CD40LG were negatively associated with risk scores (Figure 8H-N

3.7.Connections between hub ER stress genes and drug sensitivity
Hub ER stress genes were sensitively and respectively correlated to chemotherapy drugs (p < 0.05; Supplementary Table 4). For example, the increased sensitivity to LDK-378 and alectinib was positively correlated with increased ZBP1 expression (Figure 9). In contrast, increased TNFRSF21 expression was negatively associated with the sensitivity to etoposide, teniposide, melphalan, ifosfamide, tfdu, nitrogen mustard, valrubicin, uracil mustard, lomustine, epirubicin, and triethylenemelamine.

Discussion
Although increasing melanoma biomarkers have been identified using next-generation sequencing, novel markers associated with early prognosis and more close detection In the present study, 583 ER stress genes were systematically analyzed to identify Cancer stem cell-like cells (CSCs) promote cancer progression due to their invasion and self-renewal abilities 39,40 . In the present study, the ER stress gene signature was positively connected with stem cell score, confirming that this signature was a risk factor for CM. The m6A-related genes comprehend another tumor research hotspot 41 . Our ER stress gene signature could effectively predict the expression levels of m6A-related genes in CM, including HNRNPC, RBM15, ZC3H13, WTAP, YTHDF1, FTO, and YTHDF2. However, the specific mechanisms underlying these relationships need further exploration.
Despite the prognostic value of the current risk signature, this study also has some limitations. First, the results from our present retrospective study need further confirmation by prospective studies. Second, more experimental assays are needed to verify and validate the conclusions obtained from bioinformatics analyses. In the future, functional studies should be performed to gain mechanistic insights into ER stress genes and their role in CM development.

Conclusions
In the present study, a novel prognostic risk signature consisting of seven hub ER stressrelated genes was constructed and presented high predictive accuracy. This gene signature was valuable to predict parameters related to immune components, immune functions, immune cell infiltration, tumor microenvironment, stemness, m6A-related genes, and drug sensitivity in CM patients. To the best of our knowledge, this is the first ER stress-associated gene signature for CM. These results also provided a novel basis for understanding the specific effects of ER stress genes in CM. Therefore, this study comprehends a significant contribution to the literature and can contribute to improvements in the outcomes and individualized treatments for CM patients

Acknowledgements
Not applicable.

Funding
Not applicable.

Availability of data and materials
The datasets analysed during the current study are available sourced from the publicly available TCGA database (https://portal.gdc.cancer.gov) and GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Ethics approval and consent to participate
Not applicable.

Patient consent for publication
Not applicable.