Identification of a pyroptosis-related lncRNA risk model for predicting prognosis and immune response in colon adenocarcinoma

Background Colon adenocarcinoma (COAD) is one of the most common malignant tumors and is diagnosed at an advanced stage with a poor prognosis worldwide. Pyroptosis is involved in the initiation and progression of tumors. This research focused on constructing a pyroptosis-related ceRNA network to generate a reliable risk model for risk prediction and immune infiltration analysis of COAD. Methods Transcriptome data, miRNA-sequencing data, and clinical information were downloaded from the TCGA database. First, differentially expressed mRNAs (DEmRNAs), miRNAs (DEmiRNAs), and lncRNAs (DElncRNAs) were identified to construct a pyroptosis-related ceRNA network. Second, a pyroptosis-related lncRNA risk model was developed applying univariate Cox regression analysis and least absolute shrinkage and selection operator method (LASSO) regression analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses were utilized to functionally annotate RNAs contained in the ceRNA network. In addition, Kaplan-Meier analysis, receiver operating characteristic (ROC) curves, univariate and multivariate Cox regression, and nomogram were applied to validate this risk model. Finally, the relationship of this risk model with immune cells and immune checkpoint blockade (ICB)-related genes was analyzed. Results A total of 5373 DEmRNAs, 1159 DElncRNAs, and 355 DEmiRNAs were identified. A pyroptosis-related ceRNA regulatory network containing 132 lncRNAs, 7 miRNAs, and 5 mRNAs was constructed, and a ceRNA-based pyroptosis-related risk model including 11 lncRNAs was built. The tumor tissues were classified into high- and low-risk groups according to the median risk score. Kaplan-Meier analysis showed that the high-risk group had a shorter survival time; ROC analysis, independent prognostic analysis, and nomogram further indicated the risk model was a significant independent prognostic factor what had an excellent ability to predict patients’ risk. Moreover, immune infiltration analysis indicated that the risk model was related to immune infiltration cells (i.e., B cell naïve, T cell follicular helper, macrophage M1) and ICB-related genes (i.e., PD-1, CTLA4, HAVCR2). Conclusions This pyroptosis-related lncRNA risk model possessed good prognostic value, and the ability to predict the outcome of ICB immunotherapy in COAD. Supplementary Information The online version contains supplementary material available at 10.1186/s12957-022-02572-8.


