Integrative analysis to screen novel pyroptosis-related LncRNAs for predicting clinical outcome of glioma and validation in tumor tissue

Identification of prognostically significant PRLs in glioma patient samples

Graphical abstract showed the overview of this study. We evaluated 14391 lncRNAs by performing the RNA-seq of the glioma patient samples from the Cancer Genome Atlas (TCGA). We also took 11 pyroptosis related genes from the previous literatures. Then we identified 46 PRLs by carrying out Pearson correlation between the lncRNAs and the PRLs applying |R|>0.7 and P < 0.05 as the selective standards (Figure 1A, Supplementary Table 1). These were AC002553.2, AC004817.3, AC004847.1, AC005785.1, AC006369.1, AC007038.1, AC015813.1, AC015961.2, AC027682.4, AC040162.3, AC069281.2, AC073869.1, AC087672.2, AC087741.1, AC092809.4, AC109460.3, AC116049.1, AC135048.4, AC145098.1, AL008729.1, AL031705.1, AL109811.1, AL138831.3, AL355987.4, AL450384.2, AL512770.1, AP002490.1, ARHGAP27P1-BPTFP1-KPNA2P3 (ABK), CARD8-AS1, CCDC18-AS1, CYTOR, FAM13A-AS1, GUSBP11, HCP5, LINC01146, LINC01355, LINC01426, LINC01504, LINC01506, PCED1B-AS1, PDXDC2P-NPIPB14P, PSMB8-AS1, RFPL3S, SLC25A25-AS1, SSBP3-AS1, USP30-AS1. Moreover, we examined the correlation between the expression patterns of these regulators and molecular features. All PRLs were clearly distinguished in tumor grading (Figure 1B). Of the 46 PRLs, 44 revealed significant differences in groups of molecular subtypes, while AC005785.1, and PDXDC2P-NPIPB14P did not (Figure 1C). Of the 46 PRLs, 38 showed significant differences between tumor and normal tissues, while AC006369.1, AC007038.1, AC116049.1, AL008729.1, AL031705.1, AL138831.3, LINC01506 and SSBP3-AS1 did not (Figure 1D). These findings also correspond to the poor prognosis of cluster2 as shown in Figure 2E. Taken together, these results suggested that cross-talk among the PRLs plays an important role in the development of glioma.

Figure 1. Construction and validation of the 46 pyroptosis-related lncRNA (PRLs) signature in Glioma patients. (A) Correlations between PRLs and pyroptosis regulator for glioma (Pearson test). (B) The expression of 46 PRLs among grades in glioma. (C) The expression of 46 PRLs between molecular subtypes. (D) The expression of 46 PRLs between normal tissues and glioma tissues. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Identification of PRLs pattern in glioma

To conduct further studies on the expression features of PRLs in glioma, we classified glioma patients using a consensus clustering algorithm the expression profiles of 46 PRLs. We chose k = 2 for stable clustering of PRLs according to their cumulative distribution functions. Subsequently, we applied two modulation patterns by using the unsupervised clustering method, containing 484 cases in PRLs cluster 1 and 214 cases in PRLs cluster 2 (Figure 2A, 2B). Cluster 1 has a higher survival advantage than cluster 2 (Figure 2E). We further proceeded heatmap and quantitative analysis of the expression values of the 46 PRLs. Among the 46 PRLs, 14 PRLs in cluster 2 were significantly higher than in cluster 1, and the remaining 32 PRLs which were the reverse (Figure 2C, 2D, Supplementary Table 2). Then, we calculated the distribution of grades, IDH and 1p19q in the two clustering groups respectively, and found that G4, IDHwt, and non-codel accounted for a high proportion of cluster 2 (Figure 3A–3C). The strong relationship between our PRLs subgroups and clinical features further demonstrates the precision and robustness of our identification of PRLs patterns in glioma.

Figure 2. Subgroups of glioma related by PRLs. (A) The consensus score matrix of all samples in TCGA cohorts at k = 2. (B) Consensus clustering cumulative distribution function (CDF) for k = 2–6 in TCGA cohort. (C) The expression of 46 PRLs between two cluster groups. (D) The heatmap for 46 PRLs and 22 kinds of immune cells. (E) KM curves for the two cluster groups (Log-rank test, p = 0.0001).

Figure 3. The proportion of patients with grades, IDH1 mutation type, and 1p19q codel status in the high or low-risk groups in the TCGA cohort. (A) The proportion of patients with grades mutation type in the high or low-risk groups. (B) The proportion of patients with IDH1 mutation type in the high or low-risk groups. (C) The proportion of patients with 1p19q codel status in the high or low-risk groups. (D) The expression value of 8 lncRNA in glioma tissues and normal brain tissues.

