A Novel Pyroptosis-Related lncRNA Signature for Predicting the Prognosis of Skin Cutaneous Melanoma

Abstract

Purpose

To construct a prognosis model of melanoma based on pyroptosis related genes.

Methods

Melanoma RNA-sequencing data was downloaded from TCGA. First, the lnRNAs related to pyroptosis were obtained through Pearson correlation analysis. Then, the prognosis model of pyroptosis related genes was constructed by Cox regression and Lasso regression. Melanoma patients were divided into high-risk and low-risk groups by risk score, and the differences in prognosis and immune microenvironment between the two groups were explored.

Results

We found that the high-risk group had a significantly poorer prognosis, and different groups differed in immune infiltration, m6A methylation, and immune checkpoint.

Conclusion

Our prognostic model can provide a reference for the study of pyroptosis in melanoma cells and provide a new idea for melanoma treatment.

Keywords:

• melanoma
• pyroptosis
• programmed cell death
• immunotherapy
• bioinformatics

Introduction

Melanoma is a highly malignant skin tumor with strong invasiveness and poor prognosis.¹ The National Cancer Institute estimates that there will be about 100,000 new cases of melanoma and 6850 deaths in the United States in 2020.² Risk factors for melanoma include light skin color, UV exposure, nevi, genetic predisposition, and a family history of melanoma.² Since the beginning of the 21st century, research on the treatment of melanoma, especially the discovery of immune checkpoint inhibitors and BRAF inhibitors, has greatly improved the survival rate of melanoma patients.^3,4 However, although the mortality of melanoma patients is on the decline, the morbidity of melanoma is still increasing year by year.⁵ At the same time, drug resistance of PD-1/PD-L1 inhibitors and BRAF mutation inhibitors has also gradually emerged, which makes the exploration of new biomarkers urgent.⁶

Pyroptosis, an inflammatory form of programmed cell death, shares many similarities with apoptosis, but there are also significant differences.⁷ When cell pyroptosis occurs, there will be many pores on the cell membrane so that the cell will have osmotic swelling and membrane rupture, and the inflammatory substances will be discharged, causing an obvious inflammatory reaction.⁷ This is the morphological characteristic of cell pyroptosis. As early as 1992, Zychlinsky et al observed that Shigella flexneri can induce this death in host macrophages, but due to limited knowledge at the time, it was still considered to be apoptosis.⁸ It was not until 2001 that Cookson et al found that this death was different from apoptosis and had Caspase-1-dependent characteristics, which induced them to name it “pyroptosis” for the first time.⁹ “Pyro” means “fire” and “ptosis” means “fall”, which vividly describes the characteristics of this type of death. According to the differences in the mechanism of cell pyroptosis, it can be divided into the classical pyroptosis pathway mediated by Caspase-1 and the non-classical pyroptosis pathway mediated by Caspase-4/5/11.¹⁰ Activated caspase cleaves Gasdermin protein to form the N-terminal domain (GSDM-N) and C-terminal domain (GSDM-C).¹¹ GSDM-N has pore-forming activity and is “oligomerized” to form pores of varying sizes in the cell membrane, mediating osmotic swelling of the cell and the expulsion of inflammatory molecules such as IL-1β and IL-18.⁷ Pyroptosis is a double-edged sword in cancer. On the one hand, we can induce Pyroptosis of tumor cells to inhibit tumor growth. On the other hand, the inflammatory substances released in the process of cell pyroptosis form the inflammatory microenvironment of the tumor and promote tumor growth by activating different signaling pathways.¹² The correlation between pyroptosis and cancer still needs further research to clarify. Among them, the correlation between pyroptosis and melanoma is still unclear.

Now it’s been 60 years since Jacob and Monod discovered mRNA.¹³ Their momentous discovery laid the foundation for the “genetic central dogma”.¹⁴ To this day, the “genetic central dogma” remains the classic law of genetics. But in 2001, the Human Genome Project found that only 1.5% of nucleic acid sequences are transcribed into mRNA to encode proteins, while the remaining 98.5% are often considered “junk sequences”.¹⁵ With the development of research, people found that these “junk sequences” can be transcribed into a variety of non-coding RNAs. LncRNA is one of them, which usually refers to non-coding RNA transcripts with more than 200 nucleotide sequences, and plays an important role in homeostasis regulation and the occurrence and development of tumors.¹⁶ At present, many studies have confirmed the relationship between lncRNA and cell apoptosis, but the specific role network is still unclear.