Introduction
Colorectal cancer (CRC) is the most common cancer diagnosed in the world [1]. The incidence and mortality rates of CRC are in the top three of all cancers based on the American Cancer Society 2021 report [1]. Colon adenocarcinoma (COAD) is the most common histological subtype of CRC. With the advancements in the diagnosis and treatment of COAD in recent years, its incidence and mortality remain at 10.2% and 9.2%, respectively [2]. Therefore, the improvement of early diagnosis and treatment modalities for COAD patients is an urgent clinical need.
Pyroptosis is a gasdermin-mediated inflammatory programmed cell death characterized by cell swelling, pore formation, and the release of intracellular contents, such as IL-1β and IL-18 [3]. Pyroptosis is typically triggered by canonical pathways and noncanonical pathways [4,5]. In the past several years, an increasing number of studies have depicted that pyroptosis is involved in the progression of cancer. The primary therapeutic strategy of cancer is to induce cell death, and some researchers are trying to find novel targeted therapies for COAD by activating pyroptosis pathways [6].
The competitive endogenous RNA (ceRNA) hypothesis, including non-coding RNAs and mRNAs, is considered as a novel regulatory network, which reveals a novel mechanism of interaction between RNAs. These ceRNA molecules can compete to bind the same miRNA through microRNA response elements (MRE) to affect the gene expression [7]. Long non-coding RNA (lncRNA), more than 200 nucleotides in length, is defined as a non-coding RNA and has been found to be involved in diverse key biological processes, including cell proliferation and differentiation, genetic regulation of gene expression, and regulation of micro-RNAs (miRNAs) [8]. Studies have shown that lncRNAs can disrupt the balance of the ceRNA network, thereby promoting cancer progression [9,10]. For example, Ma et al. depicted that the lncRNA RP1-85F18.6 was upregulated in CRC and played major roles in tumorigenesis and repressed the pyroptosis of CRC cells [11]. In addition, several studies have described that lncR-NAs promote tumorigenesis by changing the immune microenvironment in cancers [12,13]. To date, the pyroptosis-related ceRNA networks have not been elucidated in COAD.
In this study, transcriptome and miRNA sequencing data between COAD tumor tissues and normal tissues were retrieved from The Cancer Genome Atlas (TCGA) database (https:// portal. gdc. cancer. gov/). The pyroptosis-related ceRNA network was constructed using integrated analysis. A pyroptosis-related prognostic signature was extracted from the ceRNA network. Then, we investigated the role of this pyroptosis-related lncRNA prognostic signature in immune microenvironment and immune checkpoint inhibitor treatment.

Data acquisition
RNA-seq data of COAD was derived from the TCGA-COAD dataset including 398 tumor samples and 39 normal samples. The miRNA data was downloaded from the TCGA including 380 tumor samples and 8 normal samples. Meanwhile, relevant clinical data of 385 COAD patients was also obtained from TCGA, including age, gender, survival status, stage, T, N, and M classification ( Table 1).

Identification and validation of a ceRNA-based pyroptosis-related lncRNA risk model
A total of 379 COAD patients with survival data were included. One hundred thirty-two pyroptosis-related lncRNAs (PRlncRNAs) in the above ceRNA network were analyzed by univariate Cox regression analysis to filter PRlncRNAs associated with survival. To avoid overfitting, PRlncRNAs were screened via least absolute shrinkage and selection operator (LASSO) regression analysis (R package "glmnet", P<0.05) [19,20]. Then, a PRlncRNA risk model was constructed and the risk score of each sample was calculated based on the formula below: risk score = ∑ i Xi × Yi (X: coefficients, Y: lncRNA expression level) [21]. And tumor tissues were separated into high-and low-risk groups according to the median risk score. Next, lncRNAs in the PRlncRNA risk model with associated DEmiRNAs and DEmRNAs were used to construct a prognostic ceRNA network. To assess the putative biological role of RNAs in the ceRNA network, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were conducted (R packages "clusterProfiler", "org.Hs.eg. db, " and "enrichplot") [22].

Validation of the risk model
Kaplan-Meier analysis was used to compare the overall survival (OS) time between the two risk groups using R packages "survival" and "survminer" [23]. The 1-, 3-, and 5-year receiver operating characteristic (ROC) curves were acquired by utilizing the "timeROC", "survival, " and "survminer" R packages [24]. The area under the receiver operating curve (AUC) was used to indicate forecast performance. Independent prognostic analysis showed the relationship of risk score and clinical traits through the "survival", "survminer, " and "forestplot" R packages [23].

Establishment of the nomogram
A nomogram was established to quantitatively calculate patient survival. Then, calibration curves were employed to confirm the predictive effect of nomogram (R package "rms").

Immune cell infiltration and immune checkpoint blockade (ICB) analysis
To further investigate the difference between immune cell infiltration in the two risk groups, the relationship between the immune cell proportion and risk score was estimated by the CIBERSORT algorithm (P<0.05). The Spearman correlation analysis was employed to estimate the association of the immune cells and risk score. The expression level of ICB-related genes was closely involved with the outcome of immunotherapy. Therefore, the Spearman correlation between the 11 ICBrelated genes and risk score was analyzed [25].

Statistical analysis
All statistical analyses were performed using Perl v5.32.1 and R 4.1.0 software. The ceRNA network was constructed by Cytoscape v3.8.2. The Wilcoxon test was used to compare the proportion of immune cells between the two groups. DEmRNAs, DElncRNAs, and DEmiRNAs were gained with the thresholds FDR <0.05 and |log2 fold change|>1. Statistical tests were twotailed (P<0.05).

Construction of a pyroptosis-related ceRNA network in COAD
To construct the pyroptosis-related ceRNA network that facilitates understanding of the links between DEmRNAs, DElncRNAs, and DEmiRNAs, several databases were used.

Functional enrichment analysis of pyroptosis-related ceRNA genes
To explore potential functions of these 11 lncRNAs, 7 miRNAs, and 5 mRNAs in biological processes, GO and KEGG functional enrichment analyses were performed with P value < 0.05 as the threshold. For GO analysis, these genes were primarily enriched in the cellular response to biotic stimulus and lipid catabolic process in biological processes (BPs); TORC2 complex, GATOR2 complex, and TOR complex in cellular components (CCs); and ubiquitin-like protein ligase binding in molecular function (MF) (Fig. 4E). KEGG pathway analysis suggested that these genes were markedly concentrated in the IL-17 signaling pathway, TNF signaling pathway, tuberculosis, and NOD-like receptor signaling pathway (Fig. 4F).

Validation of the prognostic PRlncRNA risk model
In the prognostic PRlncRNA risk model,   prognostic capability of the PRlncRNA risk model in predicting COAD patient overall survival, COAD patients were categorized into high-or low-risk groups in terms of the median risk value. Kaplan-Meier analysis showed that patients in the high-risk group had shorter survival times (Fig. 5B). Likewise, patients were separated into high-risk and low-risk groups based on the median risk score (Fig. 5C). As the risk score increased, the patient's survival time decreased gradually (Fig. 5D). ROC analysis was used to evaluate the predictive power of this prognostic model, which indicated that the PRlncRNA risk model was able to excellently predict the 1-year (0.744), 3-year (0.696), and 5-year (0.623) survival of COAD patients, respectively (Fig. 5E).

Independent prognostic analysis of the prognostic PRlncRNA risk model
To explore whether the prognostic PRlncRNA risk model can independently predict the prognosis of COAD, univariate and multivariate Cox regression analyses were employed. The univariate Cox regression analysis showed that age (P = 0.018), clinical stage (P<0.01), T stage (P<0.001), risk score (P<0.001), N stage (P<0.001), and M stage (P<0.028) predicted dismal OS (Fig. 6A). And the results of multivariate Cox regression analysis confirmed the independence of the prognostic PRlncRNA risk model for predicting COAD prognosis (Fig. 6B). Moreover, the heatmap of clinical characteristics implied that the survival status of patients was differentially distributed between the lowand high-risk subgroups (P<0.05, Fig. 6C).

Construction of a predictive nomogram
To further assess whether this prognostic PRlncRNA risk model had optimal predictive capabilities, we collected clinical characteristics, including age, gender, and stage, as candidate predictive biomolecular indicators to construct a nomogram (Fig. 7A). The calibration curves of 1 year, 3 years, and 5 years suggested that the nomogram was close to the actual morality and had good predictive power (Fig. 7B-D). The above findings showed a promising capacity of the PRlncRNA risk model for patient prognosis and survival prediction.

Immunoinfiltration analysis
Studies have shown that the pyroptosis of tumor cells can effectively regulate the tumor immune microenvironment (TIME) and activate a strong T cell-mediated anti-tumor immune response [26]. To further explore the relationship between this PRlncRNA risk model and TIME, the CIBERSOR algorithm was utilized to compare the proportions of 22 immune infiltrating cells between the high-risk and low-risk groups. B cell naïve, T cell CD8, T cell follicular helper, T cell regulatory (Tregs), and macrophage M1 were higher expressed while T cell CD4 memory resting, dendritic cells activated, and mast cells activated were lower expressed in the highrisk group by using Wilcoxon rank-sum test (Fig. 8A, B). These results provided strong evidence that our prognostic risk model was significantly associated with immune cells infiltration in COAD. Additionally, we performed correlation analysis between the prognostic risk model with the above eight significantly different immune infiltrating cells, and the results showed that the risk score   (Fig. 8C), indicating that this risk model might be a key part in the assessment of responsiveness to ICB immunotherapy in COAD.

Relationship of the prognostic PRlncRNA risk model and immune checkpoint blockades (ICBs)
ICB-related immunotherapy has become a promising modality for patients with COAD [27]. To investigate the response of COAD samples to immunotherapy, we examined the association of ICB-related genes (i.e., PD-1, CTLA4, HAVCR2) and risk score (Fig. 9A). The results showed that the risk score was significantly positively .042), indicating that the risk score model was useful in assessing patients' response to immunotherapy (Fig. 9B).

