Identification of pyroptosis-related signature for cervical cancer predicting prognosis

Datasets and preprocessing

The TCGA database (https://portal.gdc.cancer.gov/) and the UCSC Xena database (https://xena.ucsc.edu/) were used to obtain CC data. Our research program was shown in Supplementary Figure 1. The RNA sequence data, clinicopathological parameters, and genomic mutation data (including somatic mutations and copy number variants (CNV)) were obtained for 306 CC tissues and 3 paraneoplastic tissue samples. Based on earlier studies [9, 12, 13], 35 PRGs were included, namely AIM2, CASP1, CASP3, CASP4, CASP5, CASP6, CASP8, CASP9, ELANE, GPX4, GSDMA, GSDMB, GSDMC, GSDMD, GSDME, GZMA, GZMB, IL18, IL1B, IL6, NLRC4, NLRP1, NLRP2, NLRP3, NLRP6, NLRP7, NOD1, NOD2, PJVK, PLCG1, PRKACA, PYCARD, SCAF11, TIRAP, and TNF, as detailed in Supplementary Table 1. Comparison of mRNA expression levels of PRGs in paraneoplastic tissue and CC tissues, and visualization of somatic mutations and CNV in CC tissues by the R software ‘maftools’ and the ‘Rcircos’ package. The Kaplan Meier-Plotter online platform (https://kmplot.com/analysis/) was used to further analyze the survival and prognosis of PRGs [14]. All subsequent statistical analyses were performed via R version 4.0.2. The screening criterion for differential expression of PRGs was P<0.05.

Establishment of prognostic signature

The data for 304 CC patients were obtained from the TCGA database and used as the training dataset to screen for prognosis-related PRGs and exclude samples with no survival time. We constructed a prognostic signature for the PRGs using a univariate Cox proportional regression analysis and the least absolute shrinkage and selection operator (Lasso) regression, which was stratified according to the risk score (Risk score = EXP _{PRGs 1} × Coefficient₁ + EXP _{PRGs 2} × Coefficient₂ + ... + EXP _{PRGs n} × Coefficient_n). Subsequently, the associated risk score was calculated for each patient in the training dataset. Based on the median score, the patient’s data were grouped as low or high risk. The difference in survival between the two groups was determined using the log-rank test and the Kaplan-Meier curves, while the sensitivity and specificity of the prognostic signature were determined from the ROC curves plotted using the ‘survivalROC’ package. To determine the feasibility and reliability of the proposed prognostic signature, the training dataset of CC patients (N = 304) was randomly divided into test datasets A (N = 152) and B (N = 152) using the ‘Caret’ package. Subsequently, the same statistical methods were used for validating both test datasets.

Then, independent prognostic analyses and decision curve analysis (DCA) were performed to combine the clinical features, and a nomogram was constructed to predict the one-, three- and five-year survival rates of CC patients to assess the prognostic value of the signature.

The TCGA pan-cancer transcriptional expression data and clinical data were obtained from the Xena database that includes 33 cancer types, as detailed in Supplementary Table 2. The prognostic model was applied to the remaining 32 tumors to analyze their prognostic outcomes for the different cancers. The risk score was considered to be significantly associated with CC patient survival when p<0.05.

Tumor mutation burden, tumor microenvironment, and immune cell infiltration analysis

Tumor mutation burden (TMB) is a biomarker that can predict how a patient might respond to immunotherapy [15]. Differences in TMB between the high and low-risk groups were compared based on original data of somatic mutations in CC patients, and the risk scores and TMB were compared using the Spearman correlation analysis. An ESTIMATE algorithm was used to predict the purity of the CC [16], and the proportion of stromal cells (stromal score) and immune cells (immune score) that have infiltrated the tumor tissue as well as the measured differences in the tumor microenvironment (TME) between the two risk groups were assessed.

Immune cell infiltration data from CC were obtained from the Timer 2.0 database (http://timer.cistrome.org/) [17]. Four algorithms were eventually chosen [18–20], including TIMER, CIBERSORT, QUANTISEQ, and MCPCOUNTER, according to the PRG signature, and the correlation between the immune cell infiltration and risk score under different algorithms was analyzed using Spearman correlation and represented as a bubble plot. In addition, the differences in immune function, HLA gene expression, immune checkpoint, and m6A gene expression were compared for the high- and low-risk groups using the Wilcox test. P < 0.05 was considered to be statistically significant.

Construction of mRNA-miRNA-lncRNA interaction network

To screen for the relevant genes in the prognostic signature, we searched for protein-protein interactions (PPI) via the String (https://string-db.org/) online database. The degrees of connectivity of these PRGs were also counted. Next, the correlations of genes in the PRG signature were evaluated concerning the clinical feature ‘Grade’. The patients were divided into two subgroups, namely, G1-2 and G3-4, and the Wilcox test was used for their analysis. The genes with the highest connectivity and statistically significant correlation with ‘Grade’ in the signature were designated as the core genes.

Then, pan-cancer differential analysis, survival analysis, and TMB correlation analysis were performed on these core genes that were designated in the previous step. We used boxplots to summarize the overall levels of expression of the core gene in the 33 TCGA cancer samples. Based on 18 types chosen from more than five adjacent normal tissues, a linear mixed-effect model was used to compare and analyze the difference in gene expression between the tumor and the paracancerous tissues. The TMB scores were obtained from UCSC Xena database. Correlation analysis between the core gene expression and TMB was performed using Spearman’s method. A forest plot was drawn to determine whether the expression of the core gene correlated with the patient survival obtained by the univariate Cox analysis.

To construct the mRNA-miRNA-LncRNA interaction network for the core genes, the miRNAs targeted by the core genes were predicted using the TargetScan (http://www.targetscan.org/vert_72/) and miRTarBase (http://mirtarbase.cuhk.edu.cn/) databases. The predicted miRNAs were intersected. Subsequently, the differential expression levels and prognostic values of these miRNAs were analyzed. Based on the prognosis and differentially expressed miRNAs, LncBase Predicted v.2 (https://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2/index-predicted) and StarBase (http://starbase.sysu.edu.cn/) were used to predict the LncRNAs that the miRNAs might eventually bind. The differential expression levels and prognostic value of these LncRNAs were analyzed in detail, and on this basis, we constructed the mRNA-miRNA-LncRNA interaction networks. The statistically significant difference was P < 0.05.

Availability of data and materials

The information of this study is obtained by the TCGA, UCSC Xena, Kaplan Meier-Plotter, Timer 2.0, String, TargetScan, miRTarBase, LncBase Predicted v.2, and StarBasedatabase database. We are grateful to them for the source of data used in our study.

Genetic variation profile of PRGs in CC

First of all, we found that the AIM2, CASP3, CASP4, CASP5, CASP6, CASP8, GSDMB, GSDMC, GZMB, IL18, NLRP2, NLRP7, NOD2, PYCARD, and TNF were significantly up-regulated in CC tissues, while ELANE, NLRP1, NOD1, and PJVK were significantly down-regulated (Figure 1A). A survival prognosis analysis showed that IL 1B and PRKACA had a significant correlation with the prognosis of CC patients. Those patients with low expression of IL 1B and high expression of PRKACA have a better prognosis for survival (Figure 1B, 1C). Figure 1D shows the locations of the CNV alterations on the PRG chromosomes, and it can be observed that CNV changes were prevalent in these genes. PRKACA had the highest probability of CNV amplification and GPX4 had the highest frequency of CNV deletion (Figure 1E). The somatic mutations of these PRGs in CC were also analyzed and it was concluded that the overall mutation frequency was not high in these cases (≤4%) (Figure 1F).

Construction and evaluation of the effectiveness of prognostic signature

A survival analysis of the PRGs in the training dataset was conducted via a univariate Cox proportional regression model and screened for eight PRGs that had prognostic value (Figure 2A). Among these, GZMB and TNF were expressed in high level in CC tissue while NOD1 was significantly low (Supplementary Figure 2B).

Following this, the pyroptosis-related signature was constructed by further downscaling the prognosis-related PRGs via Lasso regression (Figure 2B, 2C). The risk score was calculated for each patient using the expression, risk score= (-0.0037×EXP_GPX4) + (0.0377×EXP_GSDME) + (-0.0043×EXP_GZMA) + (-0.026×EXP_GZMB) +(0.0017×EXP_IL1B) + (0.4741×EXP_NOD1) + (-0.1298×EXP_PRKACA) + (0.0642×EXP_TNF). According to the median, the patients were divided into low-risk and high-risk groups, and the survival prognostic analysis showed a significant survival benefit for patients classified as low-risk (P < 0.001, Figure 3A), with an AUC value of 0.789 (Figure 3D). In test datasets A and B, the difference in survival prognosis (P < 0.05, Figure 3B, 3C) and predictive efficacy (AUC values of 0.707 and 0.943, respectively, Figure 3E, 3F) were validated similarly, thereby confirming the overall accuracy and validity of the prognostic signature. Interestingly, this prognostic model indicated a better prognosis for patients in the low-risk group than in the high-risk group in glioblastoma multiforme (GBM), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), and uterine corpus endometrial carcinoma (UCEC) (Supplementary Figure 3).

Univariate and multivariate COX analyses in conjunction with the patient age, stage, and grade, showed that this pyroptosis-related signature was an independent prognostic factor for OS in CC patients (p < 0.001, Figure 4A, 4B). In addition, the risk score corresponded to the highest AUC value (AUC = 0.794, Figure 4C) compared to other clinical features. The DCA further demonstrated the clinical usefulness of the signature (Figure 4D). The Nomogram for both the signature and clinical features was stable and accurate, and it can be used to predict the one-year, three-year, and five-year survival rates in CC patients (Figure 4E).

Correlation of TMB, TME, and immune cell infiltration with prognostic signature

Previous studies have indicated the crucial role of TME in the occurrence and development of tumors [21]. As TMB predicts patient response to immunotherapy, our study found that patients in the high-risk group typically had low TMB scores, and the risk score was negatively correlated with the TMB score to a significant extent. The results of the survival analysis showed that patients that scored high on risk and low on TMB possessed a significant survival advantage (Figure 5A–5C). Compared to the high-risk group, the ESTIMATEScore, ImmuneScore, and StromalScore were all higher in the low-risk group (Figure 5D).

The relation between the immune response and risk score based on the TIMER, QUANTISEQ, MCPCOUNTER, and CIBERSORT algorithms is shown in Figure 5E. The risk score was negatively correlated to the infiltration levels of B cells, T cell CD8+, and macrophages, and positively correlated to the infiltration level of neutrophils. The analysis of immune function differences showed that mechanisms such as the antigen-presenting cell (APC) co-inhibition/ APC_co_stimulation, CC chemokine receptor (CCR), check-point, cytolytic_activity, human leukocyte antigen (HLA), Inflammation promotion, and T cell co-inhibition/ stimulation were significantly more active in the low-risk group (Figure 6A).

Considering this analysis of the immune function, the levels of HLA-related gene expression l were further studied as shown in Figure 6B. The levels of expression of HLA-DMA, HLA-DQB1, HLA HLA-DMA, HLA-DQB1, HLA-DRA, etc. were significantly higher in the low-risk group. As checkpoint inhibitors are crucial to immunotherapy, the differences in immune checkpoint expression between the two groups were analyzed in detail. Significant differences were demonstrated in the expression of PDCD1LG2, TMIGD2, CD27, CD40LG, etc. for the two groups (Figure 6C). Recent research suggests that m6A can regulate the immune response and also modify the immune microenvironment of the tumor [22]. On comparing the m6A-related expression of mRNA between the high- and low-risk groups, it was seen that YTHDF3 and LRPPRC ad significantly high levels of expression in the high-risk group (Figure 6D).

Construction of mRNA-miRNA-lncRNA interoperability network

The interrelationships between the 35 PRGs considered in this study were retrieved from the STRING database and a PPI network having a composite score of > 4.0 was constructed. A complex, close-knit network of interrelationships was found among these genes was determined (Supplementary Figure 4A). Subsequently, the degree of connectivity between these genes (Supplementary Figure 4B). We further examined the relationship between levels of mRNA expression in the pyroptosis-related signature and CC grade staging, and found that the expressions of GPX4, GZMA, and GZMB were significantly higher in the G3-4 group compared to the G1-2 group (P < 0.05, Supplementary Figure 4C–4F). Moreover, the degree of connectivity of GPX4, GZMA, and GZMB in the PPI network above was 4, 4, and 10, respectively, suggesting that GZMB may indeed be involved in tumor growth and progression in CC patients.

Pan-cancer expression analysis was performed for GZMB, which showed its significantly high expression in cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), and UCEC tissues, while significantly low expression was observed in lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) tissues (Figure 7A). TMB is a biomarker that indicates response to immunotherapy, and so, the relation between the expression of GZMB and TMB was investigated. The levels of expression of GZMB were significantly and positively correlated with TMB in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma, and endocervical adenocarcinoma (CESC orCC), CHOL, COAD, LGG, LUAD, STAD, and UCEC, Meanwhile, the correlation was significantly negative with kidney renal papillary cell carcinoma (KIRP) and thyroid carcinoma (THCA) (Figure 7B). Interestingly, the cases of CHOL, COAD, STAD, and UCEC in the high-GZMB-expression group also had high TMB values, and these are expected to be more sensitive to immunotherapy. In contrast, LUAD in the low expression group had low TMB values and so, may be resistant to cell scorch-related immunotherapy. Prognostic correlation analysis showed that the GZMB was negatively correlated with survival in BLCA, BRCA, CESC, ovarian serous cystadenocarcinoma (OV), skin cutaneous melanoma (SKCM), UCEC, and uveal melanoma (UVM). Meanwhile, it was positively correlated with GBM, KIRC, KIRP, acute myeloid leukemia (LAML), and LGG (P<0.05, Figure 7C).

To elucidate the potential molecular regulation mechanism of GZMB in CC, we further constructed an mRNA-miRNA-LncRNA interaction network. First, based on the comprehensive prediction using the TarBase and mirTarBase databases, we screened out miRNAs that can target and regulate the GZMB core gene, and subsequently, ten miRNAs were obtained by considering the intersection (Figure 8A). Further analysis showed that miR-378a was significantly upregulated in CC (Figure 8B) and patients who show a high expression have a better prognosis for survival (Figure 8C). The LncRNA targets upstream of miR-378a were obtained from the LncBase and StarBase databases, and a total of five LncRNAs (Figure 8D) were obtained by intersection to finally construct the mRNA-miRNA-LncRNA interaction network (Figure 8E). Subsequently, these five LncRNAs were subjected to differential expression analysis and survival prognosis analysis, which showed that TRIM52-AS1 was significantly downregulated in CC (Figure 8F), and its high level of expression indicated a considerable survival advantage (Figure 8G). Thus, it may be concluded that the GZMB/miR-378a/TRIM52-AS1 regulatory axis plays an important regulatory role in the development of CC.