Identification of pyroptosis-related lncRNA signature and AC005253.1 as a pyroptosis-related oncogene in prostate cancer

Background Pyroptosis and prostate cancer (PCa) are closely related. The role of pyroptosis-related long non-coding RNAs (lncRNAs) (PRLs) in PCa remains elusive. This study aimed to explore the relationship between PRL and PCa prognosis. Methods Gene expression and clinical signatures were obtained from The Cancer Genome Atlas and Gene Expression Omnibus databases. A PRL risk prediction model was established by survival random forest analysis and least absolute shrinkage and selection operator regression. Functional enrichment, immune status, immune checkpoints, genetic mutations, and drug susceptibility analyses related to risk scores were performed by the single-sample gene set enrichment analysis, gene set variation analysis, and copy number variation analysis. PRL expression was verified in PCa cells. Cell Counting Kit-8, 5-ethynyl-2′-deoxyuridine, wound healing, transwell, and Western blotting assay were used to detect the proliferation, migration, invasion, and pyroptosis of PCa cells, respectively. Results Prognostic features based on six PRL (AC129507.1, AC005253.1, AC127502.2, AC068580.3, LIMD1-AS1, and LINC01852) were constructed, and patients in the high-score group had a worse prognosis than those in the low-score group. This feature was determined to be independent by Cox regression analysis, and the area under the curve of the 1-, 3-, and 5-year receiver operating characteristic curves in the testing cohort was 1, 0.93, and 0.92, respectively. Moreover, the external cohort validation confirmed the robustness of the PRL risk prediction model. There was a clear distinction between the immune status of the two groups. The expression of multiple immune checkpoints was also reduced in the high-score group. Gene mutation proportion in the high-score group increased, and the sensitivity to drugs increased significantly. Six PRLs were upregulated in PCa cells. Silencing of AC005253.1 inhibited cell proliferation, migration, and invasion in DU145 and PC-3 cells. Moreover, silencing of AC005253.1 promoted pyroptosis and inflammasome AIM2 expression. Conclusions Overall, we constructed a prognostic model of PCa with six PRLs and identified their expression in PCa cells. The experimental verification showed that AC005253.1 could affect the proliferation, migration, and invasion abilities of PCa cells. Meanwhile, AC005253.1 may play an important role in PCa by affecting pyroptosis through the AIM2 inflammasome. This result requires further research for verification.


