Construction of a Prognosis Model of the Pyroptosis-Related Gene in Multiple Myeloma and Screening of Core Genes.

Pyroptosis is an important factor affecting the proliferation, invasion, and metastasis of tumor cells. However, in multiple myeloma (MM), there are few studies on whether the occurrence of pyroptosis is related to the occurrence and prognosis of the disease. Based on the Gene Expression Omnibus and Cancer Genome Atlas database search dataset, this study identified pyroptosis-related genes with a specific prognosis, constructed and verified the prediction model by stepwise Cox regression analysis and time receiver operating characteristic curve analysis, and predicted specific functions by single-sample gene set enrichment analysis and the Kyoto Encyclopedia of Genes and Genomes. Dataset analysis identified key genes, which were used to construct a risk scoring system for the prognosis of MM. The entire test set and external verification set verified the results. The expression levels of related genes in the clinical samples were detected using fluorescence quantitative PCR. A prognostic gene model based on six pyroptosis-related genes (CYCS, NLRP9, AIM2, NOD2, CHMP3, and GSDME) was constructed. The model has an excellent prognostic ability and can be popularized in the external validation set. The predictive prognostic nomogram integrating clinical information can effectively evaluate the risk score of each patient and predict their survival. After sample validation, our study found three potential key pyroptosis-related genes in multiple myeloma. GSDME, NOD2, and CHMP3 were significantly different between MM and healthy subjects, suggesting that they are pyroptosis-related protective genes. This study shows that the key pyroptosis-related gene in the model can be used as a marker for predicting the prognosis of myeloma, which may provide a basis for clinical individualized stratification therapy.

Introduction

Multiple myeloma (MM) is a hematologic malignancy characterized by genetic diversity and the pathological proliferation of plasma cells.[1] The incidence rate of MM worldwide is 2.1/10 million, and the incidence rate in Latin American countries is higher (10/10 million).[2−5] With the application of proteasome inhibitors, immunosuppressants, monoclonal antibodies, small-molecule targeted drugs, and targeted B cell maturation antigen(BCMA)therapy, the depth of remission of patients with MM has achieved significant results, but it is still an incurable malignant tumor and its 5- and 10-year overall survival rates remain lower than 60% and 40%, respectively.[2] Therefore, improving the survival rate is currently an urgent challenge to be addressed, and a more effective prognosis evaluation system can provide personalized treatment for patients with MM to improve the prognosis of patients.

Pyroptosis is a new type of programmed cell death mediated by gasdermin.[6] Its essence is the formation of plasma membrane pores and release of inflammatory mediators (IL-1β, IL-18, etc.).[7] First, the innate immune system initiates the assembly and activation of inflammatory bodies. Gasdermin-D (GSDMD), as a key executive molecule, mediates cell death and further expands the inflammatory response.[8,9] The non-classical pathway is activated after oligomerization with caspase-11 in pyroptosis initiated by bacterial lipopolysaccharide(LPS) and changes in cell membrane permeability and current, promoting pyroptosis.[10] The cleavage of the GSDMD precursor into active fragments by activated caspase-related proteins is the core process of cell membrane damage via two pyroptosis-related pathways.[11,12] It has been found that the caspase-3-dependent pyroptosis pathway[13] and other proteases can also mediate pyroptosis.[14]

The relationship between pyroptosis and tumors is the focus of current research. A variety of solid cancer studies have reported relevant results. Chen et al.[15] found that small nucleolar RNA host gene 7 (snhg7) in the HepG2 hepatoma cell line inhibits pyroptosis by targeting the microRNA-34a/SIRT1 (nicotinamide adenine dinucleotide-dependent enzyme) axis and promotes cell proliferation, migration, and invasion. Knockout of the snhg7 gene can relieve the cell inhibition of caspase-1-dependent pyroptosis in hepatocellular carcinoma so as to aggravate the death of hepatocellular carcinoma cells and finally inhibit the growth of hepatocellular carcinoma in vivo. Tan et al.[16] showed that GSDME-mediated pyroptosis caused the DAMP (pathogen-related molecular pattern) proinflammatory mediator to release HMGB1 (high mobility group protein B1), which induced colorectal carcinogenesis and PCNA (proliferating cell nuclear antigen) expression through the MAPK/ERK1/2 signaling pathway. Shao et al.[17] pointed out that gasdermin-induced pyroptosis plays an important role in many genetic diseases, inflammatory diseases, and cancer, which provides an important suggestion for disease treatment.

