Gene transcriptome profiling data, mutation data, and corresponding clinical information of LUAD patients were downloaded from (https://portal.gdc.cancer.gov/). Fragments per kilobase of exon model per million mapped fragments (FPKM) were converted to l o g 2   ( F P K M + 1 )   as a reflection of the gene expression level to visually display the results when constructing the figures. The dataset from TCGA served as a training cohort to construct the pyroptosis-related lncRNA prognostic model.

A total of 45 LUAD patient specimens were recruited at Shandong Provincial Hospital, Shandong, China. The biomedical research ethic committee of Shandong Provincial Hospital approved this study (SWYX: NO. 2021-433).

We identified 33 pyroptosis-related genes from the literature, as listed in Supplementary Table S1 (Deets and Vance, 2021; Liu et al., 2021; Ye et al., 2021). The R software (R, vision 4.1.0) package “limma” (Ritchie et al., 2015) was used to determine the differently expressed pyroptosis-related genes with the absolute log2-fold change (   | l o g 2 FC | ) > 1 and adjusted p value <0.05. In addition, copy number variations (CNVs), mutation frequencies, and location analyses of pyroptosis-related genes were performed using “maftools” (Mayakonda et al., 2018) and “RCircos” with R software (Zhang et al., 2013).

Pyroptosis-related lncRNAs were extracted using Pearson’s correlation analysis between pyroptosis-related genes and lncRNAs (correlation coefficient >0.60, p < 0.001, and false-discovery rate (FDR) < 0.05).

We extracted differently expressed pyroptosis-related lncRNAs with an | l o g 2 FC | > 1 and an adjusted p value <0.05 using R package “limma”. Subsequently, 17 prognostic pyroptosis-related lncRNAs were identified by applying univariate Cox regression analysis (p < 0.05). A hazard ratio (HR) < 1 represented that lncRNAs were protective factors. To explore the regulatory mechanism of the selected pyroptosis-related lncRNAs, an mRNA–lncRNA regulatory network was constructed based on the 17 lncRNAs and two co-expressed mRNAs.

The Kaplan–Meier plotter (https://kmplot.com/analysis/) was used to analyze the impact of co-expressed mRNAs, while TIMER2.0 (Li et al., 2020) (http://timer.cistrome.org/) and gene expression profiling integrative analysis (Tang et al., 2017) (GEPIA, http://gepia.cancer-pku.cn/) were used for differential expression analysis of co-expressed mRNAs.

According to the prognostic-related lncRNAs, multivariate Cox regression analysis was performed to identify the best prognostic signature, which consisted of three lncRNAs, namely, AC090559.1, AC026355.2, and AC034102.8. We calculated the risk score as follows: risk score = (−0.2603×AC090559.1) + (−0.0974×AC026355.2) + (−0.9235×AC034102.8). Based on the median of risk score, the LUAD patients from the TCGA database were divided into high- and low-risk groups for further analysis. Subsequently, risk score distribution maps, survival status maps, and lncRNA expression heat maps were plotted. Kaplan–Meier analysis was applied to compare the overall survival (OS) of the two groups. The time receiver operating characteristic (ROC) curve was plotted using the “timeROC” package with R software, which was used to evaluate the predictive capability of the risk model. The area under the curve (AUC) of the constructed risk model in predicting the OS of LUAD patients was compared with several previously published lncRNA signatures, including the ferroptosis-related lncRNA signature of Lu (Lu et al., 2021), hypoxia-associated lncRNA signature of Shao (Shao et al., 2021), autophagy-related lncRNA signature of Liu (Liu and Yang, 2021), and six-lncRNA-based prognostic signatures of Yang (Yang et al., 2021).

The prognostic significance of the risk score and other clinical characteristics was evaluated by using univariate and multivariate Cox regression analyses. A nomogram was established to predict the 1-, 3-, and 5-year OS, which consisted of the risk score, age, and stage. Calibration curves were plotted to assess the accuracy of the risk model.

GO enrichment analysis was applied to investigate the potential pathways of the differentially expressed pyroptosis genes using “org.Hs.eg.db,” “clusterProfiler,” and “enrichplot” modules within the R package. By applying the “GSVA” tool, GSVA enrichment analysis was performed based on hallmark gene sets extracted from the MSigDB database (Hänzelmann et al., 2013).

Single-sample gene set enrichment analysis (ssGSEA) was performed to compare the difference in abundance of 28 types of infiltrating immune cells, 13 immune functions, and 13 other tumor-related biological processes, which were extracted from previous articles (Şenbabaoğlu et al., 2016; Charoentong et al., 2017; Mariathasan et al., 2018). In addition, we used computational methods to assess the infiltrating immune cells, including the TIMER (Li et al., 2020), CIBERSORT (Chen et al., 2018), quanTIseq (Finotello et al., 2019), MCP-counter (Becht et al., 2016), xCell (Aran et al., 2017), EPIC (Racle et al., 2017), and TIDE (Jiang et al., 2018) algorithms. Immune score and tumor purity were calculated using the “ESTIMATE” tool within the R package (Chakraborty and Hossain, 2018). In addition, correlation analysis was performed using the correlation heat map tool in HiPlot (https://hiplot.com.cn), a comprehensive web platform for scientific data visualization.

The immunophenoscore (IPS) was obtained from The Cancer Immunome Atlas (https://tcia.at/). Information on the dysfunction and exclusion of infiltrating cytotoxic T lymphocytes (CTLs) was downloaded from the Tumor Immune Dysfunction and Exclusion (TIDE) website (http://tide.dfci.harvard.edu/). In addition, TIDE was used to evaluate patients who received a benefit or no benefit from ICB therapy through the comprehensive biomarkers of the ICB response in different groups (Fu et al., 2020).

Following the manufacturer’s protocol, total RNA was extracted from clinical specimens using AG RNAex Pro Reagent (Accurate Biotechnology (Hunan) Co., Ltd., China). The Evo M-MLVRT Master Mix kit (Accurate Biotechnology (Hunan) Co., Ltd., China) was used for reverse transcription to obtain cDNA. Relative gene expression was detected using the SYBR Premix Ex Tap kit (Accurate Biotechnology (Hunan) Co., Ltd., China) and normalized to the expression using 18S. The primers are listed in Supplementary Table S1.

Student’s t test was used to compare the differences between the two groups. Pearson’s correlation test was used for correlation analysis. Survival analysis was performed by the Kaplan–Meier method and compared with the log-rank test. All statistical analyses were performed using R software (vision 4.1.0), and p < 0.05 was considered statistically significant.

We identified 33 pyroptosis-related genes and performed differential gene expression analysis between LUAD and normal lung tissues. A total of 15 genes, namely, AIM2, CASP4, CASP8, GSDME, NLRP7, GPX4, CASP6, TIRAP, GSDMD, PLCG1, GSDMA, GSDMC, GSDMB, CASP3, and PJVK were highly expressed in LUAD, while the expression of IL1B, IL6, NLRP3, IL18, PYCARD, TNF, CASP1, PRKACA, NLRP1, CASP5, NOD1, ELANE, and NLRC4 was decreased (Figure 1A). Subsequently, we developed a panorama of the somatic mutations of pyroptosis-related genes. A total of 30.48% of the samples had mutations, and the three genes with the highest mutation rates were NLRP3 (11%), NLRP7 (5%), and NLRP2 (4%) (Figure 1B). Figures 1C and D show the frequency of the CNV alterations of the 33 pyroptosis-related genes and their locations on the chromosome. The frequency of copy number mutations in pyroptosis-related genes was greater than the frequency of copy number deletions. To clarify the functions of the pyroptosis-related genes, we conducted GO functional enrichment analyses. Apart from the pyroptosis pathway, pyroptosis-related genes were also enriched in defense response to bacterium, regulation of interleukin-1 or interleukin-1β production, and inflammasome complex (Figure 1E).

We first performed correlation analysis (correlation coefficients >0.60, p < 0.001, FDR<0.05) between 14,057 lncRNAs and 33 pyroptosis-related genes in LUAD samples, and 1,070 pyroptosis-related lncRNAs were identified. Subsequently, 320 differently expressed pyroptosis-related lncRNAs were identified and exhibited in a volcano map after differential expression analysis (Figure 2A). Finally, 17 lncRNAs related to LUAD prognosis were identified after applying univariate COX regression analysis (Table 1). A heat map was plotted to show the differential expression of the 17 lncRNAs between normal and tumor tissues (Figure 2B). The forest map indicated that the hazard ratios of these 17 lncRNAs were all <1, suggesting that they were protective factors for prognosis (Figure 2C). In addition, we used Cytoscope to construct a co-expression network for the 17 pyroptosis-related lncRNAs and two corresponding genes (Figure 2D). The correlation scores between NLRC4, SCAF11, and 17 lncRNAs are shown in Figure 2E, suggesting that these lncRNAs may perform functions through NLRC4 and SCAF11. Next, we analyzed the impact of NLRC4 and SCAF11 on survival by applying the Kaplan–Meier plotter, which indicated that high expression was significantly associated with high overall survival (Supplementary Figure S1A). TIMER2.0 and GEPIA were used for differential gene expression analysis, and the results showed a significant downregulation of NLRC4 in LUAD but no change in SCAF11 (Supplementary Figure S1B–D).

We subsequently performed multivariate COX regression analysis on the previously obtained 17 lncRNAs and three lncRNAs were identified (Figure 3A). A risk score was established based on the multivariate regression coefficients (Table 2). According to the median risk score, LUAD patients were divided into high- and low-risk groups. Most pyroptosis-related genes were significantly upregulated in the low-risk group, suggesting a more active involvement of pyroptosis (Figure 3B). By contrast, survival analysis revealed that patients in the low-risk group showed a better prognosis, indicating that pyroptosis was associated with survival advantages (Figure 3C). Principal component analysis (PCA) indicated significant distinction in transcription profiles between the two groups (Figure 3D). The risk score, survival status, and lncRNA expression are shown in Figure 3E, revealing that mortality was significantly related to risk score. Moreover, a higher risk score was significantly related to advanced stage (Figure 3F). These results indicated that a lower risk score was associated with active pyroptosis and better clinical outcome. The areas under the curves (AUCs) of the 1-year, 3-year, and 5-year ROC curve were 0.775, 0.730, and 0.705, respectively, revealing a high accuracy in the prognosis prediction of the risk model (Figure 3G). Compared with four published lncRNA signatures, the risk model we constructed had higher accuracy in predicting 1-, 3-, and 5-year survival for TCGA–LUAD patients (Figure 3H).

After incorporating the risk scores and clinical features, univariate and multivariate COX regression analyses were performed, and the results indicated that the risk score could serve as an independent factor affecting the survival of LUAD patients, similar to stage and age (Figures 4A,B). To better predict the 1-, 3-, and 5-year OS of LUAD patients, we constructed a nomogram by incorporating risk score, age, and stage (Figure 4C), and calibration curves were constructed to assess the accuracy of nomogram (Figure 4D).

To further understand the significance of the risk score, we conducted GSVA analysis. Several immune-related pathways and events, including interferon (IFN) gamma/alpha response, IL-6–JAK-STAT3 signaling, allograft rejection, and inflammatory response, were upregulated in the low-risk group, while in the high-risk group, metabolic and cancer-promoting pathways, such as oxidative phosphorylation, glycolysis, MYC signaling, E2F signaling, and MTORC1 signaling, were activated (Figure 5A).

We investigated the differences in various immune-related functions, immune-infiltrating cells, and other tumor-related functions between high- and low-risk groups using ssGSEA. As expected, the low-risk group, with a higher level of pyroptosis, was more involved in immune-related functions such as antigen-presenting cell (APC) co-stimulation/inhibition, inflammation, cytolytic activity, human leukocyte antigen (HLA) function, T cell co-stimulation/inhibition, and type I/II IFN responses, indicative of active immune functions (Figure 5B). Consistently, the higher abundance of most infiltrating immune cells existed in the low-risk group, including activated B cells, activated CD8⁺ T cells, dendritic cells, eosinophils, MDSCs, and macrophages (Figure 5C, Supplementary Figure S2). In addition, several tumor-related pathways and events, including immune checkpoint, angiogenesis, antigen processing machinery, and CD8⁺ T effector function, were also upregulated in the low-risk group. By contrast, DNA damage repair, mismatch repair, nucleotide excision repair, and DNA replication were upregulated in the high-risk group, which may play roles in genome stability and LUAD progression (Figure 5D). In addition, HLA-related genes were significantly upregulated in the low-risk group, indicative of antigen presentation (Figure 5E). Finally, tumor purity, immune score, ESTIMATE score, and stromal score were calculated using the “estimate” tool within the R package (Figure 5F). Significantly upregulated immune and stromal scores were features of the low-risk group, while the high-risk group was characterized by higher tumor purity.

Considering the significant differences in the TME landscape, we identified several ICPs and performed differential expression analysis between low- and high-risk groups. As shown in Figure 6A, all 39 selected ICPs were upregulated in the low-risk group, indicating the potential benefit of ICI therapy. Furthermore, we used IPS obtained from TIDE as a predictor of the response to anticytotoxic T lymphocyte antigen-4 (CTLA-4) and antiprogrammed cell death protein 1 (anti-PD-1) antibodies (Charoentong et al., 2017). The results showed a higher IPS level in the low-risk group, revealing a better response to combined PD1 and CTLA4 blockade therapy or PD1 monotherapy (Figures 6B,C). Considering the upregulation of immune checkpoints and infiltration of Treg cells in the low-risk group, which could suppress the effect of CD8⁺ T cells, we further investigated the status of T cells in the TME. Consistently, the low-risk group was characterized by higher T cell dysfunction and lower exclusion, suggesting a potential advantage to ICI therapy (Figures 6D,E). To further investigate the response to ICIs, we utilized the TIDE website to predict the “responder” and “nonresponders” of PD1 and CTLA4 blockade therapy in TCGA–LUAD patients, which were constructed by integrating TIDE score, IFNG, MSI, MDSC, CAFs, and other published modules. We found that 44% of patients in the low-risk group were identified as responders to PD1 and CTLA4 blockade therapy, while only 31% in the high-risk group were classified as responders (Figure 6F). In further analysis, only 4% of patients in the low-risk group showed no benefit from ICB compared with 13% in the high-risk group (Figure 6G). In addition, patients who were not responders or showed no benefit from PD1 and CTLA4 blockade therapy had significantly higher risk scores (Figures 6H,I).

In conclusion, the low-risk group could be clarified as an immune “hot” phenotype, with a high abundance of infiltrating immune cells and better efficacy for PD1 and CTLA4 blockade therapy, while the high-risk group represented the immune “cold” phenotype which was characterized by less sensitivity to ICB.

To explore whether the risk score affected immune function through NLRC4 or SCAF11, we performed a correlation analysis between NLRC4 or SCAF11 and infiltrating immune cells. The heat map showed that lncRNA AC090559.1 and its corresponding mRNA NLRC4 were closely correlated with infiltrating immune cells, suggesting that AC090559.1 may exert its regulatory function by targeting NLRC4 (Supplementary Figure S3). The scatter plots of NLRC4, SCAF11, and various infiltrating immune cells downloaded from TIMER indicated that NLRC4 was highly correlated with a variety of immune cells (Supplementary Figure S4). Finally, we predicted the sensitivity of several common chemotherapy drugs in high- and low-risk groups. The high-risk group was more sensitive to doxorubicin, sorafenib, docetaxel, and erlotinib but less sensitive to gefitinib (Supplementary Figure S5).

We detected the relative expression levels of the three lncRNAs in LUAD specimens by using qRT-PCR. The risk score was subsequently calculated, and the patients were divided into high- and low-risk groups based on the median value. The overall survival of the low-risk group was significantly better, consistent with previous results (Figure 7A). The AUC values of 3 year, 5 year and 7 year were 0.599, 0.671, and 0.704, respectively (Figure 7B). By correlation analysis of the three identified lncRNAs and immune checkpoint genes in LUAD specimens, we found a strong correlation between risk score and CTLA4, indicative of a better response to anti-CTLA4 immunotherapy (Figures 7C,D).