A New Risk Score Model Based on Pyrogenic Signatures for the Prediction of Survival and Immune Microenvironment Features in Lung Adenocarcinoma

Introduction

Currently, lung cancer remains the leading cause of cancer-related deaths worldwide.¹ Lung adenocarcinoma (LUAD), as one of the two non-small cell lung cancer (NSCLC), accounts for about 40–50% of all lung cancer cases.² The latest progress in targeted therapy and molecular pathology has brought some benefits to LUAD patients. However, LUAD is prone to metastasis in the early stage. The prognosis of LUAD patients is poor, with the average 5-year overall survival (OS) rate less than 20%.³ Further understanding of the molecular mechanisms of tumor genesis and progression in LUAD could enhance the overall prognosis and treatment outcome of LUAD.

Pyroptosis (also known as cell inflammatory necrosis) is a new type of programmed cell death,⁴ which characterized by cell swelling and multiple bubble-like swells processes. During this process, the cells first form multiple vesicles, followed by the formation of small holes on the cell membrane, ultimately leading to the rupture of the pores and the outflow of cellular material.⁵ The pyroptosis in cells depends on the inflammatory caspase family and the release of large numbers of pro-inflammatory cytokines.⁶ Pyroptosis plays important roles in the development of tumors.^7,8 This potent pro-inflammatory effect is also associated with the regulation of the tumor immune microenvironment (TIME). Xi et al indicated that expression deficiency of gasdermin D (GSDMD), one of the key pro-inflammatory cytokines of pyroptosis, was accompanied by the significant decrease in the activity and number of CD8⁺ T lymphocytes.⁹ In addition, research has confirmed that pyroptosis plays a crucial role in the antitumor function of NK cells.¹⁰ Hence, pyroptosis-related genes may be able to predict the prognosis of LUAD and provide guidance for treatment.

Pyroptosis plays an important role in the development of tumors and anti-tumor processes. Few studies have reported its specific functions in LUAD. In this study, bioinformatics methods were used to investigate the expression levels of pyroptosis-related genes in normal and tumor tissues in COAD, so as to explore the prognostic value of these genes and the correlation between pyroptosis and TIME.

Materials and Methods

Datasets and Preprocessing

Three datasets were used. The dataset in the experimental group was downloaded from UCSC Xena. RNA sequencing data (FPKM value) and related clinical information of gene expression of LUAD was obtained and annotated. Mutation data of LUAD were obtained from the Cancer Genome Atlas (TCGA) database. LUAD data (oncosg, NAT Genet 2020)¹¹ from cbioportal were downloaded for the verification analysis. The GSE61676 dataset used for the immunotherapy queue is from the Gene Expression Omnibus database (GEO) database. The following data were removed for further analysis: 1) samples without clinical follow-up information; 2) samples with unknown survival time, survived less than 0 days, and unknown survival status. The gene symbol was converted to obtain expression data. The probe linked to multiple genes was removed. An average value was used for expression of multiple gene symbols. Totally, 59 normal samples and 513 LUAD samples were included in the experimental group. One hundred and sixty-three LUAD samples were used for the validation analysis and 43 non-small cell lung cancer (NSCLC) samples were included in the immunotherapy cohort. Thirty-three pyroptosis genes were all referenced from the literature.¹²

The study was conducted in accordance with the Declaration of Helsinki and the obtained of tissue sample was approved by the ethics committee of the First Affiliated Hospital of Dalian Medical University. The informed consent was obtained from each patient and their families for clinical data using and publication.

Pyroptosis-Related Risk Score Model

Correlation analysis Pearson of pyroptosis-related genes was performed in the experimental group to obtain the risks core model. The univariate Cox algorithm was used to reduce the noise or redundant genes of pyroptosis gene set. The variables were further screened by the least absolute shrinkage and selection operator (Lasso mentioned by Tibshirani (1996))¹³ to reduce the number of genes. Finally, the multi-factor Cox regression model was applied to construct the risk score model. The calculation formula was risk scores=Σ Coef (i)*Exp (i). Patients were divided into high-risk group and low-risk group according to the median risk score. The impact of characteristic genes on survival in the risk score was analyzed by the Kaplan Meier method. The survival curve and time-dependent receiver operating characteristic (ROC) curve of the risk score were used to verify the accuracy of the signature, and further verification was conducted.