In addition to solid tumors, a relationship between pyroptosis and hematological diseases has been reported. At present, it has been found to be involved in acute and chronic leukemia, myelodysplastic syndrome, and so on. Johnson et al.[18] confirmed that small molecule inhibitors of serine dipeptidase DPP8 and DPP9 induce pyroptosis in most human acute myeloid leukemia (AML) cell lines and primary AML samples and confirmed the inhibition of AML progression in mice, providing a new potential therapeutic strategy for AML. Fenini et al.[19] found that in immune and AML cells, the anticancer drug talabostat induces the activation of CARD8 (inflammatory body) and leads to focal death dependent on caspase-1 activation. Leu et al.[20] found that adicetone, a plant-derived natural component with an alkylbenzoquinone structure, can induce programmed death, such as leukemia cell apoptosis and pyroptosis, by downregulating IAP (inhibitor of apoptosis protein) and activating the caspase pathway. Sallman et al.[21] revealed the pathogenesis of MDS as being that the mutation of alarmin S100A9 and/or founder gene induces the pyrotosis of hematopoietic stem cells/hematopoietic precursor cells by producing reactive oxygen species, resulting in the ineffective hematopoiesis of patients with MDS.

There is evidence that pyrotosis plays a role in the occurrence and development of MM. Xia et al.[22] found that PRMT5 regulates pyroptosis by silencing CASP1 in MM. Zmorzyński et al.[23] found that the NOD2/CARD15 variant (3020insc) and PS6 polymorphism affect the occurrence and prognosis of MM. New insights into the genetic characteristics by Wang et al.[24] between MM and COVID-19 showed that NOD2 and COX6C genes could be used as prognostic biomarkers for MM. These studies suggest that there is a correlation between pyroptosis-related genes and MM, which affects the occurrence and development of MM and may provide evidence for exploring the unknown mechanism of MM. Here, we collected a large number of pyroptosis-related genes to explore the heterogeneity in MM and its impact on the tumor microenvironment. At the same time, we are committed to finding potential pyroptosis-related genes with prognostic significance in patients with MM and further screened key related genes in bone marrow samples, which may provide a new basis for diagnosis and treatment selection.

Results

Identification of Pyroptosis-Related Subtypes in MM

Based on the expression of 57 pyroptosis-related genes, we performed a clustering analysis of 696 patients in the TCGA-MMRF database and identified two clusters: cluster A included 328 patients and cluster B included 368 patients (Figure A). Subsequently, we conducted a survival study on patients in both clusters, and the findings indicated that cluster B had a significant survival advantage over cluster A (Figure B). When ISS scores were visually compared between the two clusters, patients with a better prognosis in cluster B (ISS = I) had a higher proportion of lower scores, whereas patients with a worse prognosis in cluster A had a higher proportion of higher scores (ISS = III) (Figure C), which is consistent with the prognostic trend and implied heterogeneity between pyroptosis subtypes. Given that prognosis is directly related to a patient’s sensitivity to drug therapy, we evaluated the sensitivity of two widely prescribed drugs in MM using the R package ″pRRophetic″ and unearthed no significant difference between the two clusters for lenalidomide (Figure D), but a lower IC50 for bortezomib (Figure E), denoting that cluster B patients had a higher drug sensitivity to bortezomib, which correlates with patient prognosis and suggests the value of our pyroptosis-related subtypes in predicting the response to drug therapy.

Figure 1

Identification and clinical characteristics of pyroptosis-related subtypes in MM. (A) Clustering of 696 patients in the database TCGA-MMRF. Cluster A has 328 patients while Cluster B has 368. (B) Kaplan–Meier curves showing the prognosis of patients in pyroptosis-related clusters, where the blue line represents cluster A and the red line represents cluster B; (C) differences in the ISS stage between cluster A and cluster B. Blue I represents low scores, orange II represents medium scores, green III represents high scores, and the area represents the percentage occupied; (D,E) differences in sensitivity to drug treatment between cluster A and cluster B. Blue represents cluster A, red represents cluster B; (D) differences in sensitivity to lenalidomide treatment; (E) differences in sensitivity to bortezomib treatment. Ns means “not statistically significant”; *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001 (all significance designations that appear in this paper are minor criteria).

Differences in Biological Pathways between Pyroptosis-Related Clusters

