A novel pyroptosis related genes signature for predicting prognosis and estimating tumor immune microenvironment in lung adenocarcinoma.

Background: Pyroptosis is a newly found form of programmed cell death, accompanied by inflammatory response as well as immune response. Here, the specific function and prognosis predictive of pyroptosis-related genes (PRGs) were systematically explored in lung adenocarcinoma (LUAD).
Methods: The gene expression data and corresponding clinical information of LUAD patients were obtained from The Cancer Genome Atlas (TCGA), and the expression level of PRGs was identified between normal and tumor tissues. Furthermore, univariate Cox proportional hazards regression was conducted to filter the PRGs related to overall survival, and least absolute shrinkage and selection operator (LASSO) regression was subsequent employed to establish the PRGs risk model. Besides, the correlation of risk score with patients' clinical features, tumor mutational burden (TMB) as well as tumor microenvironment (TME) was also investigated.
Results: A total of 5 PRGs (NLRC4, NLRP1, NLRP3, NOD1, PLCG1, and BAK1) was used to establish the risk prognostic model. According the median value of risk score, all the patients were classified into low- and high-risk score group. Kaplan-Meier analysis indicted that the LUAD patients in low-risk group exhibited a better survival outcome compared the patients in high-risk group (P<0.001). After adjusting for age, gender, and clinical stage, the risk score was also considered as and independent risk factor affecting the overall survival of LUAD patients (HR =2.949, 95% CI: 1.762-4.937). Moreover, low-risk score group exhibited a higher Immune score and lower Tumor purity compared with high-risk score group. ssGSEA results proved that the enrichment scores of most immune cells and immune related signal pathway in low-risk score group was significant higher than that in high-risk score group. In addition, the PRGs risk score was also positive correlated with TMB in LUAD tissues. Conclusions: In this study, a novel prognostic model based on PRGs was constructed and used to predict the survival outcome of LUAD patients. In addition, the PRGs risk signature was also associated with TMB and anti-tumor immune environment. The induction of pyroptosis inside tumors might be considered a potential strategy in cancer treatments. 2022 Translational Cancer Research. All rights reserved.

Introduction

Lung cancer is the second most common cancer, and the most common deadly disease worldwide, with an estimated 2,220,000 new diagnoses, resulting in 1,800,000 deaths each year worldwide (1). Among the lung cancer patients, lung adenocarcinoma (LUAD) is the most common histologic subtype, accounting for approximately 40% of all patients. Most of LUAD patients die of metastatic, despite surgery, chemotherapy, targeted therapy and immunotherapy being used in cases treatment (2). So far, the clinical prognosis of LUAD mainly depends on tumor staging system (TNM). However, the TNM stage system is not sufficient to accurately predict the clinical prognosis of LUAD, especially to evaluate the effect of chemotherapy and immunotherapy. Therefore, it is imperative to establish a novel prognostic signature based on biomarkers for prognostic predicting and effect evaluation of cancer therapy.

Pyroptosis, also known as cellular inflammatory necrosis, is a newly found form of programmed cell death. Pyroptosis induced by various pathological stimuli, mainly mediated by the gasdermin family, accompanied by inflammatory response as well as immune response (3). Pyroptosis was initially found to be a key mechanism for combating infection, while recently a growing number of studies suggest that it also plays an important role in cancer initiation and progression (4-6). However, the exactly function of pyroptosis in malignant tumor is remain controversial.

Tumor microenvironments (TMEs) are composed of cellular components, extracellular matrix (ECM) and interstitial fluid. The long-term chronic inflammation in TMEs contributes to the progression and immunity activity of tumors. Hence, chronic inflammation induced by pyroptosis may facilitate tumor development, promoting the generation and maintenance of TMEs. In contrast, acute activation of pyroptosis results in the infiltration of immune cells, which could inhibit tumor development and progression (4). Gasdermin D (GSDMD) is the major executor of pyroptosis, which can be cleaved by caspases and leading cellular membrane pores. A recently study conducted by Gao reported that higher GSDMD expression is related to poor survival statue, as well as larger tumor size (7). Oppositely, Wang et al. reveal that low expression of GSDMD contributed to the occurrence and progression of cancer cell in gastric tumor (8). More importantly, growing evidence indicated that various tumor treatments such as chemotherapy, immunotherapy and other drugs could induce tumor pyroptosis, thereby inhibiting the malignant progression of tumors (5). Therefore, the trigger of pyroptosis might be considered a novel and effective strategy for cancer therapy.