Identification of Differentially Expressed Mutant Genes

The Wilcox test was utilized to analyze the differentially expressed characteristic genes in cancer and adjacent samples of risk score. The human protein atlas was used to find different protein levels of characteristic genes in the database. The maftools function was applied to analyze the mutation data of characteristic genes. The transcription factor (TF) subsets were obtained from Cistrome, and the corresponding TF expression values were obtained from the TCGA database. TFs associated with characteristic genes in the risk score were identified by correlation analysis (Pearson correlation coefficient >0.3, p-value <0.05).

RNA Extraction and Real-Time PCR

TRIzol method was used to extract total RNA from normal tissues and tumour tissues of 5 patients with LUAD. According to the protocol of the manufacturer (TIANGEN, BEIJING), cDNA was obtained by reverse transcription using FastKing cDNA first-strand Synthesis kit. Relative gene expression in LUAD tissues was detected using SuperReal PreMix Plus (SYBR Green). The PCR primers were as follows: NLRC4 forward 5-CTCATGAGCCAGAGGAGGTG-3, reverse 5-ATACACTGCTGCGAGGTGCT-3; NLRP1 forward 5-AGATGGGGCTGAGGTCACT-3, reverse 5-CGCCGGCAATTCATGGATC-3; NOD1 forward 5-TCCAAGTTCGTGCTGTGCT-3, reverse 5-ATTGAGGATGCCGGTGGTG-3; PLCG1 forward 5-GGAAGACCTCACGGGACTTTG-3, reverse 5-GCGTTTTCAGGCGAAATTCCA-3. GAPDH forward 5-GGAGCGAGATCCCTCCAAAAT-3, reverse 5-GGCTGTTGTCATACTTCTCATGG-3 were used as internal genes. The relative gene expression was determined by the 2^−ΔΔCT method.

Cell Infiltration Evaluation

A single sample gene set enrichment analysis (ssGSEA) algorithm was used to quantify the relative abundance of each cell infiltration in TIME. The gene set of labeled infiltrating immune cell type in TIME was obtained from Charoentong’s study. The gene set is rich in a variety of human immune cell subtypes.^14,15 Enrichment score calculated by ssGSEA analysis was applied to represent the relative abundance of infiltrating cells in each sample of TIME.

Gene Set Variation Analysis (GSVA)

GSVA is an R-packet, which is used to estimate the activity of pathways and biological processes in expression dataset samples, and is used for functional analysis of gene sets. Limma package was used to analyze pathways of differentially expressed genes. False discovery rate (FDR) <0.05 or log2 fold change (FC) >0.1 is considered as statistically significant.

Construction and Verification of Prognostic Nomogram

All independent prognostic parameters and related clinical parameters were tested by collinearity, and then were constructed by stepwise Cox regression model. The 1-year, 3-year, and 5-year OS rates of patients in the experimental group were predicted. Nomogram calibration curves were drawn to compare prognostic and observed OS. Besides, the differences in risk score in different clinical subgroups were compared.

Tumour Mutation Burden (TMB) Analysis

The maftools package in R was utilized to calculate the TMB score. The experimental group was divided into high and low TMB groups according to the median of TMB. The Kaplan Meier survival curves were drawn for analysis. To evaluate the distribution of somatic variation in driving genes, the top 20 driving genes with the highest change frequency between the two subgroups were obtained.

Immune Checkpoints Analysis and Gene Ontology (GO) Functional Analysis

Expression of three immune checkpoints PD-L1, PD-1, and CTLA4 in high- and low-risk score groups was analyzed. The Limma package was used for gene differential expression analysis, with thresholds set at FDR <0.01 and log2 (FC) >0.6. Differentially expressed genes between high- and low-risk score groups were enriched and analyzed by using the David database. The ggord function in R was used for principal components analysis (PCA) analysis and to determine the effectiveness of grouping.

Prediction of the Tumor Risk Score in Immunotherapy Benefits