The immune landscape of PRLs groups

The interplay between PRLs groups and the immune system relies on a complicated cellular cross-talk involving PRLs and immune cells [14]. Therefore, we assessed the immune status of gliomas by ssGSEA base on the infiltration of immune cells in the tumor organization, and classified glioma samples into two subgroups by consensual clustering. The results showed that we differentiated immune cell infiltration between the two cluster groups and there was significant enrichment of immune cell infiltration observed in cluster2 group, the risk score gradually increased with tumor grades (Figure 2D, Supplementary Figures 1B and 3B). Immune checkpoints (ICPs) are crucial for cancer immunotherapy, and a number of ICPs agonists and antagonists being evaluated in clinical oncology [15–17], we further performed their expression profile in distinct subtypes. Except for TNFRSF25, TNFSF18, TNFSF9, and VSIR, all the ICPs were overexpressed in the TCGA cohort high immunity group (Supplementary Figure 1A). Thus, patients in cluster2 may have a stronger immune response to tumorigenesis and tumor progression and therefore may benefit more from immune checkpoint inhibitors therapy than patients in cluster1. Then, we quantitatively analyzed the ESTIMATEscores of the two clusters and found that the Immunescore and Stromalscore of cluster 2 were significantly higher than that of cluster 1, and the tumor purity was just the contrary (Supplementary Figure 1C). The above results confirm the high immune status of cluster2 and the low immune status of cluster1.

The mutation landscape of PRLs groups

Related papers indicate that immunity state may also contact with mutation [18]. Greater tumor mutational burden (TMB) and somatic mutation rates are linked to stronger cancer immunity [19]. We performed variation rate analysis in both groups in order to investigate whether immune infiltration status was associated with mutation rate. Among the two subtypes, cluster1 had the higher mutation rate (97.28%) than that in cluster2 (87.63%) (Figure 4A, 4B). The IDH1 mutation was higher in cluster1 (77%) than cluster2 (18%), IDH1 mutation dramatically indicated the outcome of glioma patients, so the distinction in IDH1 mutation between two cluster subtypes may be one of the reasons that contribute to the survival of patients [20]. This finding is in agreement with the results of the survival analysis of the two clusters in Figure 2E. In addition, we explored the landscape of co-occurrence using the top first 25 mutation genes with the comet algorithm. Twelve pairs of cases (IDH1-IDH2, IDH1-FLG, IDH1-PTEN, IDH1-EGFR, IDH1-NF1, PIK3CA-TP53, FUBP1-TP53, FUBP1-ATRX, CIC-TP53, CIC-ATRX, IDH1-EGFR, IDH1-PTEN) compared with prevalent mutually exclusive mutations showed mutually exclusive mutations, indicating that they may have superfluous impact in the common pathway and a selected advantage of retaining the mutation copy between them (Figure 4C). After checking transcriptional changes in the two subgroups mentioned above, the presence of genomic level distinction between the two subgroups was further investigated. Somatic mutations, comprising single nucleotide polymorphisms (SNP), single nucleotide variants, insertions, and deletions, were computed and shown applying the “maftools” algorithm [21]. The SNPs and Total in cluster2 group were also exceeded by those in cluster1 group, while the major portion of genomic variants were missense mutations (60%) in the two subgroups. As for SNVs, all samples were examined, and C>T was the prevalent type in both groups. For most of the types of SNV (C>A, C>T, T>A, T>G), the mutation number was clearly higher in the cluster2 group than in cluster1 group (Figure 4E). Moreover, the samples in the cluster1 subgroup have a remarkably higher level of variant allele fractions (VAFs) than those in the cluster2 subgroup (Supplementary Figure 2B), which had been thought to be linked to cancer progression and worse results, supporting the discovery of relatively superior tumor purity and inferior heterogeneity in the cluster1 samples. This finding echoes the answer we got in Figure 4C. The proportion between transversion and transition in all SNVs were roughly 3:1 and maintained stabilized in both subgroups (Supplementary Figure 2A). More interesting is that several genes had distinct mutation rates between the two subgroups. In terms of outcomes, the top 10 genes were shown in Figure 4D. Furthermore, IDH1 is another classic example demonstrating the distinct mutation sites between two subgroups (Supplementary Figure 2C) and the plausible chain reaction of the variance in prognostic effect. Eventually, we evaluated driver gene for the two subtypes, and the findings indicated that the dominant driver genes of cluster1 subgroup were IDH1, IDH2, and KRT15, meanwhile, the driver gene of cluster2 subgroup was IDH1 (Supplementary Figure 2D). With removing germline CNV, remarkable gains and losses were observed for each cohort with a threshold FDR < 0.05. Through comparison, we detected more regional alterations in the cluster1 subgroup (Figure 5A, 5B, Supplementary Table 3). By computing the frequency of each CNV in all patients, we discovered that 7q gain and 10p loss were the most common CNVs in cluster1 subgroup; whereas 1p loss, 7q gain and 19p loss were also among the most common changes occurred in cluster2 subgroup (Figure 5C, 5D, Supplementary Table 4). We found significant differences between two subgroups in most mutation features, further demonstrating the accuracy of grouping by PRLs in predicting mutational features.