To investigate the differences between the clusters, we first analyzed the level of pyroptosis scoring in the two clusters using ssGSEA. We observed that cluster B had a higher level of pyroptosis scoring (Figure A) than cluster A, suggesting that pyroptosis may behave differently in distinct subtypes. We performed a pathway enrichment difference analysis of the two clusters to determine the possible biological implications of this difference. KEGG analysis revealed that cluster A was elevated for the ″PROTEASOME″ pathway, and DDR-related pathways, such as ″MISMATCH REPAIR″ and ″DNA REPLICATION″, which have been demonstrated to play a critical role in tumor growth, were also activated. However, cluster B exhibited activation of a wide range of immune modulation-related pathways, including ″ANTIGEN PROCESSING AND PRESENTATION″, ″TOLL-LIKE RECEPTOR SIGNALING PATHWAY″, ″CYTOKINE-CYTOKINE RECEPTOR INTERACTION″, and ″B CELL RECEPTOR SIGNALING PATHWAY″ (Figure B). Likewise, by analyzing the HallMark database, it was found that cluster A was enriched with pathways associated with cancer activation and pro-tumor proliferation, such as ″E2F TARGETS″, ″G2M checkpoint″, and ″MYC TARGETS″, while cluster B was enriched with immune-related pathways, such as ″COMPLEMENT″, ″IL6/JAK/STAT3 SIGNALING″, ″INFLAMM″, and ″INFLAMM″ (Figure C). The large biological pathway variations between clusters A and B were key contributors to the clinically positive differences between the two clusters. The predominance of tumor-friendly pathways in cluster A might give rise to a more aggressive phenotype, whereas cluster B subtypes exhibited high levels of immune activation, indicating that pyroptosis could be involved in defining the immunological microenvironment of MM.

Figure 2

Biological differences in pyroptosis-related clusters. (A) Difference between cluster A and cluster B pyroptosis scores. Blue indicates cluster A, red indicates cluster B. (B,C) GSVA enrichment analysis for cluster A and cluster B is demonstrated in a heatmap. Top bars represent different clusters, and different colored squares on the left correspond to the pathway annotations on the right. In the heatmap, a yellow color indicates up-regulated expression of the pathway and dark blue color indicates down-regulated expression of the pathway.

Differences in the Tumor Microenvironment between Pyroptosis-Related Clusters

A study of the distinctions in the pathways between the two clusters revealed that cluster B exhibited active activation of immunomodulatory-related pathways, indicating that the subcategories may have distinct immunological microenvironments. To begin, we compared the immune and stromal scores between the two clusters via ″ESTIMATE″, and, as seen in Figure , cluster B had a considerably higher immune score. Similarly, cluster B showed a large increase in stromal score, which seems to be related to the activation of pyroptosis (Figure A,B). Greater pyroptosis activation in cluster B than in cluster A evoked stromal and immunological component infiltration.

Figure 3

TME differences and immune checkpoints differences between 2 pyroptosis-related clusters. Blue represents cluster A and red represents cluster B. (A) Difference between stromal scores between cluster A and cluster B; (B) Difference between immunescore between cluster A and cluster B; (C) ssGSEA assessment of differences in the abundance of 21 infiltrating immune cells in two clusters; (D) Differential expression of immune checkpoints between the two clusters.

To further characterize immune cell infiltration, we used ssGSEA to compare the two subgroups of the 21 immune cells. Similar to the immune scoring findings, all cells were considerably more abundant in cluster B except activated CD8 T cell, CD56dim natural killer cell, and Type 2 T helper cell, which did not vary statistically between the two clusters (Figure C). Based on these findings, we concluded that cluster B featured a ″hot tumor″ profile, while cluster A exhibited a ″cool tumor″ profile. We then evaluated the expression of 10 immune checkpoints between the two clusters. CD276, CD86, CD96, CLTA4, and VTCN1 were not significantly different between the two clusters, whereas CD200, PDCD1LG2, PDCD1, CD274, and CD80 were significantly higher in cluster B (Figure D), indicating that cluster B may be a candidate group that could benefit from immune checkpoint inhibitors.

Analysis of Drug Sensitivity in Two Clusters

We used the R package ″pRRophetic″ to assess the sensitivity of many commonly used medicines in the two categories. The medications PF.02341066, PD.0332991, CHIR.99021, BX.795, and Bryostatin.1 showed a greater sensitivity in cluster B (Figure A–E), whereas the drug Embelin, DMOG A low was regarded as potentially beneficial in cluster A patients (Figure F,G). Notably, the targeted DNA repair agent AG.014699 demonstrated increased sensitivity in cluster A, with significant activation of the DDR-related pathways (Figure H).

Figure 4

Analysis for differences in drug sensitivity between two pyroptosis-related clusters. Blue represents cluster A and red represents cluster B. (A) Differences in sensitivity of the two clusters to PD.033291 treatment; (B) differences in sensitivity of the two clusters to CHIR.99021 treatment; (C) differences in sensitivity of the two clusters to BX.795 treatment; (D) differences in sensitivity of the two clusters to Bryostatin.1 treatment; (E) differences in sensitivity of the two clusters to PF.02341066 treatment; (F) differences in sensitivity of the two clusters to DMOG treatment; (G) differences in sensitivity of the two clusters to Embelin treatment; (H) differences in sensitivity of the two clusters to AG.014699 treatment.