Discussion
Colon adenocarcinoma (COAD) is the most common malignancy with increased mortality worldwide [28]. Despite significant improvements in surgery, radiotherapy, chemotherapy, and immunotherapy, the rates of 5-year survival remain low [29]. Therefore, it is crucial to identify potential biomarkers for the diagnosis and treatment of COAD. Pyroptosis was regarded as a new type of programmed cell death that played a dual function in the development of cancer [3]. In recent years, pyroptosis has become a hot topic in the field of oncology research, and an increasing number of researches have emphasized the key effects of pyroptosis in tumorigenesis and TIME [26,30,31]. ceRNA represents a new regulation mode of gene expression, which is more sophisticated and complex than the miRNA regulatory network. There is growing evidence that the ceRNA network is involved in the progression of COAD [32,33]. Nevertheless, the underlying function of the ceRNAbased pyroptosis-related risk model in prognostic prediction and tumor immunity of COAD has not yet been elucidated, and our study was designed to clarify this role.
Our study constructed a pyroptosis-related ceRNA network and PRlncRNA risk model. Moreover, we also explored the prognostic predictive ability of the risk model and its association with immune cell infiltration and assessed the reactivity of COAD patients to ICB therapy. In current research, a pyroptosis-related ceRNA regulatory network including 5 mRNAs, 7 miRNAs, and 11 lncRNAs was first constructed to investigate the potential molecular mechanism of the ceRNA network. Furthermore, GO and KEGG enrichment analyses showed that these genes were mostly enriched in the IL-17 signaling pathway and TNF signaling pathway. These results suggested that this pyroptosis-related ceRNA network can be a new tool to predict clinical results of COAD. However, these discoveries need to be verified in additional studies. Moreover, 11 pyroptosis-related lncRNAs were incorporated into developing a PRlncRNA risk model. Then, Kaplan-Meier curve, time-dependent ROC curves, Cox regression analysis, and nomogram showed that this risk model possessed excellent prediction ability and became an independent predictor of COAD prognosis.
A growing body of evidence has demonstrated that lncRNAs play an essential role in regulating immune cell infiltration [34,35], and some researches have also reported that pyroptosis of tumor cells regulates tumor-suppressed immune cells [26,31]. Furthermore, we investigated the fundamental effects of the risk score in the regulation of the TIME. Consistent with the previous reports, our research also showed that the risk score was negatively correlated with resting immune cell proportions, but positively associated with immunosuppressive cells, indicating that patients with low-risk scores were immunologically resting, while those with high-risk scores represent an immunosuppressive tumor microenvironment.
With the development of immune checkpoint inhibitors, ICB immunotherapy has generated promising therapeutic results in COAD [36,37]. To date, immunotherapy has become the fifth pillar in the foundation of COAD therapeutics. Unfortunately, the majority of COAD patients do not respond to ICB treatment [38]. Pyroptosis can alter the immune microenvironment and remodel immune cells to enhance the efficiency of tumor immunotherapies [39]. Previous researches have demonstrated that pyroptosis induction plus PD-1 improves the anti-tumor activity [26,40]. Thus, a novel PRlncRNA risk model was developed to investigate the correlation of pyroptosis and ICB-related genes and to predict COAD patients' responses to ICB immunotherapy. The current study has shown that the PRlncRNA risk model was strongly related to ICB-related genes (i.e., PD-1, PD-L1), which implied that the PRlncRNA risk model could be utilized to evaluate the response to ICB treatment of COAD patients. Meanwhile, further validation of the PRlncRNA risk model as a useful predictor of immune checkpoint therapy in COAD is needed in the future.
Our research has some limitations. All analyses were employed by using the TCGA-COAD cohort, which would be better to validate with other database cohorts. In addition, in vivo and in vitro experiments should be conducted to further verify our results. However, the novelty of our study is that for the first time, the molecular mechanism of COAD was investigated from the perspective of the pyroptosis-related ceRNA network. Additionally, 11 pyroptosis-related lncRNA prognostic biomarkers were screened based on the ceRNA network. Moreover, the PRlncRNA risk model possessed a high predictive ability for survival in COAD patients. This may provide a new idea for the study of COAD.

Conclusions
To sum up, we performed comprehensive and systematic bioinformatics analysis and constructed a PRl-ncRNA risk model for COAD patients, which could be used as a potent tool in predicting the prognosis of COAD patients. In addition, the risk model was related to TIME-and ICB-related genes. The pyroptosis-related ceRNA network might be a promising therapeutic target in COAD. Fig. 9 The connection between the prognostic risk score and ICB-related genes. A Heatmap of the correlation of the risk score and 11 ICB-related genes. C The connection between the risk score and 11 ICB-related genes, respectively (P<0.05)