Introduction
Prostate cancer (PCa) is the second most common cancer in the male population worldwide and is one of the top five causes of cancer-related death (1). Epidemiological examinations and clinical studies have found that the incidence of PCa is still rising (2). The symptoms of patients with early-stage PCa are non-specific, so PCa is often found in the middle and late stages, and surgical treatment at this time has a poor prognosis and a low survival rate (3). With several important recent discoveries in immune mechanisms and advanced molecular diagnostic platforms, immunotherapy is emerging as a viable option for PCa, especially castration-resistant PCa, to stimulate antitumor immunity (4,5). Different patient responses to the same immunotherapy have been observed in patients with different types and stages of cancer (6). Moreover, the patient response depends on multiple factors, including intratumoral heterogeneity and prior treatment history, suggesting that the need for individualized and combined therapy is an important direction for future successful immunotherapy (7,8). Diagnostic, prognostic, and predictive biomarkers enable patient-specific management of PCa (9). Specific biomarkers to facilitate the clinical selection of immunotherapy patients include programmed death ligand 1 (PD-L1) and prostatic acid phosphatase (PAP), but these approaches are limited by tumor heterogeneity or small percentage populations (10,11). Therefore, there is an urgent need to identify new and effective biomarkers to establish a prognostic model of PCa.
Pyroptosis, distinct from apoptosis, is a type of programmed cell death induced by the inflammasome and carried out by gasdermin proteins (12). It is characterized by cell rupture and many pro-inflammatory factors being discharged (13,14). Pyroptosis affects tumor cell invasion, multiplication, and migration, affecting cancer prognosis (15). The association between pyroptosis and cancer is highly intricate because as a way of cell death, pyroptosis could inhibit cancer occurrence and development. Meanwhile, the release of inflammatory mediators and various signaling pathways in pyroptosis is associated with tumorigenesis and resistance to chemotherapy (16,17). Due to the close association between pyroptosis and cancer progress and prognosis, various prognostic biomarker studies based on pyroptosis genes have been identified and used to construct gene signatures with predictive power. For example, risk signatures based on five pyroptosis-related genes (PRGs) were biomarkers to predict the immunological condition and the outcome of lung adenocarcinoma (18). Furthermore, the expression of four PRG features strongly predicted a breast cancer patient's prognosis (19). Nonetheless, the predictive merit of pyroptosis gene signatures in the prostate has not been completely clarified.
Long non-coding RNA (lncRNA) does not have the proteincoding capacity, and its length exceeds 200 nucleotides (20). LncRNA has been shown to play key functions in a variety of biological and disease processes, including cancer (21,22). Growing evidence supports the involvement of lncRNA in PCa progression, including cell proliferation, apoptosis, metastasis, and invasion (23, 24). For example, overexpression of lncRNA PCAT14 inhibits the invasion of PCa cells and correlates with a good prognosis of PCa, which can be a diagnostic marker (25,26). However, the role of pyroptosisrelated lncRNAs (PRLs) in PCa still requires further exploration. Thus, exploring lncRNA biomarkers associated with pyroptosis in PCa is of clinical importance.
Machine learning is a branch of artificial intelligence that has been rapidly developed and applied in the field of medicine (27). Predictive models of diseases based on machine learning have been extensively mined (28). For example, Wu et al. used an ensemble of machine learning to develop a novel pyroptosis scoring system based on six lncRNAs to predict the prognosis of patients with low-grade glioma (29). In the present study, we used a machine learning approach to construct a PRL risk signature for PCa prognosis. Then, we investigated the associations of risk score models and clinical features, immune microenvironment, immune checkpoints, genetic mutations, and drug sensitivity to provide potential diagnostic and prognostic biomarkers for PCa. This study may help to understand the effect of PRL on PCa prognosis.

Data sets and pretreatments
The Cancer Genome Atlas (TCGA)-PRAD (Prostate Adenocarcinoma) dataset was downloaded from UCSC Xena (https://xenabrowser.net/). The data processing was performed to obtain FPKM data directly from TCGA and convert it into TPM value, normalized by log 2 (TPM+1). The GSE116918 (GPL25318) dataset was from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/gds/). For this validation cohort, 248 localized/locally advanced PCa patients commencing radical radiotherapy (with androgen deprivation therapy (ADT)) were included. The Affymetrix platform was utilized to generate raw data from the GSE116918 (GPL25318) dataset. The robust multi-chip averaging (RMA) algorithm was used to achieve background correction and normalization.

Pathway and immune-infiltration evaluation
The single-sample gene set enrichment analysis (ssGSEA) algorithm was used to quantify the abundance of 28 immune  (38) in PRAD and to compare immune infiltration with prognostic scores. Stromal score, Immune Score, ESTIMATE Score, and Tumor Purity were assessed with the ESTIMATE package, and the relationship between these scores and prognostic scores was compared. The GSVA package was used for GSVA analysis of Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). The correlation analysis was performed with prognostic and functional enrichment pathway scores. Correlation analysis was performed with prognostic scores and all genes, and then gene set enrichment analysis (GSEA) was performed by the clusterProfiler package.

Gene mutation and copy number variation analysis
The gene mutations of the two groups were compared, and the gene mutation patterns were checked. The somatic mutations and somatic copy number variation (CNV) profiles were gathered from the TCGA-PRAD datasets. The Genomic Identification of Significant Targets in Cancer (GISTIC) analysis was performed to evaluate the genomic features. The CNV landscape based on ITGA5 levels and the copy number gains or losses at the amplified or deleted peaks were assessed by GISTIC 2.0 analysis (https://gatk.broadinstitute.org).

Drug prediction
Information on the sensitivity of tumor cell lines to potential drugs was downloaded from Cancer Therapeutics Response Portal 2 (CTRP v2) and Profiling Relative Inhibition Simultaneously in Mixtures (PRISM). The lower the area under the curve (AUC) of the cell line, the higher the sensitivity to the potential drug. Expression of cancer cell lines was downloaded from Cancer Cell Line Encyclopedia (CCLE). Predictions were made using the R package of pRRophetic.
Total RNA was obtained using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) from RWPE1, PC-3, and DU145 cells. RNA was reverse transcribed into cDNA using an mRNA reverse transcription kit (CW2569, CWBIO, Beijing, China). LncRNA expression was detected using the SYBR method (CW2601, CWBIO, China) and quantitative reverse transcription PCR (RT-qPCR) analysis with GAPDH as an internal reference. The primer sequences are shown in Table 2. The relative expression levels of genes were investigated by 2 −DDCt .