Identification of Prognostic Genes Associated with Pyroptosis

We utilized the GSE24080 dataset as a training set and retrieved the expression of all 57 genes as well as survival time and status. First, we used univariate Cox regression to determine the association between gene expression and prognosis and then screened 16 genes associated with prognosis based on a threshold of p < 0.05. We then performed “lasso + stepwise Cox regression” to identify six genes with an independent prognosis: ″CYCS″, ″NLRP9″, ″AIM2″, ″NOD2″, ″CHMP3″, and ″GSDME″ (Figure ). ″CYCS″ and ″AIM2″ were considered to correlate unfavorably with prognosis, while ″NLRP9″, ″NOD2″, ″CHMP3″, and ″GSDME″ were recognized as a positive prognostic correlate.

Figure 5

Forest map of six gene multivariate Cox regression construction and validation of a prognostic model for pyroptosis-related genes.

In the training set GSE24080, we constructed a genetic prognostic model by screening six pyroptosis-related genes with independent prognostic value. We calculated independent prognostic scores for each patient using the following equation:

Construction and validation of a prognostic model for pyroptosis-related genes

Patients were then separated into two groups according to the best cutoff value, with those over the threshold being classified as high-risk (n = 189) and those below the threshold as low-risk (n = 365). We re-ran the survival study for GSE24080 using the updated risk classification and found that high-risk patients had poorer outcomes than low-risk patients (p < 0.0001) (Figure A). Figure B,C depicts the evolution of patient risk stratification and survival status as the risk scoring increases. We found that the prognostic factors ″CYCS″ and ″AIM2″ were expressed more in the high-risk group, while the protective genes ″NLRP9″, ″NOD2″, ″CHMP3″, and ″GSDME″ were expressed higher in the low-risk group (Figure D), which proved the reliability of our selected prognostic genes. A time-dependent ROC curve revealed prospective areas under the ROC curve (AUCs) of 0.66, 0.7, 0.76, and 0.69 at 2, 4, 6, and 8 years, respectively, indicating the strong predictive performance (Figure E). Additionally, we conducted external validation, and TCGA-MMRF found that the high-risk group had considerably lower survival rates than the low-risk group (Figure F), showing the six-gene model’s reliability.

Figure 6

Construction and validation of six pyroptosis-related genes prognosis models. (A) Survival differences between high- and low-risk groups in dataset GSE24080; (B) distribution of risk scores in MM patients in GSE24080; (C) survival status distribution of MM patients in GSE24080; (D) expression differences of six prognostic genes in high and low risk groups; (E) time–ROC curve analysis of the signature in training dataset; (F) Kaplan–Meier analysis of the validation dataset.

Construction and Validation of the Predictive Prognostic Nomogram

To further illustrate the predictive efficacy of the genetic prognostic model as a clinical prognostic tool, we combined it with clinical parameters in the form of dichotomous variables (high-risk and low-risk). The variables were screened by utilizing lasso regression, and the best prognostic factors, including ″AGE″, ″B2M″, ″LDH″, ″ALB″, and ″Risk group″, were identified by stepwise Cox regression based on the minimum AIC, and a prognostic nomogram was constructed to predict the 2, 4, 6, and 8 year survival rates based on the included factors (Figure A). AUCs for 2, 4, 6, and 8 year OS predictions were 0.77, 0.76, 0.81, and 0.87, respectively (Figure B). Additionally, we employed calibration curves to evaluate the model’s predictive capacity, as shown in (Figure C). The predictions nearly perfectly matched the best predictive performance, indicating that the model is capable of predicting patient prognosis. The alluvial diagram demonstrated the variance in terms of ″ISS stage″, ″risk group″, and cluster. As can be observed, the majority of high-risk patients constituted a larger proportion of cluster A patients, whereas patients in the lower-risk group constituted a higher proportion of cluster B patients (Figure D).

Figure 7

Construction of nomogram for predicting prognosis. (A) Nomogram for predicting survival after 2, 4, 6, and 8 years in MM patients; (B) time-dependent curve analysis for testing the ability of nomogram to predict the survival rate of MM patients after 2, 4, 6, and 8 years; (C) calibration plots of predicted nomogram. (D) Alluvial for observing the relationship between high- and low-risk groups, ISS stage, and pyroptosis-related clusters.

Clinical Sample Validation of Prognostic Genes