To explore the prediction of the tumor risk score in the benefit of immunotherapy, we analyzed the response of patients in the high-risk score and low-risk score on the immunotherapy cohort groups. The difference of risk score between patients with response to immunotherapy and non-response was analyzed.

Results

Establishment of Pyroptosis-Related Risk-Score Model

A total of 33 pyroptosis-related genes were included in univariate Cox regression analysis. Six genes were identified with the P-value of 0.05 (Figure 1A). Then, LASSO Cox regression analysis showed that 5 genes could be used to construct prognostic risk models in the optimum λ value (Figure 1B and C). The formula for calculating the risk score is as follows: Risk score = (0.23539672 * CASP9) + (−0.251204085 * NLRC4) + (−0.204675619 * NLRP1) + (−0.198883793* NOD1) + (−0.129123952 * PLCG1).

Figure 1 The risk-score model in LUAD patients. (A) Hazard ratio of univariate Cox analysis for pyroptosis-related genes; (B) Coefficients for 5 genes analyzed by LASSO; (C) Distribution of LASSO coefficients for 5 genes. (D and E) Survival curves of experimental group and verification group (F and G) ROV curve of experimental group and verification group.

Based on the median risk score, the patients were divided into low-risk and high-risk groups. As shown in Figure 1D, the overall survival (OS) was worse in the high-risk group. In the validation set, the risk model based on pyroptosis-related genes can effectively predict the prognosis of patients (Figure 1E). The risk-score model could serve well to predict the OS both in experimental and validation group (Figure 1F and G). In the experimental group, the AUC of ROC curves in 1-year, 3-year, and 5-year were 0.657, 0.596, and 0.619, while in the validation group the AUC of ROC curves in 1-year, 3-year, and 5-year were 0.613, 0.557, and 0.661.

Expressions of Characteristics Genes

The expression analysis showed that the expression of NLRC4, NLRP1, and NOD1 was significantly reduced compared to normal samples, while the expression of PLCG1 was significantly increased (p-value <0.05), and there was no significant difference in the expression of CASP9 (Supplementary Figure 1A). qRT-PCR results showed that the expression of NLRC4 in tumor tissues was significantly decreased compared with normal tissue samples, which was consistent with the gene expression in the experimental group. Although there were no significant differences in NLRP1, NOD1 and PLCG1 among the tissue samples, the gene expression trend was consistent with that in the experimental group (Figure 2A). The mutation probability of characteristic genes was not high (Figure 2B). Only 6.52% of tumor samples have gene mutations, and the mutation proportions of the 5 characteristic genes were 2%, 2%, 1%, 1%, and 1%, respectively. A total of 140 transcription factors (TFs) were screened with the Pearson correlation coefficient > 0.3 and p-value <0.05. Most TFs were involved in the expression of PLCG1 (n = 66) and NLRP1 (n = 58) in the 5 genes (Figure 2C).

Figure 2 Expressions of characteristics genes and degree of immune cell infiltration in LUAD patients. (A) Expression of characteristic genes in tissue samples; (B) Mutant waterfall map of characteristic genes; (C) Transcription factors of characteristic genes; (D) Difference in infiltration degree of 23 immune cells in high and low-risk score patients. *P < 0.05; **P < 0.01; *** P < 0.001; **** P < 0.0001. Abbreviations: ns, not significant.

Immune Infiltration and GSVA Analysis

The ssGSEA was utilized to evaluate the state of 23 immune cell infiltrations in the experimental group. The infiltration of immune cells in the low-risk score group was significantly higher than that in the high-risk group (Figure 2D). The immune score, stromal score, tumor purity, and estimate score of LUAD patients in the low-risk group were significantly higher than those in the high-risk group, while the tumor purity was lower than that in the high-risk group (Supplementary Figure 1B). The differences of metabolic pathways between high-risk and low-risk groups were analyzed by GSVA and calculated with R package Limma. A total of 141 pathways with FDR <0.05 and log2FC >0.1 were obtained. The changes in metabolic pathways may have a significant impact on the occurrence and development of LUAD. The top 20 results are listed in Supplementary Figure 2A.