5-Ethynyl-2′-deoxyuridine assay
The 5-ethynyl-2′-deoxyuridine (EDU) assay kit (Guangzhou RiboBio, Guangzhou, China) was used to monitor cell proliferation. The cells were inoculated into 96-well plates (1 × 10 4 /well) until 80% confluence. Each well was added with 100 ml of EDU solution and incubated for 2 h. The cells were incubated with 4% paraformaldehyde at room temperature for 30 min. Then, the cells were treated with 100 ml of 1× Apollo ® staining reaction solution for 30 min. Next, 100 ml of Hoechst 33342 reaction solution was added to each well and incubated for 30 min. A microscope (DSZ2000X, Beijing Cnmicro Instrument Co., Ltd., Beijing, China) was used to observe and take pictures.

Wound healing assay
Cells (1 × 10 5 /well) were plated in 6-well plates until they achieved about 90% confluence. A 1-ml pipette tip was used to create scratch wounds, and photographs of the wounds (time 0 h) were immediately taken. Then, the cells were cultured in a serum-free medium. After incubation with 5% CO 2 at 37°C for 48 h, photographs were taken again.

Transwell assay
The invasion ability of cells was evaluated using the transwell assay. Transwell chambers (3428, Corning, New York, NY, USA) were pre-cooled overnight at 4°C one day in advance. Then, 100 ml of Matrigel dilute in serum-free medium was added to each well. The transfected cells were suspended in a serum-free medium and added to the transwell chamber; 500 ml of 10% fetal bovine serum complete medium was placed into the lower chamber. The cells were incubated at 37°C for 48 h. The upper chamber was removed and washed three times with phosphatebuffered saline (PBS), and the cells in the upper layer of the membrane were wiped off with a cotton ball. Cells were fixed with 4% paraformaldehyde for 20 min and stained with 0.1% crystal violet solution for 5 min. Cells were observed under an inverted microscope, and three fields of view were taken. The number of invasive cells was recorded.