We detected the relative expression of genes in the experimental group and the control group (2—ΔΔCt), and the results show that the expression of GSDME, NOD2, and CHMP3 indicated a good prognosis in myeloma (p values were 0.0210, 0.0329, and 0.0039, respectively), while the expression of NLRP9 increased in MM (p = 0.0074). For the two genes with a poor prognosis, there was no difference in the expression of AIM (p = 0.0787) compared with the control group, whereas the expression of CYCs in MM decreased p = 0.0091) (Figure ).

Figure 8

Expression of real-time fluorescence quantitative PCR. Black represents MM, and red represents nomal control. (A–F) Expression levels of AIM2, CHMP3, CYCS, GSDME, NOD2, and NLRP9, respectively5. P values less than 0.05 were statistically significant.

Discussion

MM is a highly heterogeneous hematologic malignant disease with an extremely complex pathogenesis, and the mechanisms underlying its development remain largely unclear. Therefore, there is a need to explore disease-related biomarkers to differentiate patients with different prognoses and develop individualized treatment plans to ultimately improve patient prognosis. Pyroptosis is defined as cell death dependent on the pore-forming toxicity of gasdermin family proteins and is frequently accomplished by caspase activation, but not always.[25] In addition to classical and non-classical pathways, recent studies have focused on caspase-3-dependent pyroptosis pathways as well as on the direct induction of GSDME-dependent cancer cell pyroptosis by other substances.[14]

Transcriptomic analysis is a viable method for identifying tumor heterogeneity.[26] Using clustering analysis, we were able to identify two pyroptosis subtypes in MM, one of which had high levels of pyroptosis and the other had relatively low levels of the disease. We also discovered that the two subgroups had significantly different ISS scores, medication sensitivity, and survival rates, and that cluster B with high pyroptosis levels tended to have favorable clinical characteristics and prognosis, whereas cluster A demonstrated the reverse. According to these findings, pyroptosis in patients with MM has a distinct mechanism of action, which might exert an impact on medication responsiveness, hence the ultimate prognosis of patients.

Along with considerable activation of DDR-related pathways in cluster A, E2F TARGETS, the G2M checkpoint, MYC TARGETS, and other pro-cancer pathways were also upregulated. Dysregulation of DNA repair pathways in MM cancer cells has been shown to promote tumor resistance.[27] It has also been reported that MM cells overexpressing RECQ1, a decapping enzyme that plays a role in the repair of damaged replication forks, DNA damage response, and homologous recombination, were able to more efficiently repair DNA fragmentation caused by genotoxic agents, resulting in drug resistance.[28] Cluster A patients had several activations that supported tumorigenic growth and avoided drug killing-related pathways, which may engender more aggressive tumor cells and a primary drug-resistant phenotype as well as accelerated resistance to existing treatments. We observed that a substantial number of immune regulation-related pathways were activated in cluster B. In MM cells, the Toll-like receptor (TLR) is heterogeneous and often regarded as a ’double-edged sword’ in malignancy. Although TLRs have been shown to increase host protection against malignancies, cancer cells can exploit this pathway for their own benefit.[29,30] This effect can also be applied to MM. TLR activation is often associated with NF-kB activation, IFN secretion, and proinflammatory cytokine secretion, all of which promote the survival of patients with MM. However, several investigations have indicated that[31] TLR might initiate pyroptosis by binding to PAMP and activating the NF-kB pathway, thereby triggering the transcription of inflammasome components and downstream target substances. In our study, the activation of TLR was shown to be associated with elevated levels of pyroptosis in patients with cluster B disease, which seemed to be a protective mechanism. TLR-induced pyroptosis is involved in tumor suppression and improves patient prognosis.

In addition, we focused on the role of pyroptosis in the MM microenvironment. Immune evasion is a critical feature of cancer,[32] and we found that cluster A has a suppressive tumor microenvironment with a low abundance of immune and stromal components invading the tumor, which we believe represented ″cold tumor″. Cluster B, on the other hand, had a higher immune score, with an abundant microenvironment of CD4+ T, NK, NKT, and other immune cells, as well as increased expression of immune checkpoints CD200, PDCDLG2, PDCD1, CD274, and CD80, which tended to have a ″hot tumor″ phenotype. Without sensitization, NK cells can target infected and malignant cells, effectively eliminating tumor cells.[33] PD-L1 expression was elevated in T and NK cells isolated from patients with MM compared to healthy donors,[34] indicating that NK and T cells, which play a tumor-killing function in the myeloma microenvironment, might be repressed. In our study, cluster B with high pyroptosis levels had a significant infiltration of NK and T cells, whereas high CD4 (PD-L1) expression implied that these infiltrated cells could be dysfunctional. It is reasonable to speculate that the combined application of PD-1/PD-L1 inhibitors might be effective in cluster B, enhancing the antitumor effect and improving prognosis. CD4+ T cells can trigger efficient anticancer immune responses by interacting with antigen-presenting cells in the TME. Haabeth et al.[35] demonstrated that idiotype-specific CD4+ T cells induced therapeutic responses against MM in a mouse model, and subsequent studies revealed that MM created an environment in which CD4+ T cells interact with antigen-presenting cells, resulting in tumor cell death. Cluster B, which had a higher amount of CD4+ T-cell infiltration than cluster A, showed this interaction and exerted an active antitumor effect in our investigation. As a result, we hypothesized that the high level of pyroptosis in cluster B resulted in an inflammatory microenvironment capable of generating immune cells associated with a favorable patient prognosis, such as NK and CD4+ T cells, and that cluster B could serve as a potential target population for immunotherapy, which requires further investigation.