Construction and Validation of the Nomogram Based on Risk Score and Clinical Characteristics

Multivariate Cox analysis showed that risk score (HR: 3.01, 95% CI: 1.81–5.02) was an independent prognostic factor different from age, stage, M stage, N stage, T stage, and gender, while N stage and T stage also had independent prognostic significance (Figure 3A). The nomogram was created by combining independent prognostic indicators of risk score, T stage, and N stages in multivariate Cox analysis (Figure 3B), and the probability of 1-year, 3-year, and 5-year OS was predicted. The calibration curve showed that the nomogram demonstrated high accuracy in 1-year OS (Figure 3C), indicating that the nomogram constructed using the combined model predicts the short-term survival of patients with LUAD is better than using univariable analysis. The risk score in different clinical subgroups is demonstrated in Figure 3D. Additionally, Kaplan Meier analysis of the impact of risk scores on survival in different clinical groups (stage 1–2, stage 3–4, stage 1–2, stage 3–4) showed that in all clinical groups, the OS of patients in the high-risk score group was significantly lower than that in the low-risk score group (P < 0.05). Notably, the risk score is more effective in predicting of survival in LUAD patients with stage T 3–4, the AUC of ROC curves in 1-year, 3-year, and 5-year were 0.75, 0.772, and 0.613 (Figure 4).

Figure 3 Relationship between risk score and clinical characteristics. (A) Multivariate Cox analysis of clinical features and risk scores; (B) The nomogram for predicting the probability of OS; (C) calibration diagram of nomograph in 1, 3, and 5-year in patients with different TMB; (D) risk scores among different clinical subgroups. *P < 0.05; **P < 0.01; *** P < 0.001. Abbreviations: ns, not significant.

Figure 4 The impact of risk score on survival in different clinical groups. (A) The impact of risk score on survival in Stage 1–2; (B) The impact of risk score on survival in Stage 3–4; (C) The impact of risk score on survival in T 1–2; (D) The impact of risk score on survival in T 3–4.

Difference of TMB Between High and Low-Risk Score Groups

The tumor samples in the experimental group were divided into high and low TMB score groups according to the median TMB. The TMB of patients in the subgroup with the high-risk score was significantly higher than those in low-risk score (Figure 5A and B). The survival of patients with low TMB score was significantly higher than that of patients with high TMB score (Figure 5C). To further evaluate the distribution of somatic variation in driver genes between low- and high-risk score subgroups, the top 20 driver genes with the highest change frequency were compared (Figure 5D and E). The analysis results showed significant differences in the mutation spectrum between high- and low-risk score subgroups, which may provide a new insight for mechanism research of pyroptosis and gene mutation in immune checkpoints.

Figure 5 Characters of patients with different risk score and TMB. A and (B) Bar chart and violin plot of TMB distribution; (C) survival curve of patients with high and low TMB; (D and E) Waterfall map of gene mutation in patients with high-risk score group and low-risk score group.

Correlation Between Risk Score, Characteristic Genes with TMB, and Immune Checkpoint

In the current study, the expression of three immune checkpoints, PD-L1, PD-1, and CTLA4 in the high- and low-risk score groups was evaluated. The risk score was positively correlated with the TMB and negatively associated with the three immune checkpoints (Figure 6A). The characteristic gene CASP9 is not correlated with TMB, but is negatively associated with three immune checkpoints (Figure 6B). The characteristic genes NLRP1 (Figure 6C) and NLRC4 (Figure 6D) were negatively correlated with the TMB and positively correlated with the three immune checkpoints. Moreover, the characteristic gene NOD1 was negatively correlated with the TMB and positively correlated with the immune checkpoint CTLA4 but did not correlate with PD-L1 and PD-1 (Figure 7A). PLCG1 was positively correlated with immune checkpoint CTLA4 and PD-1 and did not correlate with TMB and immune checkpoint PD-L1 (Figure 7B). Three immune checkpoints were highly expressed in the low-risk score group (Figure 7C). Comparing risk score, characteristic genes with the results of TMB and immune checkpoints, four combinations were found to have high correlations, which is characteristic genes NLRC4 and CTLA4 (0.4), NLRC4 and PD-L1 (0.44), NLRP1, and CTLA4 (0.42), and risk score and CTLA4 (−0.42). Therefore, it can be concluded that the correlation between risk score, characteristic genes, and immune checkpoint CTLA4 is stronger.