Figure 4. The mutation landscape of two cluster groups. (A, B) Waterfall diagram displays the mutation landscape of the top 20 most commonly mutation genes. (C) The heatmap analyzes the mutual co-occurrence and exclusion mutations of the top 25 commonly mutation genes. (D) The forest plot shows the top 10 most distinctively mutation genes between the two groups. (E) Boxplots displaying the comparisons of mutation frequencies of each mutation type classified by effects, SNV, DEL, INS, and SNP.

Figure 5. Comprehensive analyses of copy number variation between two cluster groups. (A, B) Detection and comparison of the percentage of significant scores between the two groups. (C, D) Detection and comparison of the percentage of significant gains and losses between the two groups.

Identification hub PRLs in glioma and construction of pyroptosis related LncRNA prognostic model

To further obtain more accurate PRLs, we conducted univariate and multi-variate cox regressions from 46 PRLs and yielded 10 hub PRLs (Figure 6D). The 10 hub PRLs were AC004817.3 (lnc-MAP3K9-10), AC007038.1, AC040162.3, AC145098.1 (PRR7-AS1), AL450384.2 (NSUN6), CYTOR, FAM13A-AS1, HCP5, PSMB8-AS1, USP30-AS1. To verify the reliability of these 10 hub PRLs for glioma patients, we validated these 10 PRLs in human glioma samples. AC007038.1 and AC040162.3 were not found with corresponding sequences due to their relatively short discovery time, and we performed PCR analysis on the remaining 8 lncRNAs, which showed consistent expression with the TCGA database (Figure 3D). Seeing that this has the stability of lncRNA expression, it lays the foundation for our further analysis. We constructed pyroptosis related lncRNA prognostic model (PRLPM) by lasso regression for 10 hub PRLs (Figure 6A, 6B). Moreover, we built the riskscore system according to PRLPM to classify patients into high and low risk groups based on PRLPM (Supplementary Table 5). Kaplan Meier survival analyses revealed that overall survival (OS) was worse in the high-risk cohort in comparison to the low-risk cohort (Figure 6E), the same result was obtained in Chinese Glioma Genome Atlas (CGGA) cohort (Supplementary Figure 3F). The survival preference in the low-risk group was superior to the high-risk group, whatever treatment was used (Figure 6F, 6G). Meanwhile, the timeROC showed the mean AUC values of prognosis predictions on TCGA cohort reached 0.86, 0.87, 0.82, 0.82 and 0.75 (Figure 6H), the AUC values on CGGA were 0.7, 0.7, 0.69, 0.71, and 0.79 (Supplementary Figure 3E). Distribution of PRLPM riskscore and OS also verified the survival analyses, the findings suggested that the expression of 10 hub PRLs was remarkably positively associated with PRLPM riskscore and OS (Figure 6C). From the above results, we found that PRPLPM has high accuracy in predicting the risk profile and prognosis of glioma patients.

Figure 6. Predicting patient prognosis in the TCGA cohort based on PRLPM. (A) Regression coefficient profiles of identified pyroptosis immune regulators in the TCGA cohort. (B) Ten-time cross-validation for tuning parameter selection in the TCGA cohort. (C) Patients were divided into high and low-risk subgroups according to median level of PRLPM riskscore in train set; heatmap of 9 PRLs. (D) multivariate cox regression analyses of the association between PRLs and OS of patients in the TCGA cohort. (E) KM curve plot of OS for patients in high and low-risk subgroups. (F, G) Survival analyses for subgroup patients stratified by both PRLPM riskscore and treatment with radiotherapy (F) and pharmacological chemotherapy (G) in the TCGA cohort. (H) The timeROC curve to evaluate the prognostic value of PRLPM riskscore in TCGA cohort.