To aid the clinical translation of pyroptosis, we evaluated the sensitivity profiles of several FDA-approved anticancer medicines used in MM. Cluster B was not only susceptible to bortezomib but also to medications such as PF.02341066, PD.0332991, CHIR.99021, BX.795, and bryostatin.1. A phase II clinical trial was conducted to assess the safety and effectiveness of paboxib (PD.0332991) in combination with bortezomib and dexamethasone in patients with RMM, and the investigators noticed an objective response in 20% of patients and stable disease in 44% of patients. Cluster B patients with higher pyroptosis levels had greater sensitivity to paboxib in our study, and this group of patients might be potential beneficiaries of paboxib in combination with bortezomib and dexamethasone in sequential therapy. Additionally, we observed that the targeted DNA repair medication AG.014699 demonstrated increased sensitivity in cluster A while activating DDR-related pathways. AG014699 is a PARR inhibitor, and combination treatment with bortezomib and PARP1 inhibitors resulted in MM cell death.[36] Subsequently, a study[27] showed that these PRR inhibitors might be utilized in combination with DNA damage drugs to treat MM cells and reverse drug resistance. PARP inhibitors have been reported to be effective in reversing melphalan resistance in human myeloma cell lines, and their combination with melphalan exerted a synergistic effect on FA and homologous recombination DNA repair mechanisms. These results might bolster the evidence for the use of PARP inhibitors in patients with MM, and it is plausible to infer that patients with MM with low pyroptosis levels would have DDR activation, making them prospective PARP drug beneficiaries.

Additionally, we developed a prognostic gene model based on six pyroptosis-related genes (CYCS, NLRP9, AIM2, NOD2, CHMP3, and GSDME), which demonstrated high predictive prognostic performance and was generalizable via validation in an external validation dataset. A predictive prognostic nomogram, which incorporates clinical data, can efficiently analyze and forecast each patient’s survival.

Our study identified three potentially key pyroptosis-related genes in MM. GSDME is a member of the gasdermin family of proteins. It is expressed in many human tissues and organs and is closely related to tumors. Studies have shown that the growth and invasion of tumor cells increase after GSDME knockdown.[37] Studies have found that GSDME can lead to epigenetic silencing due to gene promoter methylation in many cancer cells.[38] Increasing evidence shows that epigenetic regulatory abnormalities (including abnormal expression of miRNAs, DNA methylation, histone methylation, acetylation, etc.) play an important role in the occurrence and development of MM.[39] In this study, it was found that GSDME was significantly different between MM and healthy individuals. Combined with biological information, it was suggested that GSDME is a protective gene, making it the first report on the relationship between GSDME and MM. NOD2 is an intracellular pattern recognition receptor that activates NF after activation of the κB pathway and triggers the innate immune system,[40] which plays a role as a suppressor gene in a variety of cancers.[41−44] In addition to inducing apoptosis and a focal response, NOD2 also promotes apoptosis.[45] Zmorzyński et al.[23] showed that the 3020 insc variant of the NOD2/CARD15 gene can lead to the upregulation of proinflammatory cytokines in patients with MM, and it was also found to be a positive prognostic biomarker of MM. In this study, NOD2 was significantly different between MM and healthy subjects. Combined with biological information, these results suggest that NOD2 is a protective gene. This conclusion is similar to that reported by Ski et al. CHMP3 is a member of the endosomal sorting complex required for the transport family. As a tumor susceptibility gene, CHMP3 participates in the EMT process of epithelial-mesenchymal transformation.[46] In this study, there were significant differences between the MM and healthy human CHMP3 gene. Combined with biological information, these results suggest that it is a protective gene, which has never been studied before.

With the deepening of research on pyroptosis and blood system diseases, we have acquired a deeper understanding of pyroptosis. It causes programmed cell death through a variety of activation pathways that are related to the occurrence and development of many diseases. Next, we will continue to further explore the biological functions of GSDME, NOD2, and CHMP3 as well as the possible mechanisms of the three prognostic genes, which will provide new ideas for the diagnosis and treatment of MM.