In the current study, we aim to explore the correlation of PRGs with clinical prognostic and tumor immune microenvironment in LUAD. Thus, our study provides a novel understanding of the role of pyroptosis in LUAD, suggesting that pyroptosis might be considered as a potential strategy for clinical prognostic and effect of immunotherapy assessment. We present the following article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-327/rc).

Methods

Data acquisition and processing

Transcriptome profiling data (FPKM) of LUAD patients were obtained from the TCGA data portal on 18 Sep 2021 (https://portal.gdc.cancer.gov/repository). A total of 497 lung cancer tissues and 54 normal lung tissues were included in the present study. The corresponding clinical information of LUAD patients was also obtained from the TCGA data portal, and the patients without survival data were removed from further analysis. To identify the different expression level of PRGs between cancer tissues and adjacent normal tissues, the “limma” package was employed with a P value <0.05 considered statistically significant. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification of the Prognostic Gene Signature based on PRGs

A total of 39 PRGs were extracted from previous research (9-11), which is shown in Table S1. To further establish the clinical prognostic model based on PRGs, the “survival” package of R language software was employed. The univariate Cox proportional hazards regression analysis was initially used to screen the PRGs, which were significantly correlated to the survival status of LUAD sample. Furthermore, LASSO Cox regression analysis was subsequently conducted to establish the risk model using the “glmnet” package of R language. Ultimately, the 5 genes and their coefficients were retained. The risk scores of each LUAD patient were calculated using the following equations: risk score= [X: coefficients (retained by LASSO Cox regression analysis), Y: gene expression level]. According the median value of the risk score calculated by above formula, all the samples were classified into high- and low-risk score groups in TCGA cohort. Moreover, ROC curve was used to access the sensitivity and specificity of the risk score model performance using the R “survival” package. In addition, the correlation of risk score and overall survival was also investigated by the univariate and multivariate Cox proportional hazard regression analysis using the “survival” package of the R language. The risk factors of Cox proportional hazard regression analysis included the risk score, age, gender, and tumor clinical stage.

Functional enrichment analysis of the DEGs between the low- and high-risk groups

The entire LUAD samples were classified into high- and low- risk group, according to the median value of risk score in TCGA LUAD cohort. First, the differentially expressed genes (DEGs) between low- and high-risk score group were screened by “limma” package of R language with a criteria of (|log2FC| ≥0.585 and FDR <0.05). Second, Kyoto Encyclopedia of Gene and Genome (KEGG) and gene ontology (GO) analysis was subsequent conducted by “clusterProfiler” package of R language software. Ultimately, gene set enrichment analysis (GSEA) was performed to clarity the signal pathway, which were significantly changed between high- and low-risk score groups.

Tumor mutational burden and tumor immune microenvironment analysis

The mutation data of LUCA patients were downloaded form the TCGA data portal and analysed by “maftools” package of R language software. The correlation of TMB and risk score was estimated by Spearman correlation test using “ggplot” package of R language software. Furthermore, The ESTIMATE score, immune score, stromal score and tumor purity of each LUCA samples were calculated using R package “estimate”. Besides, CIBERSORT was also performed to evaluate relative infiltrated proportion of immune cells by the ‘CIBERSORT’ R package. In addition, the enrichment score of 16 immune cells and 13 immune-related pathways were also calculated using the single-sample gene set enrichment analysis (ssGSEA) method in the Gene Set Variation Analysis (GSVA) package of R language software.

Statistical analysis

All the analyses were conducted using the R language software (R version 4.1.0 for Mac). Moreover, Wilcoxon test was conducted for two groups comparison, while Kruskal-Wallis test was conducted for more than two groups. The risk factors affecting the overall survival of LUCA patients were evaluated by the univariate and multivariate Cox proportional hazard regression analysis using the “survival” package of the R language. The area under the ROC curve (AUC) was used to test the capability of prognostic accuracy. In all analyses, P<0.05 were considered statistically significant.

Results

Identification of DEGs between normal and tumor lung tissues

First, the different expression level of 39 PRGs between 54 normal and 497 tumor tissues in TCGA LUAD samples was compared. As illustrated in , it was showed that 33 PRGs were significantly differentially expressed with P value <0.05 between normal and tumor tissues (). Among the 33 PRGs, the expression level of NLRC4, CASP5, IL6, IL1B, ELANE, NLRP3, CASP1, TNF, NLRP1, IRF1, IL18, NOD1, PYCARD, GZMB, HMGB1, and IRF2 were down-regulated, while the expression of GSDMD, GPX4, CASP4, PLCG1, TIRAP, BAX, CASP8, CASP3, BAK1, GSDME, CASP6, PJVK, GSDMA, GSDMB, NLRP7, AIM2 and GSDMC were up-regulated in tumor samples. The expression pattern of the PRGs was presented in .

Figure 1

Different expression level of PRGs between normal and lung tumor tissues. (A) Different expression level of PRGs between normal tissues (n=54) and LUAD tissues (n=497). P values were notated as follow: *, if P<0.05, **, if P<0.01, and ***, if P<0.001. (B) Heatmap of different expression pattern of PRGs in normal and LUAD tissues. PRGs, pyroptosis related gene; LUAD, lung adenocarcinoma.

Establishment of prognostic gene model based on PRGs

To investigate the association between pyroptosis and clinical prognostic in LUAD sample, univariate Cox regression analyses were initially performed to filter the PRGs significantly related to overall survival outcomes. As a result, the expression level of NLRC4, NLRP1, NLRP3, NOD1, PLCG1, and BAK1 were found to have significant association with overall survival for LUAD patients. Interesting, all the 6 PRGs, except BAK1 (hazard ratio >1), were considered as protecting factors (hazard ratio <1) (). Furthermore, LASSO Cox regression analysis with minimized lambda was subsequently conducted to establish the risk score model. The risk scores of each LUAD patient were calculated using the follow formula: risk score=(NLRC4*-0.2827)+(NLRP1*-0.1227)+(NOD1*-0.1570)+ (PLCG1* -0.1030)+ (BAK1* 0.2021) (, ). According the median value of risk score calculated by the above formula, all the patients were classified into low- and high-risk score group. The AUCs of PRGs risk score for 1, 3 and 5 years were 0.678, 0.610 and 0.634 respectively (). Moreover, Kaplan-Meier analysis indicted that the LUAD patients in low-risk group exhibited a better survival outcome compared the patients in high-risk group (, P<0.001). In addition, as shown in , the high-risk group exhibited a higher death rate compared with low-risk group. Consisted with the above result, the principal component analysis (PCA) showed that patients with different risk score were well separated into two clusters ().

Figure 2

Construct of risk model based on PRGs in TCGA LUAD cohort. (A) Forest plots for the results of the univariate Cox regression analysis between overall survival and PRGs expression in LUAD patients. (B) LASSO regression of the 5 overall survival related PRGs. (C) ROC curves tested the specificity and sensitivity of the risk score model for 1, 3, and 5 years. (D) The survival curves for LUAD patients in high- and low-risk score group. PRGs, pyroptosis related gene; LUAD, lung adenocarcinoma; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; AUC, area under the curve; TCGA, The Cancer Genome Atlas;.

Table 1

Coefficients in the LASSO Cox regression model

No.	Gene	Coef
1	NLRC4	−0.2827
2	NLRP1	−0.1227
3	NOD1	−0.1570
4	PLCG1	−0.1030
5	BAK1	0.2021

LASSO, least absolute shrinkage and selection operator.

Figure 3

Correlation of risk score and survival status in LUAD patients. (A) Distribution of risk score for each LUAD patients with different risk score. (B) Distribution of survival status for each LUAD sample with different risk score. (C) Heatmap for the expression pattern of 5 PRGs and clinical features between high- and low-risk score group. (D) PCA plot for each LUAD samples with different risk score. (E) Univariate Cox regression analysis for the overall survival of LUAD patients. (F) Multivariate Cox regression analysis for the overall survival of LUAD patients. LUAD, lung adenocarcinoma; PRGs, pyroptosis related gene; PCA, principal component analysis.

Correlations of risk score and LUAD samples’ clinical features

To further explore the role of PRGs risk score on the prognosis of LUAD patients, univariate Cox regression and multivariable Cox regression analysis was employed to identify the risk factors that affecting the patients’ overall survival. As shown in , univariate Cox regression analysis results indicated that clinical stage and risk score were potential factors affecting the patients’ prognosis. More important, after adjusting for age, gender, and clinical stage, the risk score was also considered as an independent risk factor affecting the overall survival of LUAD patients (HR =2.949, 95% CI: 1.762–4.937), which was shown in . Moreover, the PRGs risk score was also positive correlated with clinical stage, T and N stage of LUAD patients, while risk score was not significant related to patients’ Age, Gender and M stage ().

Figure 4

Boxplots for the correlation between PRGs risk score and clinical features in LUAD patients. (A) Patients’ age at diagnosis. (B) Patients’ gender. (C) Tumor clinical stage. (D) Tumor T stage. (E) Tumor N stage. (F) Tumor M stage. PRGs, pyroptosis related gene; LUAD, lung adenocarcinoma.

Functional analyses of DEGs between high- and low- risk score groups

To further explore the differences in the gene functions and pathways between low- and high-risk score group, GO, KEGG, and GSEA were employed. With the screen criteria of (|log2FC| ≥0.585 and FDR <0.05), the DEGs between high- and low-risk score group were mainly enriched in “neutrophil activation”, “neutrophil mediated immunity”, “neutrophil degranulation”, “neutrophil activation involved in immune response” and “mitochondrial inner membrane”, etc. in GO analysis (P<0.05, ). In addition, KEGG pathway analyses indicated that the DEGs were highly enriched in “Phagosome”, “Staphylococcus aureus infection”, “Hematopoietic cell lineage”, “Cell adhesion molecules”, “Tuberculosis”, etc. (P<0.05, ). Moreover, with the GSEA results for the DEGs enrichment, the gene sets of the high-risk score group were enriched in “Alzheimers disease”, “Huntingtons disease”, “Oxidative phosphorylation”, “Parkinsons disease” and “Ribosome” (), while the gene sets of the low-risk group were enriched in “Allograftrejection”, “Asthma”, “Cytokine cytokine receptor interaction”, “Intestinal immune network for IgA production”, “Systemic lupus erythematosus” ().

Figure 5

Functional analysis in low- and high-risk score group. (A) The top 5 enrichment in biology process, cellular components and molecular functions for DEGs between high- and low-risk score group. (B) The top 5 enriched KEGG pathways for DEGs between high- and low-risk score group. (C) Gene sets enriched pathways of GSEA in high-risk score group. (D) Gene sets enriched pathways of GSEA in low-risk score group. DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Gene and Genome; GSEA, gene set enrichment analysis; BP, biological progress; CC, cellular component; MF, molecular function.

Correlations of risk score and TMB

As a well-known evaluation marker of tumor immunotherapy efficacy, we also estimate the correlation between the PRGs risk score and TMB in LUAD patients. The gene mutation for LUAD patients between high- and low-risk score group was shown in , with the top 20 most frequently mutated genes being followed: TP53, TTN, MUC16, RYR2, CSMD3, LRP1B, ZFHX4, USH2A, KRAS, XIRP2, FLG, SPTA1, NAV3, ZNF536, COL11A1, FAT3, PCDH15, CSMD1, ANK2, and KEAP1. Moreover, as illustrated in , TMB was found to be significant higher in high-risk score group than low-risk score group (P=0.0015), In addition, the risk scores were also positive correlated with TMB in LUAD tissues (, r=0.19, P<0.01).

Figure 6

Correlation analysis of PRGs risk score and TMB in LUAD patients. (A) Waterfall showing the top 20 mutated genes in low-risk score group. (B) Waterfall showing the top 20 mutated genes in high-risk score group. (C) Boxplots for the TMB of LUAD tissue in low- and high-risk score group. (D) The correlation with risk score and TMB in TCGA LUAD cohort. PRGs, pyroptosis related gene; TMB, tumor mutational burden; LUAD, lung adenocarcinoma; TCGA, The Cancer Genome Atlas.

Correlations of risk signature with TME and immune activity LUAD

As pyroptosis is often accompanied by inflammatory release and immune response, the correlation between the PRGs risk score and TME in LUAD patients was also identified. As illustrated in , our results revealed that low-risk score group displayed a higher ESTIMATE score, Immune score, and Stroma score compared with high-risk score group (, all P<0.01), while the patients with high-risk score exhibited a higher Tumor purity compared with low-risk score group (, P<0.01). Furthermore, the correlation of PRGs risk score and tumor immune microenvironment was also investigated by ssGSEA. As illustrated in , except NK cells and Th2 cells, the other enrichment scores of 16 types of immune cells in low-risk score was significant higher than that in high-risk score group. Meanwhile, with the enrichment scores of immune function analysis, except MHC class I, the other enrichment scores in low-risk score group was significant higher than that in high-risk score group.

Figure 7

Violin plots for the difference of (A) ESTIMATE score, (B) immune score, (C) stromal score, and (D) tumor purity in low- and high-risk score group in LUAD patients. LUAD, lung adenocarcinoma.

Figure 8

Correlation analysis of risk score and ssGSEA enrichment in LUAD patients. (A) Comparison for the ssGSEA scores of immune cells between high- and low- risk score group in TCGA LUAD patients. (B) Comparison of the ssGSEA scores of immune-related pathways between high- and low-risk score group TCGA LUAD patients. (C) The heatmap for the ESTIMATE score, Immune score, Stromal score, tumor purity and ssGSEA scores of immune cells between low- and high-risk score group in TCGA LUAD patients. (D) The heatmap for the ESTIMATE score, Immune score, Stromal score, tumor purity and ssGSEA scores of immune-related pathways between low- and high-risk score group in TCGA LUAD patients. *, P<0.05, **, P<0.01, and ***, P<0.001. TCGA, The Cancer Genome Atlas; LUAD, lung adenocarcinoma.

In addition, the relative infiltration of 22 types of immune cell between high- and low-risk score group was also explored by CIBERSORT. As illustrated in , The infiltrated level of T CD4 memory resting, Monocytes, and Dendritic resting cells in low-risk score group was significant higher than that in high-risk score group. Meanwhile, high infiltration levels of T CD4 memory resting, Monocytes and Dendritic resting cells were significantly correlated with better OS outcomes in LUAD patients.

Figure 9

Immune cell infiltration analysis in TCGA LUAD tissues with different risk score. (A) Overall view of relative infiltrations for 22 types of immune signatures. (B) Boxplots for different immune cell infiltrations in high- and low-risk score groups. (C) Kaplan-Meier curves for overall survival between the high- and low-risk group in immune cells of T CD4 memory resting, Monocytes, and Dendritic resting cells. TCGA, The Cancer Genome Atlas; LUAD, lung adenocarcinoma.

Discussion

Pyroptosis is a novel type of PCD that plays a critical role in both septic shock and immune defenses, discovered after apoptosis and necrosis. Pyroptosis is characterized by cellular swelling, pore formation in membrane, cell lysis and release of pro-inflammatory mediators, including IL-1β, IL-18, and HMGB14, which induce inflammatory responses (12). At present, pyroptosis is reported to participate in the occurrence and development of various diseases, especially its dual role in promoting and inhibiting tumor formation and TME (6). It was found that the expression of GSDMB, an executor of pyroptosis, was down-regulated in normal tissue as compared to BLCA tissue, while overexpressed GSDMB facilitated tumor progression. It can be explained by the fact that multiple signaling pathways and inflammatory mediators are released during pyroptosis, which contributing to tumor growth and progression (13). In contrast, a recent study conducted by Wang WJ revealed that GSDMD may work as a treatment for gastric cancer by inhibiting cell proliferation, and it also may be used as a diagnostic tool (14). It can be explained that pyroptosis suppress tumor incidence and progression, which could become a potential treatment for drug-resistant cancers in the future (15). However, the specific role of pyroptosis on the development and progression in LUAD is remained unclear.

In the current study, the different expression level of PRGs between adjacent normal and tumor tissues was identified. The results indicated that most of PRGs expression levels were significant different between LUAD and adjacent normal tissues in the TCGA LUAD cohort, thus indicated the process of the process of pyroptosis may involve in the oncogenesis, but the explicit molecular mechanism remained need to further study.

Since pyroptosis may play an opposite role on promoting and inhibiting cancer progress in different type of cancer, we further identify the exactly role on the clinical prognosis of LUAD patients. As a result, we identified six PRGs that related to the LUAD patients clinical prognosis outcome. Interesting, except BAK1, the other PRGs including NLRC4, NLRP1, NLRP3, NOD1 and PLCG1 were considered as protective factors. The NLRC4, NLRP1 and NLRP3 belong to Nod-like receptor (NLR) family, consisted of inflammasome complex to recruit caspase-1 and promote its proximity-induced activation to induce pyroptosis (16). Among the above inflammasomes, NLRP3 is the mostly characterized inflammasome, activated by bacterial toxins and secretion system components, pathogenic crystals, and altered cellular components (6). It was reported that the NLRP3 inflammasome also participated in the innate immune response to some types of tumors (17). Consistent with our research, a recent study conducted by Wei et al. found that the expression of NLRP3 in HCC tissues was significantly down regulated or even completely absent, and its expression was negatively correlated with the pathological grade and clinical stage of HCC, indicating that the NLRP3 inflammasome was involved in the progression of HCC (18). Moreover, Nadatani et al. reported that the expression of NLRP3 in Barrett’s esophageal cancer cells treated with LPS was increased, and pyroptosis were also increased. These results indicated that the NLRP3 inflammasome activation of caspase-1 induces the secretion of pro inflammatory factors and pyroptosis (19). It was reported that blocking NLRP3 activation using MCC950, a molecular inhibitor of NLRP3, suppressed tumor growth in head and neck squamous cell carcinoma, accompanied by decreased immunosuppressive cell accumulation and increased the number of effector T cells. Thus, inhibition of the TME through the NLRP3 inflammasome might provide a novel approach for tumor therapy (20). However, to knowledge, there was no related study has been reported on the effect of MCC950 in the treatment of LUAD, and further research is needed.

In contrast to NLRP3, NLRP1 has not been precisely analyzed. Most studies have proven that NLRP1 prevent the interaction of apoptotic caspases with other activating proteins, thus negatively regulating apoptosis and promoting tumor progression (21). Oppositely, a study conducted by Chen et al. proved that lower expression level of NLRP1 were correlated with higher clinical stage and shorter patient survival period (22), which consist with our study. In macrophages, NLRC4 inflammasomes could be activated by inflammatory infections and endogenous stimuli, thus inducing pyroptosis (23). The effects of the NLRC4 inflammasome on tumor occurrence and progression are still controvertible. It was reported that NLRC4 could stimulates the activation of macrophages and enhance the production of IFN-γ in T cells, suggesting that NLRC4 might lead to the increase tumor immune microenvironments, thus playing an crucial role in tumor inhibition (24).

TMEs have been suggested to play an important role in tumor progression. Growing studies revealed that the immune reaction and inflammatory cytokines released during pyroptosis, including IL-1β, IL-18, ATP, and HMGB1, could be exert an important influence on the TME (16,25). It was reported that IL-1β signaling could induce DC maturation and monocyte differentiation into DCs and inflammatory macrophages. IL-18 exerts an important effect on natural killer (NK) cell recruitment and activation, as well as Th-1 polarization. Zhang et al. (26) also found that tumors with high expression of the wild type GSDME display increased levels of immune cell infiltration, including CD8+ T cells and natural killer (NK) cells. Consist with the previous study (25), our results also indicated pyroptosis induction in LUAD tissues increased the infiltration of immune cells and activity of immune related signal pathway, thus repressing tumor progression and development.

However, there are several limitations in our study. It is important to note that all data in the current study were obtained from online public databases and did not involve any experiments. In addition, the performance of the PRGs model was not verified in another independent cohort.

In conclusion, our study revealed that pyroptosis was closely related to the clinical prognosis of LUAD patients. Moreover, a novel prognostic model based on 5 PRGs was constructed for survival predicting. Besides, the PRGs risk score was closed related to anti-tumor immune environment. The induction of pyroptosis inside tumors might be considered a potential strategy in cancer treatments. However, the beneficial and detrimental effects of pyroptosis need to be further investigated in clinical settings.

The article’s supplementary files as