In this study, we explored pyroptosis-related lncRNAs in melanoma and constructed a prognostic model, which showed that patients with a high-risk score had a worse prognosis. Our conclusions can provide some new ideas for the prognosis and treatment of melanoma.

Methods

Information Extraction of Datasets

Melanoma RNA sequencing data downloaded from TCGA (https://portal.gdc.cancer.gov/) in the database. Inclusion criteria were: (1) patients diagnosed with melanoma; (2) Patients had detailed lncRNA and clinical information. A total of 457 melanoma patients met the inclusion criteria, and their lncRNA and clinical information were downloaded (those with a follow-up of fewer than 30 days were excluded).

Identification of lncRNAs Associated with Cell Pyroptosis

Thirty-three apoptosis-related genes were obtained from previous literature. The Pearson correlation method was used to identify lncRNAs related to cell pyroptosis. The correlation coefficients R²>0.3 and P<0.001 were considered to be pyroptosis-related lncRNAs. We used Cytoscape software to visualize lncRNA co-expression networks.

Identification of Prognostic Pyroptosis-Related lncRNAs

First, we performed univariate Cox regression on pyroptosis-related lncRNAs. We used the package “forestplot” in R software for COX regression analysis. Then Lasso regression was performed for lncRNAs with P<0.05, and multivariate Cox regression was performed for the results obtained. According to the final results, risk scores were calculated. Using the regression coefficient (β), the scoring formula is: . Using a median risk score, patients were divided into high-risk and low-risk groups.

Construction of the Prognostic Model

After patients were divided into high-risk and low-risk groups, Kaplan-Meier survival analysis was performed for both groups to identify the prognostic value of the model. The accuracy of the established prognosis model was verified by the calibration curve and ROC curve. We used the “plotROC” package in R software to draw ROC curves.

Clinical Value of the Established Prognostic Model

We used univariate and multivariate Cox regressions, respectively, to test whether risk score and clinical characteristics (age, sex, stage) were effective prognostic indicators for melanoma patients. Nomograms were used to demonstrate clinical characteristics and risk scores concerning 1, 3, and 5-years survival.

Correlation Analysis of Immune Checkpoint and M6A-Methylation

We then collected 47 immune checkpoint-related genes (ICPs) and 25 M6A-methylated related genes by reading the published literature and investigated the expression of these genes in the high- and low-risk groups in two data sets. Differentially expressed ICP and ICD genes were retained (P <0.05) and displayed as a bar chart. Rank sum test was used for comparison between groups.

Model to Predict Potential Drug Candidates for Melanoma

By using the “pRRophetic” package and the expression matrix of melanoma patients, we predicted the minimum drug inhibition concentration (IC50) of drugs in melanoma patients, and finally obtained some drugs that are related to the model and may become candidates for the treatment of melanoma.

Gene Set Enrichment Analysis

Gene Set Enrichment Analysis (GSEA) is a computational method that determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states (eg phenotypes). We used GSEA analysis to explore KEGG pathways differences in gene function between the high-risk and low-risk groups.

Results

The Flow Chart

The flow chart is shown in Figure 1. First, we downloaded melanoma data from TCGA database. Then, we constructed a prognostic model based on 9 pyroptosis related lncRNAs. Subsequently, we investigated differences in immune infiltration and m6A methylation between the high-risk and low-risk groups. Finally, we performed GSEA analysis.

Figure 1 The flow chart of our study.

Nine Pyroptosis-Related lncRNAs Were Identified to Calculate the Risk Score

Firstly, the lncRNA closely associated with pyroptosis-related genes in melanoma was identified by Pearson-correlation test. Cox regression and Lasso regression were performed on the pyroptosis-related lncRNAs of melanoma in the TCGA database, and 9 pyroptosis-related lncRNAs with prognostic value were finally obtained (Table 1). The risk score was calculated by these 9 lncRNAs. Based on the median risk score, 457 melanoma patients in TCGA were divided into high-risk and low-risk groups (Figure 2A).

Table 1 9 Pyroptosis-Related lncRNAs with Prognostic Value

LncRNA	Coef	p-value
AL121603.2	−0.783701523	<0.05
AC107464.2	−0.303872015	<0.05
AC245128.3	−0.67438924	<0.05
AC092171.5	0.798948263	<0.05
AC242842.1	−0.147658944	<0.05
IRF2-DT	−0.174729664	<0.05
HLA-DQB1-AS1	−0.043901003	<0.05
AC004585.1	0.190691242	<0.05
LINC00582	−0.257021589	<0.05

Figure 2 (A) The risk score was calculated by these 9 lncRNAs. Based on the median risk score, 457 melanoma patients in TCGA were divided into high-risk and low-risk groups. (B) Prognostic differences between the high-risk and low-risk groups were investigated. High risk group was found to indicate poor prognosis (P <0.001). (C) In the form of dot plot, we can intuitively see that the survival time of patients in the high-risk group is concentrated at the bottom, that is, the survival time is shorter. (D) Heat map reflects the relationship between the expression status of 9 pyroptosis-related lncRNAs and clinical features such as TNM stage, sex and age.

The Group Model Based on Risk Score Can Effectively Predict the Prognosis of Melanoma Patients

After classifying melanoma patients into high-risk and low-risk groups, we performed a survival analysis. As can be seen from the Figure 2B, the survival rate of patients in the high-risk group decreased significantly, and the difference was statistically significant (P <0.001). In Figure 2C, in the form of dot plot, we can intuitively see that the survival time of patients in the high-risk group is concentrated at the bottom, that is, the survival time is shorter. Figure 2D heat map reflects the relationship between the expression status of 9 pyroptosis-related lncRNAs and clinical features such as TNM stage, sex and age. Heatmap showed significant differences in T staging between patients in the high- and low-risk groups. There was no significant difference in other clinical features.

Risk Score is an Independent Factor Affecting Patient Survival with High Accuracy

In order to explore whether risk score is an independent influencing factor for melanoma patients, univariate and multivariate Cox regression were performed for risk score, age, gender, TNM stage and other factors. The clinical characteristics of the patients in the database are as follows (Table 2). Univariate Cox regression showed that age, stage, T stage, N stage and risk score were independent influencing factors for the prognosis of melanoma patients (Figure 3A). Multivariate Cox regression showed that T stage, N stage and risk score were independent influencing factors for the prognosis of melanoma patients (Figure 3B). The ROC curve showed that the AUC value of the risk score was 0.821, demonstrating its sensitivity (Figure 3C). Nomogram can see that the risk score is associated with 1, 3, and 5 year survival (Figure 3D). Calibration shows that Nomogram can accurately predict 1, 3, and 5-year survival of patients with melanoma (Figure 3E).

Table 2 The Clinical Characteristics of the Patients in the Database

Clinical	Group	n	Mean	SD	t	P
Age	>65	120	1.445	0.897	1.895582261	0.059
	≤65	222	1.265	0.71
Gender	Female	131	1.316	0.781	−0.221189375	0.825
	Male	211	1.336	0.788
Stage	1–2	188	1.299	0.746	−0.762915393	0.446
	3–4	154	1.364	0.83
T	0–2	128	1.079	0.659	−4.928465071	0
	3–4	214	1.477	0.816
M	0	331	1.329	0.792	0.205587713	0.841
	1	11	1.295	0.526
N	0	196	1.294	0.737	−0.909634201	0.364
	1–3	146	1.374	0.844

Figure 3 (A) Univariate Cox regression showed that age, stage, T stage, N stage and risk score were independent influencing factors for the prognosis of melanoma patients. (B) Multivariate Cox regression showed that T stage, N stage and risk score were independent influencing factors for the prognosis of melanoma patients. (C) Among the ROC curves of patient prognosis and clinical characteristics, the risk score AUC was the highest (0.821). (D) Nomogram was used to predict 1, 3, 5-year survival. (E) Calibration shows that Nomogram can accurately predict 1, 3, and 5-year survival of patients with melanoma.

Correlation Analysis of Immune Checkpoint and M6A Methylation in Different Groups

Many immune checkpoint genes were expressed differently in samples in the two groups and tended to be higher in the low-risk group (***P<0.001, Figure 4A). Among the samples from the high-risk and low-risk groups, the expression of 7 M6A-methylation related genes differed, and all of them showed a tendency to be upregulated in the low-risk group (***P<0.001; * P < 0.05; Figure 4B).

Figure 4 (A) Immune checkpoint genes that were differentially expressed in the high-risk and low-risk groups were selected and presented in a boxplot. (***p < 0.001) (B) Methylation-related genes that were differentially expressed in the high-risk and low-risk groups were selected and presented in a boxplot. (*p < 0.05, **p < 0.01, ***p < 0.001) (C) Heat maps of immune cells that differ between the high-risk and low-risk groups.

Analysis of Immune Cell Infiltration Levels in Different Groups

The levels of immune infiltration in high-risk and low-risk groups were calculated by algorithm such as TIMER, CIBERSORT, Quantiseq and other methods. The heat map shows the level of each immune cell infiltration in the tumor microenvironment of the two groups of patients under different algorithms (Figure 4C). We can see that the low-risk group tended to have higher levels of immune cells in the tumor microenvironment.

Analysis of Sensitivity Difference of Antitumor Drugs in Different Groups

The study of the sensitivity of different groups of patients to anti-tumor drugs can provide help for the formulation of future treatment regimens. The high-risk group was more sensitive to certain drugs, such as Bexarotene, Bryostatin and Docetaxel (Figure 5A). The low-risk group was more sensitive to Bortezomib, Bosutinib, Camptothecin, etc (Figure 5B). This has guiding significance for us to select drugs according to the group situation.

Figure 5 (A) The high-risk group was more sensitive to certain drugs, such as Bexarotene, Bryostatin and Docetaxel. (B) The low-risk group was more sensitive to Bortezomib, Bosutinib, Camptothecin, etc.

Gene Set Enrichment Analysis

To explore the functional differences in gene expression between the high-risk and low-risk groups, we performed GSEA analysis. The results showed that gene expression in high-risk group was mainly related to ECM- receptor interaction, cancer signaling pathway and autophagy (Figures 6A and B).

Figure 6 (A) GSEA analysis in high-risk group:Gene expression in high-risk group was mainly related to ECM- receptor interaction, cancer signaling pathway and autophagy. (B) GSEA analysis in low-risk group.

Discussion

The mode of cell death has always been the focus of research. Since cancer cells have the ability to proliferate indefinitely, inducing cancer cell death has become a potential therapeutic direction.¹⁷ Programmed cell death (PCD) includes apoptosis, pyroptosis, and ferroptosis. Among them, pyroptosis has great potential in the treatment of cancer.¹⁸ Cell pyroptosis plays a dual role in cancer. On the one hand, we can reduce tumor load by inducing cancer cell pyroptosis; on the other hand, cell pyroptosis of cancer cells can lead to the expulsion of inflammatory substances, thus constructing the inflammatory microenvironment of tumor and providing conditions for tumor growth.⁷ Currently, the role of cell apoptosis in the physiology of cancer, including melanoma, remains unclear.

LncRNA is a class of non-coding RNAs with nucleotide numbers greater than 200. Known studies have revealed that lncRNA plays an important role in inflammation, tumor formation, and the acquisition of drug resistance.¹⁹ Many studies have confirmed that lncRNAs are involved in cell apoptosis, but the mechanism of their role in cell apoptosis in melanoma is poorly understood.

In this study, we explored the lncRNAs related to cell pyroptosis in melanoma and identified those lncRNAs with prognostic values. Through Cox regression and Lasso regression, we constructed a prognostic model consisting of 9 lncRNAs. By calculating the risk scores, 457 patients were assigned to a high-risk or low-risk group. Compared with the low-risk group, the survival rate was significantly lower in the high-risk group, with a statistically significant difference (P <0.001). ROC curves and calibration curves proved the accuracy of the prognostic model. And not only that, we found that the expression levels of immune checkpoint-related genes and M6A methylation-related genes were higher in the low-risk group.

For melanoma, a tumor with a poor prognosis, the establishment of an effective prognostic model and the exploration of its mechanism of action can improve the outcome of treatment. At present, immune checkpoint inhibitors and targeted inhibitors have achieved good therapeutic effects in melanoma. But these treatments can gradually build up resistance and bring regret to patients.²⁰ Cell pyroptosis is a promising new direction for the treatment of melanoma. GSDME, a key protein in the pyroptosis pathway, is not expressed in many cancer cells but is highly expressed in melanoma, which becomes the Achilles’ heel of melanoma.²¹ Induction of melanoma pyroptosis may play a potential role in improving the prognosis of advanced melanoma.²² Yu et al found that autophagy in melanoma is closely related to pyroptosis, a process mediated by EEF-2K.²³ Erkes et al also found that BRAF and MEK inhibitors can change the tumor microenvironment by inducing melanoma cell pyroptosis, thus improving the prognosis.²⁴ Therefore, it can be seen that there is a close correlation between cell apoptosis and melanoma, which has a great exploration prospect.

Melanoma is a highly immune-related tumor, and the discovery of immune checkpoint inhibitors is a milestone.²⁵ Clinical use of CTLA-4 and PD-1 inhibitors has benefited many melanoma patients, greatly improving disease-free survival and overall survival in melanoma patients.²⁶ Exploring the crosstalk between pyroptosis and immune checkpoints can further help us understand the clinical significance of apoptosis. According to our results, the expression of immune checkpoint-related genes was higher in the low-risk group, suggesting that immune checkpoint inhibitor therapy may be more beneficial in the low-risk group.

Epigenetics exists in a wide range of organisms, and posttranscriptional modification is a common regulatory pathway.²⁷ M6A-methylation is the most common post-transcriptional modification in eukaryotes.²⁸ This process is mediated by methylase (writer), demethylase (eraser), and methylated reading protein (reader), and is involved in the homeostasis regulation of organisms and the occurrence and development of diseases.²⁹ The regulation of M6A methylation in cells has also become a new idea for the treatment of diseases, especially tumors.²⁹ Through the correlation analysis of pyroptosis-related lncRNAs and M6A-methylation-related genes, we found that the methylation level in the low-risk group tended to be higher, which provided ideas for the future treatment of M6A methylation in melanoma. At the same time, it indirectly indicates that there is a certain degree of crosstalk between pyrosis and M6A methylation.

All in all, our study established a prognostic model of pyroptosis-associated lncRNA in melanoma patients and studied its mechanism of action. This is not only meaningful for guiding prognosis but also can facilitate the study of tumor physiology of melanoma, such as immune checkpoints and M6A methylation. However, our study has some limitations. We lack in vivo and in vitro experiments to verify the expression and function of lncRNAs in the prognostic model, and we will improve this aspect in the future.

Conclusion

Our study established a prognostic model of pyroptosis-related lncRNAs in melanoma patients to provide a reference for the study of pyroptosis in melanoma cells, and provide a new idea for melanoma treatment. But our study has some limitations due to the lack of validation.

Ethics Approval

All procedures performed were in accordance with the declaration of the ethical standards of the institutional research committee and with the 1964 Helsinki 387 Declaration and its later amendments. The ethics committee has approved this study of the First Affiliated Hospital of Nanjing Medical University (No. 2021-SR-418).

Acknowledgments

We are very grateful for data provided by databases such as TCGA, GEO. Jiaheng Xie, Haobo Li, Liang Chen, and Yuan Cao are co-first authors for this study.

Disclosure

All authors declare that no conflict of interest exists.

Additional information

Funding

There is no funding to report.