Statistical analysis
The data were mainly visualized using the R package ggplot2. The Shapiro-Wilk normality test was used to evaluate for normality of variables. For normally distributed variables, significant quantitative differences were determined by twotailed t-tests or one-way ANOVA. For non-normally distributed variables, significant quantitative differences were determined by the Wilcoxon test or the Kruskal-Wallis test. The Benjamini-Hochberg method was used, which converts pvalues to false discovery rate (FDR) to identify significant genes. The log-rank test was used to determine the statistical differences in each dataset. The Survminer R package was used to generate survival curves. Receiver operating characteristic (ROC) curves were drawn using the pROC package. All heatmaps were generated based on pheatmap. All statistical analyses were performed in R (https://www.r-project.org/). p < 0.05 was considered statistically significant.

Construction of pyroptosis-related lncRNA signatures
Forty-five pyroptosis genes were obtained from the literature, and the pyroptosis score was calculated by the GSVA method. The correlation between lncRNAs and pyroptosis scores was analyzed by Spearman's correlation test. The lncRNAs with |correlation coefficient| > 0.3 and p < 0.05 were selected, and 553 lncRNAs were obtained. These lncRNAs were subjected to univariate analysis, and 27 lncRNAs (p < 0.05) were screened ( Figure 1A). The single-factor meaningful genes were selected for survival random forest analysis, and six important lncRNAs were screened ( Figure 1B, Figure S2A). These important gene variables were then used for Lasso analysis to build a risk scoring model ( Figure 1C). The risk score was 0.4242 × AC005253.1 + 0.8249 × LIMD1-AS1 − 1.6456 × LINC01852 + 0.3408 × AC127502.2 + 1.5518 × AC068580.3 − 1.7451 × AC129507.1. With the median risk score as the cutoff value, PCa patients were categorized into high-score and low-score groups. PCa patients in the high-score group had poorer overall survival than those in the low-score group ( Figure 1D). Furthermore, the AUC values of the 1-, 3-, and 5-year ROC curves were 1, 0.93, and 0.92, respectively ( Figure 1E). The clinical characteristic score showed that the risk score of patients aged ≥45 was higher than that of patients aged <45, but there was no significant difference (p = 0.15). Patients with stage N1 had a higher risk score than patients with stage N0 (p = 0.00039). T4 stage patients had higher risk scores than T3 stage and T2 stage patients, and T3 stage patients had higher risk scores than T2 stage patients (p = 0.0014). Moreover, patients with the status alive had a lower risk score than patients with the status dead (p = 0.001) ( Figure 1F). A publicly available dataset (GSE116918) was used to validate the reliability of the constructed risk scoring model. Consistent with the findings from the TCGA-PRAD cohort, survival analysis showed that patients in the high-score group had lower survival rates than those in the low-score group ( Figure S2B). In addition, the AUC values of the 1-, 3-, and 5-year ROC curves were 0.7, 0.71, and 0.77, respectively ( Figure S2C). To investigate whether the constructed risk scoring model was independent of clinicopathological parameters, univariate and multivariate Cox regression analyses were performed on age, T stage, N stage, and risk score. Risk score was the parameter independently predicting overall survival ( Figures S2D, E). The predictive model could be considered an independent prognostic factor in PCa patients.

Correlation of risk scores with pyroptosis genes and immune infiltration
We surveyed the relationship between model genes and risk scores. Risk scores were favorably associated with LIMD1-AS1, AC127502.2, AC005253.1, and AC06850.3 and negatively associated with LINC01852 and AC129507.1 (Figures 2A, B). Next, we constructed a heatmap of risk scores and pyroptosis gene correlations. The results showed that ARRDC1-AS1, GPX4, GSDMD, GSDME, and NLRP3 were substantially associated with the risk score ( Figure 2C). We used the ESTIMATE package and ssGSEA algorithm to evaluate immune infiltration. Among them, Stromal score, Immune Score, and ESTIMATE Score were negatively associated with risk score, while Tumor Purity was positively associated with risk score ( Figure 2D). B cell, T cell, Macrophage, and Natural killer cell were significantly associated with risk score ( Figure 2D). For example, Activated B cell, Activated CD8 T cell, and Activated dendritic cell were negatively associated with risk scores. Activated CD4 T cell was positively associated with risk score.

Immune checkpoint
We determined the expression of seven classes of immune checkpoint molecules in low-and high-risk-scoring populations. As shown in the Antigen present classification, HLA-A, HLA-B, HLA-DPA1, HLA-DPB1, HLA-DQB2, HLA-DRB1, and MICA were expressed at a low level in the high-score group. In Cell adhesion, SELP was expressed at a high level in the low-score group. In Ligand, CCL5, CX3CL1, and TGFB1 were expressed at a high level in the low-score group. In Receptor, CD27, CD40, EDNRB, and TLR4 were expressed at a high level in the lowscore group. In Co-inhibitor, CD276 was expressed at a low level in the low-score group, and PDCD1LG2 and VTCN1 were expressed at a high level in the low-score group. In addition, HMGB1 was expressed at a low level in the low-score group, while ENTPD1 and PRF1 were expressed at a high level in the low-score group (Figure 3).

Functional analysis of risk score
We used the GSVA package for GO and KEGG enrichment analyses. Most samples were enriched for pathways closely related to tumorigenesis. Examples included DNA replication, cell cycle, and mTOR signaling pathway. These pathways were positively associated with risk scores ( Figure 4A). Correlation analysis with risk score and functional enrichment pathway score was performed. In addition to the immune checkpoint, risk scores were significantly associated with 16 other pathways ( Figure 4B). Risk scores were positively correlated with cell cycle, DNA replication, DNA damage repair, and WNT target while negatively correlated with CD8 T effector and antigen processing machinery. We conducted a relevant analysis with the risk score and all genes and then performed a GSEA with the clusterProfiler package ( Figure S3). The Hippo signaling pathway was downregulated.

Drug sensitivity analysis of risk scores in two groups
We downloaded information on the susceptibility of tumor cell lines to potential drugs from the CTRP v2 and PRISM. The lower the AUC of the cell line, the higher the sensitivity to the potential drug. These data revealed that among anticancer drugs, including ML258, 16-beta-bromoandrosterone, VU0155056, BRD-K02251932, BRD-K85133207, imiquimod, temoporfin, SGI-1027, and eptifibatide, the sensitivity of patients in the high-score group to the drugs was significantly increased than in the low-score group ( Figure 6).

Silencing of AC005253.1 promoted pyroptosis of prostate cancer cells
We further tested the effect of AC005253.1 on the pyroptosis of PC-3 and DU145 cells by Western blotting experiment. The results showed that after silencing AC005253.1, the expression of the inflammasomes (AIM2, NLRC4, and NLRP3) was altered ( Figure 9A). Compared with the si-NC group, the expression of Immune checkpoint. Immune checkpoint molecule expression in low-and high-score groups. *p < 0.05. **p < 0.01. ***p < 0.001 ****p < 0.0001. ns, not significant.

Discussion
PCa is one of the most common tumors in men. Due to its heterogeneity and progressive nature, it remains incurable (41). Valid prognostic models based on specific biomarkers can accurately predict survival outcomes for the effective management of PCa patients (42). Pyroptosis-related lncRNA risk prediction models have been reported to be expected to assist in the treatment and management of various tumors (43,44). Our study used a novel risk model of six PRLs (AC129507.1, AC005253.1, AC127502.2, AC068580.3, LIMD1-AS1, and LINC01852) developed using Lasso analysis. AC129507.1, AC005253.1, AC127502.2, AC068580.3, LIMD1-AS1, and LINC01852 were identified as PRLs for the first time. The results showed that the risk score model had the best ability to distinguish clinical characteristics between the high-risk and low-score groups significantly. In our proposed model, the AUC values of the ROC curves for 1-, 3-, and 5-year PCa were 1, 0.93, and 0.92, respectively. In addition, the 1-, 3-, and 5-year AUC values in the test set also had desirable results. Our risk model had excellent predictive power compared to other published pyroptosis-based prognostic models in PCa (45,46). Different from the direct use of Lasso to build a prognostic model of eight pyroptosis-related genes in the study of A B FIGURE 4 Functional analysis of risk score. (A) Heatmap for GO and KEGG analyses using GSVA package. (B) Correlation of risk score with functional enrichment pathways. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSVA, gene set variation analysis. *p < 0.05. **p < 0.01. ***p < 0.001 ****p < 0.0001. A B FIGURE 5 CNV and mutation analysis in high-and low-score groups. (A) CNV maps for groups with high and low score. (B) Top 30 gene mutation frequencies in two groups. CNV, copy number variation. *p < 0.05 Wang et al. (47), we used random forest dimensionality reduction and screening methods and further used Lasso analysis to build a predictive model. This analytical approach may help improve the predictive accuracy of the risk model signature.
It has been reported that pyroptosis is closely related to the tumor immune microenvironment. The release of inflammatory factors is caused by pyroptosis triggers powerful antitumor immunity (48). The infiltration of CD8+ T cells and natural killer cells in the pyroptosis-activated immune microenvironment can promote pyroptosis and form a positive feedback loop (49). Previous studies have shown that increased CD8+ T-cell infiltration is independently associated with improved survival after radical prostatectomy (50). Men with more CD4+ T cells in the prostate tumor environment have an increased risk of dying from PCa (51). B-cell activation is thought to be a driver of the PCa immune response and improves postoperative survival (52). Overall, the poor prognosis and outcome of PCa are closely related to pyroptosis-triggered immune cell infiltration, which is consistent with our results. Our results showed that the risk score was negatively correlated with activated B cells, CD8 T cells, and dendritic cells while positively correlated with activated CD4 T cells. Infiltration levels of B cells, CD8 T cells, and dendritic cells decreased with increasing risk scores, consistent with shorter survival times in patients with high scores. Drug sensitivity analysis of risk scores in two groups. Information on the sensitivity of tumor cell lines to potential drugs was downloaded from CTRP v2 and PRISM. Prediction of drug susceptibility in two groups using pRRophetic. The lower the AUC of the cell line, the higher the sensitivity to the potential drug. AUC, area under the curve. ***p < 0.001 Risk score gene expression identification. RT-qPCR to detect the AC129507.1, AC005253.1, AC127502.2, AC068580.3, LIMD1-AS1, and LINC01852 expression in RWPE1, PC-3, and DU145 cell. *p < 0.05, vs. RWPE1.
We compared the expression of immune checkpoints in the high-risk and low-risk groups and found that some immune checkpoints such as HLA and MICA were expressed at a high level in the low-risk group. PCa downregulated the expression of the HLA-1 antigen processing machinery (APM) and had defects in the antigen presentation pathway (53). Low expression of MICA is associated with poorer overall survival in PCa and is associated with aggressiveness (54). This suggests that our signature could effectively identify the status of immune checkpoints in different PCa patients, providing new ideas for their treatment.
CNVs were regions of the genome with integer copy number changes, including amplifications and deletions of DNA sequences, that could drive cancer's rapid adaptive evolution and progression (55). The CNV results uncovered significant differences in mutation status between the high-score and lowscore groups. In the high-score group, the gene with the highest mutation frequency was TP53 (17%), while in the low-risk group, it was SPOP (11%). TP53 mutation was the most common genetic alteration that played a major role in the pathogenesis of PCa (56,57). SPOP mutations were associated with improved overall survival, whereas TP53 mutations were associated with poorer survival in secondary metastatic hormone-sensitive PCa (58). These data implicated that the high-score group might have more tumorigenic gene mutations.
To better assess the risk model's clinical feasibility, we analyzed information on the sensitivity of tumor cell lines to potential drugs. The results showed that the high-risk group cell lines were significantly less sensitive to drugs such as importazole and imiquimod. Importazole, a specific inhibitor that alters the interaction of KPNB1 with RanGTP, has a good inhibitory effect on PCa progression (59). Imiquimod (also known as a TLR7 agonist) inhibits the growth of mouse (TRAMP C2) and human PCa cells and can be used as an alternative therapy for locally generated PCa (60). Our findings suggest that the low-risk group is more likely to benefit from these drugs. Collectively, these findings may provide prospective treatment options for PCa patients.
Notably, RT-qPCR analysis confirmed high expression of AC129507.1, AC005253.1, AC127502. Frontiers in Oncology frontiersin.org AC068580.3, and LIMD1-AS1 were thought to play important roles in different cancers, while AC005253.1, AC127502.2, and LINC01852 were identified for the first time. AC129507.1 was identified as an on-risk gene in risk models for prognosis patients with gastric adenocarcinoma (61). AC068580.3 was identified as an autophagy-related lncRNA as an indication of prognosis for colon adenocarcinoma (62). LIMD1-AS1 inhibited lung cancer progression by inhibiting cell multiplication and promoting apoptosis (63). Our study demonstrated that inhibiting the expression of AC005253.1 could inhibit cell viability, migration, and invasion. We further detected the expression of the inflammasome (AIM2, NLRC4, and NLRP3) and found that silencing of AC005253.1 could significantly increase the expression of the AIM2 inflammasome. Activation of the AIM2 inflammasome can promote pyroptosis (64). We also found that inhibition of AC005253.1 could promote pyroptosis in PCa cells. Therefore, we speculate that AC005253.1 may affect pyroptosis through the AIM2 inflammasome in PCa. In the present study, we report for the first time the relationship between AC005253.1 and pyroptosis in PCa.
In conclusion, we successfully established an efficient forecast PCa model based on six PRLs, including AC129507.1, AC005253.1, AC127502.2, AC068580.3, LIMD1-AS1, and LINC01852. This well-validated model built on these six PRLs will provide new insights into identifying PCa prognosis. Through in vitro experiments, we verified that silencing of AC005253.1 could inhibit the proliferation, migration, and invasion of PCa cells. In addition, silencing of AC005253.1 might promote pyroptosis by affecting the expression of AIM2 in PCa.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions
JL conceived the project. RT performed the experiments. JY analyzed the data. JY and JL wrote the paper. All authors contributed to the article and approved the submitted version.