The relationship between riskscore and clinical traits

We ordered the samples on the basis of their riskscore (rank low to high riskscore) and examined whether any demographic/molecular/clinical characteristics were associated with PRLPM in two public cohorts (Supplementary Figures 4A and 3A, Supplementary Table 6). We further quantitatively analyzed the relevance between risk scores and molecular/clinical traits. In the TCGA cohort, the results revealed that the PRLPM riskscore gradually increased with increasing grade, with lower PRLPM riskscore in patients with IDH1mut than that in IDH1wt patients and significantly higher PRLPM riskscore in 1p19q non-codel patients than that in 1p19q codel patients (Supplementary Figure 4B–4D). In the CGGA, the findings revealed that the PRLPM riskscore gradually increased with increasing grade, with lower PRLPM riskscore in patients with IDH1mut than that in IDH1wt patients and significantly higher PRLPM riskscore in recurrent and secondary patients than that in primary patients (Supplementary Figure 3B–3D). The proportion of tumor grading within the two subgroups showed that G4 and G3 were only found in the high-risk group. The proportion of IDH1mut in both subgroups revealed that IDH1mut was in predominant in the low-risk group. The rate of 1p19q codel was considerably higher in the low-risk than in the high-risk group (Supplementary Figure 4E–4G). These results demonstrate that PRLPM riskscores can respond to the demographic/molecular/clinical feature of glioma.

The relationship between riskscore and tumor stemness indices

Stems differed between PRLPM groups; specifically, the dedifferentiated phenotype was significant in the high-risk group, while the differentiated phenotype was significant in the low risk group (Supplementary Figure 5A, 5B). Then, we calculated the TMB and found that the TMB in the high-risk group was obviously higher than that in the low-risk group. Moreover, mRNAsi, EREG-mDNAsi, mDNAsi, EREG-mRNAsi, ENHsi, TMB, and DMPsi were actively and significantly related to the PRLPM riskscore, and the relevance were significant at 0.45, 0.49, 0.5, 0.018, 0.5, 0.45, 0.39, respectively (Supplementary Figures 5C–5H, 6C, Supplementary Table 7). Considering that the PRLs characterized an immune activation phenotype of glioma, we applied the algorithms to evaluate the effectiveness of the PRLs signatures in forecasting ICB responsiveness in glioma and used SubMap to compare the prediction results [22, 23]. Consequently, our PRLs had comparable performance in predicting the glioma response to anti-PD1 and anti-CTLA4 treatments (Supplementary Figure 6D), indicating that stratification based on PRLs has the possible to recognize ICB reactors. In conclusion, there were distinct differences in the extent of tumor differentiation, and TMB, between the PRLPM groups.

Comparison of differentially expressed genes and sensitivity to chemotherapeutic agents in high and low-risk groups

We examined the differential analysis for the both groups. Enrichment analysis was performed to forecast the potential functions of differentially expressed genes (DEGs) between high and low-risk groups, and as expected, the DEGs were concentrated in carcinogenic and immune pathways, for example neutrophil degranulation, neutrophil activation involved in immune response, wnt signaling pathway, p53 signaling pathway, Apoptosis, HALLMARK_APOPTOSIS, PI3K_AKT_MTOR_SIGNALING, HALLMARK_COMPLEMENT (Supplementary Figure 7A–7D).

Different PRLPM subgroups should guide clinical therapy. Therefore, we made a comparison of the sensitivity of 30 anticancer drugs between the both groups to identify potential glioma treatment approaches. A total of 22 chemotherapy drugs had remarkably distinct IC50 estimates between two groups (Supplementary Figure 6A). Patients in the low-risk group may be sensitive to these drugs. In this context, these agents have the potential to be used in the future for the treatment of low-risk PRLPM. Moreover, we find that the TMB in high-risk group was obviously higher than low-risk group (Supplementary Figure 6B). This result indicated high risk had better therapeutic effect than low risk in immunotherapy.

The efficacy of the PRLPM riskcore across tumor types