Figure 6 Correlation between characteristic genes and immune checkpoints. (A) TMB, immune checkpoints with risk score; (B) TMB, immune checkpoints with CASP9; (C) TMB, immune checkpoints with NLRP1; (D) TMB, immune checkpoints with NLRC4. Y-axis from left to right: TMB, PD-1, PD-L1, CTLA4.

Figure 7 Correlation between characteristic genes, immune checkpoints, and risk score. (A) TMB, immune checkpoints with NOD1, Y-axis from left to right: TMB, PD-1, PD-L1, CTLA4; (B) TMB, immune checkpoints with PLCG1, Y-axis from left to right: TMB, PD-1, PD-L1, CTLA4; (C) expression of immune checkpoints in high and low-risk score groups. Y-axis from left to right: PD-1, PD-L1, CTLA4.

Differential Expression Analysis in High and Low-Risk Groups and GO Enrichment of Differential Genes

The genes with differential expression were analyzed by Limma package of R. Totally, 432 differentially expressed genes were identified, 88 genes were up-regulated, and 344 genes were down-regulated (Supplementary Figure 2B). Enrichment analysis of differentially expressed genes between high- and low-risk score groups was performed using the David database. These differentially expressed genes were mainly enriched in protein binding, plasma membrane, extracellular exosome, immune response, etc. (Figure 8A). Meanwhile, PCA analysis was performed by using the function of R to determine the effectiveness of grouping (Figure 8B).

Figure 8 GO enrichment analysis immunotherapy responses in high and low-risk score groups. (A) GO enrichment analysis results; (B) PCA plot; (C) the proportion of patients responding to immunotherapy; (D) risk score between different clinical response groups.

Prognostic of the Tumor Risk Score in the Benefit of Immunotherapy

The relevant analysis was performed based on the immunotherapy cohort GSE61676 data set in the GEO database. Patients receiving immunotherapy in the dataset were divided into response and non-response groups. The samples were divided into high-risk and low-risk groups according to the calculated median risk score. The analysis results showed the objective response rate of immunotherapy in the low-risk score group was significantly higher than that in the high-risk group (Figure 8C). Furthermore, the risk score of patients who responded to immunotherapy was significantly lower than that in the non-response group (Figure 8D). Overall, these data suggest that the risk score constructed by pyroptosis-related genes may be related to the response to immunotherapy.

Discussion

Although minimally invasive surgery, chemotherapy, and targeted therapy have made some progress in the treatment of NSCLC, the majority of cases diagnosed as advanced still make the prognosis of NSCLC not optimistic, with a 5-year survival rate of NSCLC patients is less than 25%.¹⁶ Pyroptosis is a newly discovered programmed cell death mode different from apoptosis and necrosis, which depends on the activation of caspase-1 and caspase-4/5/11, accompanied by inflammation. In addition to being closely associated with infection and other diseases, pyroptosis is also closely related to the occurrence and progress of tumors, including NSCLC. The current study included data from three cohorts of LUAD patients in UCSC Xena, TCGA, and GSE databases, constructed the risk score model related to pyroptosis-related genes to evaluate gene expression, mutation data, and immunotherapy response in LUAD patients and explored the prognostic value of these genes and their correlation with tumor immune microenvironment. Using LASSO Cox regression analysis, a prognostic gene model was constructed based on 5 pyroptosis-related genes signature: NLRC4, NLRP1, NOD1, PLCG1and CASP9, and clarified the expression and prognostic value of these genes in LUAD.

Prognosis analysis suggested a poor survival rate in LUAD patients with low expression of NLRP1 and PLCG1, which is consistent with previous reports.^17,18 The expression analysis of transcription factors related to characteristic genes also showed that most transcription factors were involved in the expression of PLCG1 and NLRP1. Our findings also demonstrated that the risk score constructed according to pyroptosis-related genes has a good ability to predict the OS of the experimental group. The analysis of the validation group also indicated that the risk score could predict the OS, which could provide more choices for predictive survival in LUAD and has a certain value in clinical practice.

