Identification of Pyroptosis-Relevant Signature in Tumor Immune Microenvironment and Prognosis in Skin Cutaneous Melanoma Using Network Analysis

Abstract

Background: Pyroptosis is closely related to the programmed death of cancer cells as well as the tumor immune microenvironment (TIME) via the host-tumor crosstalk. However, the role of pyroptosis-related genes as prognosis and TIME-related biomarkers in skin cutaneous melanoma (SKCM) patients remains unknown.

Methods: We evaluated the expression profiles, copy number variations, and somatic mutations (CNVs) of 27 genes obtained from MSigDB database regulating pyroptosis among TCGA-SKCM patients. Thereafter, we conducted single-sample gene set enrichment analysis (ssGSEA) for evaluating pyroptosis-associated expression patterns among cases and for exploring the associations with clinicopathological factors and prognostic outcome. In addition, a prognostic pyroptosis-related signature (PPRS) model was constructed by performing Cox regression, weighted gene coexpression network analysis (WGCNA), and least absolute shrinkage and selection operator (LASSO) analysis to score SKCM patients. On the other hand, we plotted the ROC and survival curves for model evaluation and verified the robustness of the model through external test sets (GSE22153, GSE54467, and GSE65904). Meanwhile, we examined the relations of clinical characteristics, oncogene mutations, biological processes (BPs), tumor stemness, immune infiltration degrees, immune checkpoints (ICs), and treatment response with PPRS via multiple methods, including immunophenoscore (IPS) analysis, gene set variation analysis (GSVA), ESTIMATE, and CIBERSORT. Finally, we constructed a nomogram incorporating PPRS and clinical characteristics to improve risk evaluation of SKCM.

Results: Many pyroptosis-regulated genes showed abnormal expression within SKCM. TP53, TP63, IL1B, IL18, IRF2, CASP5, CHMP4C, CHMP7, CASP1, and GSDME were detected with somatic mutations, among which, a majority displayed CNVs at high frequencies. Pyroptosis-associated profiles established based on pyroptosis-regulated genes showed markedly negative relation to low stage and superior prognostic outcome. Blue module was found to be highly positively correlated with pyroptosis. Later, this study established PPRS based on the expression of 8 PAGs (namely, GBP2, HPDL, FCGR2A, IFITM1, HAPLN3, CCL8, TRIM34, and GRIPAP1), which was highly associated with OS, oncogene mutations, tumor stemness, immune infiltration degrees, IC levels, treatment responses, and multiple biological processes (including cell cycle and immunoinflammatory response) in training and test set samples.

Conclusions: Based on our observations, analyzing modification patterns associated with pyroptosis among diverse cancer samples via PPRS is important, which can provide more insights into TIME infiltration features and facilitate immunotherapeutic development as well as prognosis prediction.

Figure 1

Landscape showing the genetic variations for 27 pyroptosis-regulated genes within SKCM. (a) Mutation frequency of pyroptosis-regulated genes in TCGA-SKCM samples. (b) Difference in prognosis between patients with the pyroptosis-regulated gene mutations and wild-type (WT) patients. (c) Enriched pathways in pyroptosis-regulated gene mutations and WT patients. (d) CNV frequency of pyroptosis-regulated genes in TCGA-SKCM samples. (e) Pyroptosis-regulated genes harboring CNV amplification displayed high expression within SKCM. (f) The differential expression in pyroptosis-regulated genes between metastatic and primary tumor patients.

Figure 2

Establishment of pyroptosis-associated profiles among SKCM cases as well as the relation with additional clinical variables and prognostic outcome. (a) Correlation analysis between the scores of pyroptosis-associated profiles and clinicopathological features of melanoma patients in TCGA-SKCM cohort. (b, c) Univariate (b) and multivariate (c) Cox regressions were used to analyze the independent prognostic factors for melanoma. (d) KM analysis between samples with high and low scores of pyroptosis-associated profiles. (e) The difference in scores of pyroptosis-associated profiles between different clinicopathological feature groups.

Figure 3

WGCNA on genes related to the pyroptosis-associated profiles among melanoma cases in TCGA-SKCM cohorts. (a) Clustering tree of each sample. (b) Analysis of the scale-free fit index for various soft-thresholding powers (β). (c) Analysis of the mean connectivity for various soft-thresholding powers. (d) Dendrogram of all differentially expressed genes clustered based on a dissimilarity measure (1-TOM). (e) Gene numbers within those 34 coexpression gene modules. (f) Heat map presenting the relations of modules with scores of pyroptosis-associated profiles. (g) Scatter diagram for module membership vs. gene significance in the blue module.