Considering the strong correlation between PRLPM and immunity described above, we further examined the validity of the PRLPM scoring system in different tumor types. The percentage of 22 immune cells was characterized as immune infiltration. We computed the relevance between the 22 immune cells percentage and the PRLPM riskcores and revealed the distinct trends in the relevance for pan-cancer, except PCPG and READ. The percentage of T cells CD8, T cells CD4 memory active, T cells CD4 memory resting, Macrophages M1, B cells naive were relevant to the PRLPM riskcores of most cancer types. Interestingly, these cells belonged to the antitumor type, to some extent, suggesting that the PRLs promote tumor immunity (Figure 7A). Furthermore, we identified a correlation between PRLPM riskscores and stem cell indices for all tumors except for ACC, KICH, SKCM, and UCS (Figure 7B). We found the significant relationship between the ESTIMATEscore and the PRLPM riskcores in addition to PCPG, and READ (Figure 7C). Such cancers with high mDNAsi were expected to be hypersensitive to immune checkpoint treatment. We thoroughly assessed the correlation of the PRLPM scoring system and the pan-cancer, and deeply explored the outstanding features of the PRLPM riskscore in pan-cancer, and laid the basis for the wide utilization of PRLPM riskscore in tumor.

Figure 7. Performance of PIPM riskscore across tumor types. (A) Association between the PIPM riskscore and immune cells for each cancer type. (B) Association between the PIPM riskscore and stemness indices for each cancer type. (C) Correlations between the PIPM riskscore and ESTIMATEscore for each cancer type.

Data extraction

The RNA_seq and clinical data of 698 glioma samples were obtained from TCGA, clinical data of 1018 glioma samples were downloaded from CGGA. Mutation rate and CNV frequency were gathered from cBioPortal. Clinical data of glioma patients, including age, sex, grades, OS and survival status was also collected for subsequent analyses.

Recognizing pyroptosis regulators

Access to relevant literature, we collected 11 genes (CASP1, CASP3, CASP4, CASP5, CASP8, GSDMB, GSDMC, GSDMD, GSDME, GZMA, GZMB), which are closely related to pyroptosis and serves as pyroptosis genes [35–37].

Consensus clustering

Consensual clustering utilizes the k-means technique to determine specific pyroptosis patterns associated with the expression of pyroptosis regulators. Both stability and number were decided by the consensus clustering algorithm applying the “ConsensusClusterPlus” package [38].

Immune infiltration

The relative infiltration of 28 immune cell types in TME was described using ssGSEA. Signature gene panels were obtained for each immune cell type from a recent paper [39].

Analyses of mutation subtypes

Somatic mutation and CNV profiles were downloaded from the TCGA data portal. Analysis of somatic mutation data sorted by mutation annotation format was performed using the R package “maftools” [40, 41].

PRLs prognostic model predicts effective response to chemoradiotherapy, the analyses of relationship between clinical characteristics, tumor stemness indices, tumor mutation burden, and prognostic model riskscore

The prognostic model was developed by lasso regression for PRLs, and we analyzed the survival, ROC, risk factor, stemness, the response to chemoradiotherapy. The Kruskal-Wallis test was used to evaluate the difference between TMB, tumor stemness index and the prognostic model generated by lasso algorithm, and Pearson correlation was used to evaluate the correlation between tumor stemness index, TMB and prognostic model riskscore. Furthermore, the dedication to OS were computed by Kaplan-Meier algorithm. We also continued to analyze the association of prognostic model riskscores with clinical characteristics separately.

Significance of the pyroptosis related lncRNA prognostic model in drug sensitivity

To estimate PRLPM in the clinical treatment of glioma, we used the method developed by Geeleher et al. [42] and the corresponding R package ‘prophetic’ to compute TCGA items for the IC50 of widely applied chemotherapy agents in the glioma dataset [43, 44]. The AJCC guidance suggested 30 common anti-tumor medications, like Cisplatin, Vinblastine, Imatinib, Adriamycin and/or cancer therapy. Wilcoxon signed rank method was employed to compare the variation in IC50 values of common antineoplastic agents between PRLPM high- and low-risk groups, and the results are presented as box plots.

Real-time quantitative PCR

The RNA of 6 human glioma tissues and corresponding peritumoral tissues (Normal tissue) were reverse-transcribed by RT reagent Kit gDNA Eraser (TaKaRa), and SYBR-Green (TaKaRa) detection was performed. The PCR primers were shown in Supplementary Table 8.

Statistical analysis

All data were performed with R version 4.0.5 and its corresponding packages. ESTIMATEscore was computed by using the ‘estimate’ package [45]. The lasso cox regression was carried out using the ‘glmnet’ package [46]. The data were analyzed using proper statistical criterion. Multiple-test adjustment was applied utilizing the FDR method.

Availability of data and materials

We organized the original data and uploaded it to Github. This is the URL we uploaded: https://github.com/shuaima1991/pyroptosis-related-lncRNA.git.