Conclusions

In summary, we explored two distinct pyroptosis subtypes by integrating the pyroptosis-related genes. The subtype with a low pyroptosis response showed low drug sensitivity to current conventional treatment. Compared to the subtype with a high pyroptosis response, it had more tumor-promoting effects, activation of related pathways, and worse prognosis. There were significant differences in the immune microenvironment between the two subtypes. The immune microenvironment of the subtype with a low pyroptosis level was poor, whereas the high pyroptosis response induced strong immune infiltration, manifesting as a hot tumor, which is a potential population that can benefit from immunotherapy. Simultaneously, we constructed a gene model to evaluate the risk score of a single patient and accurately stratify patient prognosis. After sample validation, we identified three potential key pyroptosis-related genes in MM. GSDME, NOD2, and CHMP3 were significantly different between MM and healthy subjects, suggesting that they are pyroptosis-related protective genes.

Materials and Methods

Data Download and Processing

R software version 4.1 was used to conduct all analyses in this study. A list of 57 pyroptosis-related genes was compiled based on previous studies (Supplementary Information 1). The Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases were employed to scan the dataset for pyroptosis function in MM. Ultimately, three datasets, GSE6477, GSE24080, and TCGA-MMRF-CoMMpass, were used. Raw CEL files were retrieved from the GEO database and the robust multiarray average method was used to normalize the two datasets. The GSE6477 dataset was obtained using the GPL96 platform (HG-U133A Affymetrix Human Genome U133A array). The dataset, which consisted of 15 healthy subjects, 22 MGUS (Monoclonal Gammopathy Of Undetermined Significance) patients, 24 patients with smoldering MM, 73 patients with newly diagnosed MM (NDMM), and 28 patients with refractory relapsed myeloma (RMM), was used to assess the expression of pyroptosis-related molecules in tumors versus non-tumors and in various stages of myeloma. GSE24080 was acquired using the GPL570 platform (HG-U133 Plus 2 Affymetrix Human Genome U133 Plus 2.0 array). The dataset included bone marrow samples from 563 patients with NDMM and a variety of clinical data, including survival time and status, ISS stage, sex, and age, as well as albumin, LDH, CRP, HGB, and B2M levels. The data in this dataset were cleaned using the following criteria: (1) patients with a survival time of less than 30 days were excluded; (2) patients with insufficient data on the aforementioned critical clinical features were excluded; and (3) 10 year survival was monitored, and if the OS time exceeded 10 years, we treated it as default survival and designated it as 10 years of OS time. Finally, 554 patients were included and, because of the wealth of clinical data included in this dataset, it was utilized as the training set while building the model.

The MMRF-CoMMpass dataset was retrieved from TCGA database and contained sequencing data for 716 individuals with NDMM as well as information on their survival time, survival status, and ISS scores. Additionally, the data in this dataset were cleaned using the following criteria: (1) samples must encompass entire survival data, including survival status and OS time, where death must be tumor-related and OS time must exceed 30 days; and (2) samples must contain complete R-ISS or ISS data. The final dataset, which contained 689 patients, was utilized as a subtype typing dataset and validation set for the genetic model.

Identification of the Pyroptosis-Related Clusters

To investigate the function of pyroptosis in MM, the expression of 57 pyroptosis-related genes was extracted from the TCGA-MMRF dataset and the R package ″ConsensusClusterPlus″ was applied to cluster 689 patients in TCGA-MMRF. A total of 1000 resamples were conducted to confirm classification reliability. Medication sensitivity was analyzed using the package ″pRRophetic,″ and then the patients were regrouped according to the clustering results, comparing patient survival, ISS staging, and drug sensitivity between groups.

Pathway Difference Analysis between Pyroptosis-Related Clusters