Figure 4

Construction of the PPRS for melanoma patients. (a) 984 prognostic genes were obtained by univariate Cox regression analyses. (b) The changing trajectory of each independent variable. The horizontal axis represents the log value of the independent variable lambda, and the vertical axis represents the coefficient of the independent variable. (c) Confidence intervals for each lambda. (d) Distribution of LASSO coefficients of PPRS.

Figure 5

Association of patients' survival with PPRS from both training and test cohorts. (a) 8 PAG expression levels, PPRS score, status, and survival time in each case from TCGA-SKCM dataset. (b) Performance of PPRS in predicting prognosis analyzed by 1-, 3-, and 5-year ROC curves in TCGA-SKCM dataset. (c) KM analysis between high and low PPRS score patients from TCGA-SKCM dataset. (d) Performance of PPRS in predicting prognosis in GSE65904 dataset. (e) Performance of PPRS in predicting prognosis in GSE54467 dataset. (f) Performance of PPRS in predicting prognosis in GSE22153 dataset.

Figure 6

The distribution of PPRS levels in samples with different clinicopathological features from TCGA-SKCM (a), GSE65904 (b), GSE54467 (c), and GSE22153 (d) cohorts.

Figure 7

The relation between pathways and PPRS. (a) Correlation analysis between PPRS and KEGG pathways through ssGSEA. (b) A heat map demonstrating normalized enrichment scores (NESs) of pathways in MSigDB calculated by comparing PPRS-high with PPRS-low groups (with a false discovery rate (FDR) of <0.05).

Figure 8

Relation between ability of tumor stemness and PPRS. (a) Differences in levels of mRNAsi between samples with diver risk scores. (b) Differences in levels of mRNAsi between samples in different pyroptosis statuses. (c) Correlation between mRNAsi and pyroptosis statuses. (d) Correlation between the levels of risk scores, GBP2, HPDL, FCGR2A, IFITM1, HAPLN3, CCL8, TRIM34, GRIPAP1, and tumor stemness.

Figure 9

Immune profiles of melanoma samples with distinct PPRS score. (a, b) Immune infiltration degrees within TIME for melanoma samples from TCGA-SKCM (a) and GSE65904 (b) datasets with high and low PPRS scores by CIBERSORT approach. (c, d) ESTIMATE/immune/stromal scores for melanoma samples from TCGA-SKCM (c) and GSE65904 (d) datasets with high and low PPRS scores. (e) The relation between immune infiltrations and PPRS.

Figure 10

Diverse immunotherapeutic responses among melanoma samples from PPRS-high and PPRS-low groups. (a–d) IC levels within melanoma samples from TCGA-SKCM (a), GSE65904 (b), GSE54467 (c), and GSE22153 (d) cohorts with high and low PPRS scores. (e) A bar plot demonstrating frequencies of immune checkpoints upregulated in PPRS-low cancer patients across the four cohorts. The y-axis indicates the names of immune checkpoints, and the x-axis represents the number of cohort. (f) The boxplots indicate the average immunophenoscore values IPS across the two PPRS subgroups in SKCM tumors. Overall, PPRS-low tumors that could be treated with combined anti-PD-1 and anti-CTLA-4 checkpoint blockade or with anti-PD-1 alone had significantly higher IPS, which is indicative of a better response to these immunotherapies. (g, h) The distribution of PPRS levels in samples responsive and resistant to treatment from GSE78220 (g) and GSE91061 (h) datasets and its performance in predicting samples prognosis. (i) Box plots exhibiting the estimated IC₅₀ values of temozolomide, paclitaxel, and cisplatin within melanoma samples from TCGA-SKCM dataset with high and low PPRS scores.

Figure 11

The nomogram and survival decision tree were produced for improving risk stratification and predicting the survival probability. (a) Cases who had complete annotations such as T stage, N stage, PPRS, and age were applied in constructing the survival decision tree for optimizing risk stratification. (b–d) Differences in OS (b), PPRS scores (c), and living states (d) were significant across the 6 risk groups. (e) Detailed information of the nomogram. (f) Our constructed PPRS and nomogram were highly accurate based on calibration analysis. (g) Decision-making curve of the nomogram. (h) Relative to additional clinicopathological factors, our as-constructed PPRS and nomogram performed well in predicting survival.