Tumor tissue is composed of various types of cells, which constitute the TIME. The immune cells in the immune system penetrate TIME and help to regulate tumor progression. The immune infiltrates in TIME are the key to improve the response rate and developing new immunotherapy strategies. Our study showed that the pyroptosis-related risk scores were significantly associated with immune infiltration. We analyzed the infiltration states of 23 immune cells and found that the infiltration of most immune cells in the low-risk score group was significantly higher than that in the high-risk group. The further estimate algorithm score also showed that the immune score, matrix score and estimated score of the low-risk group were significantly higher than those of the high-risk group, while the tumor purity was lower than that of the high-risk group. Shen et al indicated that in LUAD with low NLRP1 expression, low immune cell infiltration was associated with a poorer prognosis.¹⁸ The current study confirms that pyroptosis plays an important role in the tumor immune microenvironment.

Immunotherapy has been introduced into the clinical treatment of NSCLC with the deepening of the research on the mechanism of tumor and immune resistance. However, the complexity of human immune pathway determines the complexity of immunotherapy effect. GO enrichment analysis showed that the differentially expressed genes were mainly involved in the regulation of the protein binding, plasma membrane, extracellular exosome, and immune response. Under physiological conditions, PD-1 combined with PD-L1 can down-regulate T cell activity and induce T cell apoptosis and then inhibit the autoimmune response. Nevertheless, mutant tumor cells may inactivate the tumor-killing effect mediated by T cells and deplete T cells by overexpression of PD-L1, thus making tumor cells to escape the cellular immune killing of the body and promoting tumor growth.¹⁹ CTLA4 is the earliest discovered immune checkpoint receptor, which is expressed on the surface of CD8⁺ effector T cells and inhibits the initiation and activation of effector T cells through competitive binding ligands and receptors. The expressions of these three most common immune checkpoints were included in our study for analysis. The current study found that the risk score was negatively associated with the three immune checkpoints, and they were highly expressed in the low-risk score group. Considering the lack of standardization of PD-1 detection and the inducibility of expression,^20,21 the results of risk score, a characteristic gene with TMB, and expressions of immune checkpoints were horizontally compared, and concluded that the correlation between risk score, characteristic gene, and immune checkpoint CTLA4 was stronger.

Cancer is characterized by genomic instability and structural changes that occur in tumor progression. This genomic variation can generate tumor antigen, which may be recognized by the immune system as non-autoantigen and trigger a cellular immune response. The emergence and clinical application of immune checkpoint inhibitors brought new hope for immunotherapy. The effective treatment response prognosis can promote a more efficient and safe clinical application of immune checkpoint inhibitors. NSCLC is a tumor with a high level of TMB.²² Clinical studies have shown that TMB has a good correlation with the efficacy of PD-1/PD-L1 inhibitors.²³ We found that the risk score was positively associated with the TMB in LUAD. The analysis of treatment benefit showed that the objective response rate of immunotherapy in the low-risk score group was significantly higher than that in the high-risk group. These results show that the established risk score model is correct and effective in evaluating the response to immunosuppressive therapy.

Our study has some limitations. Firstly, three analysis cohorts included are from UCSC Xena, TCGA, and GSE databases, which impair the consistency of analysis results. Secondly, due to a relatively small sample size of patients for q-PCR validation, the difference was not very significant, but the trend was consistent with the experimental group. In the future, we will enhance the reliability of our research by expanding the sample size and continuously evaluating gene expression, clinical features, and treatment responses of the same group of patients.

Conclusion

In conclusion, a pyroptosis-related risk-score model containing five genes (NLRC4, NLRP1, NOD1, PLCG1, and CASP9) for LUAD patients was constructed, and some of the genes in the risk model were validated in 5 patients with LUAD. The risk score model has good prognostic ability in survival evaluating and response to immunotherapy. Our results showed pyroptosis also played important roles in the immunotherapeutic response of LUAD.