The changes in the pathways between clusters were investigated to better elucidate the heterogeneity of pyroptosis in MM. First, the KEGG (Kyoto Encyclopedia of Genes and Genomes) and Hallmark gene datasets were collected using MSigDB (Molecular Signature Database, http://www.gsea-msigdb.org/gsea/msigdb/index.jsp), and single-sample Gene Set Enrichment Analysis (ssGSEA) was then used to compute the score for each pathway in each patient. The R package ″limma″ was used to examine the differences in pathways between clusters, with the threshold values ″adj.P value 0.05″ and logFC absolute value of >0.1. The findings were shown using the R package ″pheatmap.″.

Immune Microenvironment Analysis

Since previous studies have shown the role of pyroptosis in forming the immune microenvironment of tumors, the immune microenvironment of patients with MM was examined in this study. To begin, the ″ESTIMATE″ package was used to examine the patients’ Immune Scores and Stromal Scores, and the discrepancies between the different pyroptosis subtypes were further compared to measure the overall immune and stromal infiltration between clusters. Immune cell infiltration was then further evaluated via ssGSEA. To eliminate interference, B, plasma, and associated cells were removed, leaving 19 immune cell types for investigation. After scoring the infiltration of each cell in each patient, the variations in infiltration between groups were compared using the pyroptosis subtype classification and are shown as a box line plot. Finally, we assessed the expression of 10 immune checkpoint molecules across the clusters(CD200, CD276, PDCD1LG2, CD86, PDCD1, CD96, CD274, CTLA4, VTCN1, and CD80).

Construction of a Prognostic Model for Pyroptosis-Related Genes

Survival data and 57 gene expression levels were gathered from the GSE24080 dataset after rearranging them according to predefined criteria. Risk ratios (HR) were calculated for 57 genes using univariate Cox regression, and a threshold of p < 0.05 was used to keep genes as possibly prognostically significant. Lasso-penalized Cox regression and stepwise regression analyses were performed to further filter for genes with the best prognosis. The aforementioned potential genes were incorporated in the construction of a multivariate Cox regression model. After computing the regression coefficients for each gene, the following equation was used to produce a risk score based on each gene’s expression:

(Exp gene 1 coefficient gene 1) + (Exp gene 2 coefficient gene 2) + ... + (Exp gene N coefficient gene N)

After calculating the risk score for each patient, the R package ″survminer″ was used to determine the optimal cutoff value, over which patients were classified as high-risk and those below as low-risk. Kaplan–Meier curves were used to compare the OS of both risk categories. The R package ″TimeROC″ was used to evaluate the predicted prognostic capacity of the models at different time points. Regression coefficients from the training set were utilized in conjunction with the full clinical data from TCGA-MMRF dataset to calculate the risk score of each patient for external validation.

Construction and Validation of a Prognosis Nomogram

To ascertain the independent prognostic value of the pyroptosis-related genetic prognostic model, lasso regression was utilized to screen prognostic factors by integrating widely available clinician prognostic features, and the factors identified by ″lambda.min″ were then subjected to stepwise Cox regression to determine the modeling factors based on the minimum AIC value. The prognostic performance of the nomogram was evaluated using the time-receiver operating characteristic (TIME ROC) and calibration curves.

Case Collection

To further screen hub genes, we verified the expression of six prognostic genes in patient bone marrow samples. By collecting the bone marrow tissue of 16 patients with MM treated in the Hematology Department of the Second Hospital of Shanxi Medical University from 2021 to 2022, eight bone marrow samples from patients diagnosed with iron deficiency anemia in the same period were selected as the control group. All patients with MM met the following enrollment criteria:[47] each patient was informed before collecting tissue samples and provided a signed written informed consent form. This study was approved by the hospital ethics committee of the Second Clinical Medical College of the Shanxi Medical University (code:(2020)YX(076)).

Experimental Verification of the Real-Time Fluorescence Quantitative Technique (qRT-PCR)

For quantitative real-time PCR, total RNA was extracted from the tissues using the TRIzol reagent (Invitrogen, Carlsbad, CA, United States). RNA was further reverse-transcribed into cDNA using Primer-Script RT Master Mix (RR036A, TaKaRa). qRT-PCR was performed using a SYBR green real-time PCR kit (Takara, Dalian, China) and FastStart Universal SYBR Green Master (ROX) (Roche, Germany). Primer synthesis was performed by Sangon Biotech (Shanghai, China). The employed primer sequences are presented in Table .

Table 1

Employed Primer Sequences

	primer-F	Tm	primer-R	Tm
CYCS	CTTTGGGCGGAAGACAGGTC	62.8	TTATTGGCGGCTGTGTAAGAG	60.1
AIM2	TCAAGCTGAAATGAGTCCTGC	60.4	CTTGGGTCTCAAACGTGAAGG	60.2
NLRP9	TTTGGCTTGTTGTGGTATCTGAA	60.3	CTGGGTAATGTTTGTCCAGCA	60.8
NOD2	CACCGTCTGGAATAAGGGTACT	60.9	TTCATACTGGCTGACGAAACC	60
CHMP3	ACCATGAGGGAGTTGTCCAAA	61.0	ACATCTCCTCTATGATCCCAGC	60.5
GSDME	TGCCTACGGTGTCATTGAGTT	61.4	TCTGGCATGTCTATGAATGCAAA	60.3

Statistical Analysis

The description of the data belongs to “mean ± standard error”, and the comparison between samples adopts the nonparametric test. The statistical significance was p < 0.05.