Identification of Pyroptosis Genes in cervical cancer and construction of a prognostic model for pyroptosis-associated mRNA

doi:10.21203/rs.3.rs-2110804/v1

Download PDF

Research Article

Identification of Pyroptosis Genes in cervical cancer and construction of a prognostic model for pyroptosis-associated mRNA

https://doi.org/10.21203/rs.3.rs-2110804/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Cervical squamous cell carcinoma and endocervical adenocarcinoma(CESC) is one of the more common tumors in women worldwide and has a higher mortality rate. However, there is a paucity of information about specific biomarkers that assist in the diagnosis and prognosis of CESC. The development of a specific prognostic model is important if we are to improve treatment strategies. Pyroptosis is a form of programmed cell death, and its different elements are related to the occurrence, invasion and metastasis of tumors. However, the role of pyroptosis in CESC progression has not been clarified. The focus of this study is to use comprehensive bioinformatics to develop pyroptosis prognostic characteristics of CESC, so as to delineate the relationship among this characteristic, tumor microenvironment and immune response of patients. In combination with clinical characteristics, risk score is an independent predictor of OS in patients with CESC. Pyroptosis Genes(PRG) score was significantly correlated with immune score, immune infiltration, immune microenvironment, cancer stem cell (CSC) index, and chemotherapeutic drug sensitivity. These findings may improve our understanding of PRGs in CESC and provide new avenues for assessing prognosis and developing more effective immunotherapeutic strategies.

Cervical cancer

Pyroptosis

Recurrence-free survival

TCGA and GEO

Tumor microenvironment

Immune

Cervical cancer (CC) is the fourth most commonly diagnosed cancer and the fourth leading cause of cancer death in women, and is a major global health problem. While the incidence and mortality rates of cervical cancer have declined significantly in high-income countries, the overall prognosis for CC patients in low-income countries remains poor[1, 2]. Despite the push for human papillomavirus (HPV) vaccination and screening programs, there are an estimated 604,000 new cases and 342,000 deaths from CC worldwide in 2020, most of which occur in developing and less developed countries[3]. Persistent infection with high-risk HPV is a major causative factor in the development of CC[4]. The main treatments for CC include surgery, chemotherapy, and radiation therapy. While patients with early-stage tumors can often be successfully treated with modern therapies, their condition may worsen and have a poor prognosis due to recurrent or metastatic disease[1]. Due to the lack of precise diagnostic markers and limited treatment options, new prognostic models urgently need innovative biomarkers to guide targeted therapy and improve clinical outcomes in patients with CC.

Pyroptosis, also known as inflammatory necrosis of cells, is a new form of programmed cell death[5]. Pyroptotic cells are characterized by cell swelling and many bubble-like protrusions. Electron microscopy shows that pyroptosis cells first form a large number of vesicles. After these vesicles form, pores form in the cell membrane that rupture and the contents flow out[6]. The gasdermin family is a major performer of cell pyroptosis and includes gasdermin-A to gasdermin-E and pejvakin (PJVK or DFNB59). Gasdermin family proteins can be cleaved and polymerized, resulting in cleavage of N-terminal and C-terminal junction domains and release of activated N-terminal regions; these regions bind to membrane lipids, phosphatidylinositol, and cardiolipin and localize to cell membranes in the pores[7]. Cellular gasdermin family proteins form 10 to 20 nm pores in the cell membrane, and cell contents are slowly released through the membrane pores and trigger amplified inflammatory responses. Cells gradually flatten, producing 1–5 µm apoptotic vesicle-like protrusions (scorched vesicles), and cells gradually swell until the plasma membrane ruptures, with features such as nuclear condensation and chromatin DNA breakage[8, 9]. Pyroptosis was initially identified as a key mechanism in the fight against infection, and a growing body of research suggests that it also plays an important role in tumor development. Inflammative vesicles, gasdermin proteins and pro-inflammatory cytokines have been reported to be key components of cell pyroptosis and are associated with tumorigenesis, invasion and metastasis[10]. Dupaul-Chicoine et al[11]. Several studies have shown that the potent pro-inflammatory effect of cell pyroptosis is related to the regulation of the tumor immune microenvironment. Defective GSDMD expression was found to be accompanied by a significant decrease in the number and activity of CD8 + T lymphocytes[12]. A recent study also confirmed the key role of cell pyroptosis in the anti-tumor function of NK cells[13]. To date, a number of prognostic pyroptosis-related features have been developed in different tumor types, including ovarian cancer[14], gastric cancer[15], lung adenocarcinomas[16], and head and neck squamous cell carcinomas[17];however, the role of pyroptosis in cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) requires further study.

Therefore, we conducted a systematic study to identify and investigate the genes associated with pyroptosis in normal cervical tissues and CESC tissues, and then we developed a novel prognostic feature by analyzing two pyroptosis related subsets and investigated its association with clinicopathological features and overall survival (OS) in patients with CESC. In addition, signaling pathways, tumor immune microenvironment (TIME), tumor mutation burden (TMB) and immunotherapy response related to CSCC characteristics were also explored. This study provides a better understanding of the mechanisms of pyroptosis and more precise prognostic targets for patients with CESC.

Sources of Cervix Cancer Datasets and Preprocessing

RNA sequence transcriptome data and relevant clinical information from 304 patients with CESC and 3 normal adjacent tissue samples were obtained from The Cancer Genome Atlas (TCGA) database(https://portal.gdc.cancer.gov/). Moreover, somatic datasets and copy number variation (CNV) data for CESC were also downloaded from TCGA and the University of California, Santa Cruz (UCSC) Xena website(https://genome.ucsc.edu/),respectively. Data analysis was performed with the R (version 4.1.0) and R Bioconductor packages. The expression data were normalised to transcripts per kilobase million (TPM) values before further analysis. RNA sequencing data and clinical information from 48 patients with CESC were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/, GSE30759). The clinical information of the CESC patients is shown in Table S1.

Mutation Analysis of PRGs and PRG_score

Based on the findings of previous studies[18–21], 52 Pyroptosis Genes (PRGs) were extracted and the mutation frequency and oncoplot waterfall plot of 52 PRGs in CESC patients were generated by the “maftools”[22] package. The location of CNV alteration of 52 PRGs on 23 chromosomes was drawn using the “RCircos” package in R. The same, to determine somatic mutations between high-risk and low-risk groups of cervical cancer patients, a mutation annotation format (MAF) was generated from the TCGA database using the "maftools" R package. We also calculated the tumor mutation burden (TMB) score for each cervical cancer patient in both groups and analyze the correlation between the two.

Identification of Differentially Expressed Pyroptosis Genes and Functional Analysis

And their mRNA expression levels were compared between CESC and normal adjacent tissue samples in TCGA cohort. Differentially expressed PRGs were identified using the “limma” and “reshape2”package, with P < 0.05 as the significance threshold, and presented as a boxplot. A protein-protein interaction (PPI) network of DEGs was constructed using Search Tool for the Retrieval of Interacting Genes (STRING), version 11.5 (https://string-db.org/). In addition, Gene Ontology (GO), including the biological process (BP), cellular component (CC), and molecular function (MF) categories, was conducted with the “clusterProfiler” package in R software about these DEGs. Similarly, this package was also utilised to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.

Analysis of Pyroptosis Genes in Prognosis

Combine the expression matrix of TCGA database and GEO database by the “limma” and remove batch effect by using the "sva" packages. Kaplan–Meier analysis was used for screening out Pyroptosis Genes which were related to the overall survival (OS) by "survminer" package that can realize visualization. And then, construct the prognosis network of PRGs using "igraph", "psych", "reshape2" and "RColorBrewer" packages.

Tumor Classification Based on the PRGs

The “ConsensuClusterPlus” R package (parameters: clustering algorithm, k-means; distance, Euclidean) was applied, with 1000 iterations[23]. The optimal k value was confirmed, according to cumulative distribution function and delta area values, and patients with CESC then classified into different clusters, based on differentially expressed PRGs.

Identification of Differentially Expressed Genes Among Clusters

The “limma” R package was used to identify DEGs between different pyroptosis-related clusters, with a threshold of |log₂FC| ≥ 1 and P < 0.05. Intersecting DEGs were assembled in both TCGA and GEO cohorts.

Consensus Clustering of DEGs

After univariate Cox regression analysis of DEGs, consensus clustering was applied to identify DEGs by the k-means method. The number of clusters, and their stability, were determined by the consensus clustering algorithm using the “ConsensuClusterPlus” package. We performed 1,000 times repetitions to guarantee the stability of our classification.

Development and Validation of the Gene Prognostic Model Based on PRGs

The data for 350 CESC patients were obtained from the TCGA and GEO database and used as the training dataset to screen for prognosis-related genes and exclude samples with no survival time. The prognostic genes accepted Least absolute shrinkage and selection operator (LASSO) Cox Regression using the R package “glmnet”[24], narrowing the number of mRNAs. Then,multivariate analyses were employed to probe the involvement of the candidate mRNAs for the prognosis of patient OS. On the basis of the least Akaike information criterion (AIC) value, 5 robust mRNAs linked to pyroptosis were identified. This entailed the computation of the risk score for each patient employing the expression level of these mRNAs subjected to normalization and corresponding regression coefficients that is risk score = βgene(1) × EXPgene(1) + βgene(2) × EXPgene(2) +⋯ + βgene(n) × EXPgene(n). Here, EXPgene is representative of the expression level of the identified mRNA while β is the coefficient of that particular mRNA in multivariate cox regression analysis. The median score was employed to categorize the patients into low-risk (< median number) or high-risk (≥ median number) groups. The signature was corroborated by the following tests: Kaplan-Meier analysis employing the “survival” and “survminer” R package was done for inter-group overall survival. The “survival ROC” R package was employed to score the accuracy of prediction by the developed signature using time- dependent ROC curve analysis. To determine the feasibility and reliability of the proposed prognostic signature, all dataset of CESC patients (N = 350) was randomly divided into train datasets(N = 175) and test (N = 175) using the ‘Caret’ package. Subsequently, the same statistical methods were used for validating both test datasets.

Immune Cell Infiltration, Tumor Microenvironment

In order to evaluate the proportion of tumor immune infiltration cells(TIICs) in TME, The correlation between risk score and 22 immune infiltrating cells was analyzed by Cibersort[25]. We explored the associations between the fractions of 22 infiltrating immune cells and five genes in the PRG_score. Furthermore, the correlations between PRG_score and Immune scores, Stromal scores and ESTIMATE scores were also included[26]. We also used box plot to show the differential expression level of immune checkpoints between high-risk and low-risk groups. Furthermore, we analyzed the relationships between the two risk groups and MSI.

Drug Susceptibility Analysis

In order to explore the difference in the efficacy of chemotherapeutic drugs between the two groups, we calculated the half inhibitory concentration (IC50) of common chemotherapeutic drugs for CESC using the pRRophetic software package.

The predictive nomogram

A hybrid nomogram[27] was built with the “rms” R package encompassing the developed mRNA signature and clinicopathological attributes to predict CESC patient OS (1 year, 2 years, and 3 years). Calibration curve was used to check the predictive power of nomogram.

Statistical analysis

Spearman was used to calculate the correlation coefficient of PRG_score and immune cell. Wilcoxon test was used to compare the differences. Receiver operating characteristic (ROC) curve was used to verify the validity of the model. Survminer package was used to determine the cutoff point of survival information for each dataset. The “surv-cutpoint” function was used to dichotomy PRG_score, and all potential cutting points were repeatedly tested to find the maximum rank statistic, and then the patients were divided into the PRG_score-high group and the PRG_score-low group according to the maximum selected log-rank statistics, so as to reduce the calculated batch effect. Survival curves for prognostic analysis were generated using the Kaplan-Meier method, and the log-rank test was used to determine the significance of the differences. Univariate Cox regression model was used to calculate the hazard ratio (HR) between genes and CESC patients. All statistical analysis was two-side and considered p < 0.05 as statistical significance.

Genetic variation profile of PRGs in CESC

We summarised the incidence of copy number variations and somatic mutations of 52 PRGs(Table S2) in CESC.As shown in Fig. 1A, 82 of 289 (28.37%) CESC samples demonstrated genetic mutations. Of these,TP53 showed the highest frequency of mutations. In Fig. 1B, we analyzed the correlation between 52 PRGS in CESC. Figure 1C p r e s e n t s the location of CNV alterations of these 52 PRGs on chromosomes. We also investigated CNV alteration frequency, which revealed that these 52 PRGs showed prevalent CNV alterations in Fig. 1D. More than half of the 52 PRGs had copy number amplification, while the CNV deletion frequencies of CASP9, CASP8.BAX.ELANE, GPX4, CASP4.CASP5.CASP1.GZMA.IRF2.CHMP2B.HMGB1.NOD2.IL18, NLRP1, CHMP7.TIRAP and TP53 were widespread (Fig. 1D).

Functional enrichment analysis of the Differentially Expressed Pyroptosis Genes of PRGs and the prognostic value of PRGs

We then screened 26 PRGs that were differentially expressed in normal tissues and CESC from 52 PRGS. More specifically, the expression of BAK1.BAX.CASP3.CASP4.CASP5.CHMP4C.CYCS.GZMB.IL18.IL1A.TP53.TP63.AIM2.CASP6.CASP8.GSDMB.GSDMC.NLRP2.NLRP7.NOD2.PYCARD and TNF was increased, while the expression of CHMP6,ELANE,NLRP1 and NOD1 was decreased in CESC compared with normal tissues (Fig. 2A). To further explore the interactions of these Pyroptosis Genes, we conducted a protein–protein interaction (PPI) analysis, and the results are shown in Fig. 2B. To clarify the function of DEGs, the pathways were analysed using GO and KEGG databases. We found that these 26 PRGs were mainly involved in the positive regulation of cysteine-type endopeptidase activity involved in apoptotic process, interleukin-1 production, regulation of inflammatory response, pyroptosis, inflammasome complex, and cysteine-type endopeptidase activity in GO analysis (Fig. 2C). Moreover, KEGG pathway analysis suggested that 26 PRGs were mainly involved in the NOD-like receptor signalling pathway, salmonelia infection, apoptosis, influenza A, pathogenic Escherichia coli infection, and lipid and atherosclerosis (Fig. 2D). The comprehensive landscape of PRG interactions, regulator connections, and their prognostic value in patients with CESC patients was demonstrated in a pyroptosis network (Fig. 2E and Table S4).The results of univariate Cox regression and Kaplan–Meier analysis revealed the prognostic values of 34 PRGs in patients with CESC (Table S3 and Fig. 2F), and P < .05 was selected as the threshold for filtering.

Molecular Subtypes Mediated by 34 Pyroptosis Regulators

To explore the connections between the expression of the 34 PRGs(PRGs that were not expressed in all samples were excluded) and CESC subtypes, we performed a consensus clustering analysis with all 365 CESC patients. By increasing the clustering variable (k) from 2 to 10, we found that when k = 3, the intragroup correlations were the highest and the inter-group correlations were low, indicating that the 365 CESC patients could be well divided into 3 clusters based on the 34 PRGs (Fig. 3A).The overall survival (OS) time which combined TCGA and GEO database was also compared between the 3 clusters, which display is significant (P = 0.051, Fig. 3B). The gene expression profile and the clinical features including the degree of tumor differentiation (G1-G3), age (≤ 65 or > 65 years) and project(TCGA or GSE30759) are presented in a heatmap, and we found there’s little differences in clinical features between the 3 clusters but it is different between TCGA cohort and GSE30759 cohort(Fig. 3C).

GSVA, ssGSEA, PCA Analyse in Distinct Subtypes

GSVA enrichment analysis showed that subtype A was significantly enriched in immune fully-activated pathways, including Chemokine signaling pathway, T cell receptor signaling pathway, B cell receptor signaling pathway, JAK-STAT signaling pathway, Leishmania infection, Natural killer cell mediated cytotoxicity, Toll like receptor signaling pathway, Cytosolic DNA sensing pathway, NOD like receptor signaling pathway, Antigen processing and Presentation, autoimmune thyroid disease, Graft versus host disease, Allograft rejection, type I diabetes mellitus, Hematopoietic cell lineage, Cviral myocarditis, Primary immunodeficiency, Systemic lupus erythematosus, and Asthma compared with B. Besides, subtype A was significantly enriched in immune fully-activated pathways, including Fatty acid metabolism, Pathogenic escherichia coli infection, small cell lung cancer, Regulation of actin cytoskeleton, Cytosolic DNA sensing pathway, Rig I like receptor signaling pathway, Apoptosis, NOD like receptor signaling pathway, JAK-STAT signaling pathway, Toll like receptor signaling pathway, Leishmania infection, Natural killer cell mediated cytotoxicity, Viral myocarditis, Antigen processing and presentation, Graft versus host disease, Type I diabetes mellitus, Prion diseases, Chemokine signaling pathway, Cytokine receptor interaction and Hematopoietic cell lineage compared with B.(Fig. 4A-C; Table S5).

To investigate the role of PRGs in the TME of CESC, we assessed the correlations between the three subtypes and 23 human immune cell subsets of every CESC sample(Table S6). We observed significant differences in the infiltration of most immune cells between the three subtypes by ssGSEA [28]. The infiltration levels of Activated B cells, Activated CD8 T cells, Activated dendritic cells, Immature B cells, MDSCs and T follicular helper cells were obviously higher in the subtype A, while CD56 bright natural killer cells, Gamma delta T cells, Macrophagena, Natural killer T cells, Natural killer cells, Neutrophilna, Regulatory T cells and Type 2 T helper cells had significantly higher infiltration in subtype C. All immune cells were obviously highest in three subtypes(Fig. 4D).PCA analysis revealed significant differences in the pyroptosis expression profiles between the three subtypes (Fig. 4E).

Identification of gene subtypes based on DEGs

To explore the potential biological behavior of each pyroptosis pattern, we identified 126 pyroptosis subtype-related DEGs using the R package “limma,” and performed functional enrichment(GO and KEGG) analysis (Fig. 5A-C; Table S7). These pyroptosis subtype-related genes were significantly enriched in biological processes that were correlated with immunity including Extracellular matrix, ERK1 and ERK2 cascade, Chemotaxis of cytokines, Inflammatory response, etc.(Fig. 5B). KEGG analysis indicated enrichment of immune and cancer-related pathways such as Cytokine − cytokine receptor interaction, NOD − like receptor signaling pathway, Chemokine signaling pathway and IL − 17 signaling pathway (Fig. 5C),suggesting that pyroptosis plays a vital role in the immune regulation of the TME. We then conducted univariate Cox regression analysis to identify the prognostic value of 126 subtype-related genes and screened out 32 genes related to RFS time (p < .05), which were used in the subsequent analysis (Table S8).

To further validate this regulation mechanism, a consensus clustering algorithm was used to divide patients into four genomic subtypes based on prognostic genes; namely, gene subtypes A–D. Kaplan-Meier curves showed that patients with gene subtype B had the worst RFS, whereas patients in gene cluster A showed a favorable RFS (log-rank test, P = 0.003; Fig. 5D). In addition, we found there’s distinct differences between this four clusters(Fig. 5E). The four pyroptosis gene subtypes showed significant differences in PRG expression, consistent with the expected results of the pyroptosis patterns (Fig. 5F).In addition, we investigated the associations between immune checkpoints and this four subtypes. Figure 5G shows that 37 immune checkpoints were differentially expressed in the four groups, including PD-1(PDCD1), PD-L1(CD274), and CTLA-4.

Construction the prognostic PRG_score

The PRG_score was established based on the subtype-related DEGs. LASSO and multivariate Cox analyses for 32 pyroptosis subtype-related prognostic DEGs were performed to further select optimum prognostic signature. Followed by LASSO regression analysis, 8 RFS-associated genes remained according to the minimum partial likelihood deviance (Fig. 6A-B). We then performed multivariate Cox regression analysis on 8 RFS-associated genes based on the Akaike information criterion (AIC) value to finally obtain five (CD7,FAM3B,CCL19,GCNT2,CXCL3).According to the results of the multivariate Cox regression analysis, the PRG_score was constructed as follows:

Risk score = (-0.3035* expression of CD7) + (-0.2267*expression of FAM3B) + (-0.1983* expression of CCL19) + (-0.2817* expression of GCNT2) + (0.3199* expression of CXCL3).

Patients with a PRG_score lower than the median risk score were categorized into the lowrisk group (n = 167), whereas those with a PRG_score greater than the median risk score were placed in the high-risk group (n = 183) (Table S9).We observed a significant difference in PRG_score between pyroptosis gene subtypes. The PRG_score of subtype A was the lowest, while that of subtype B was the highest, indicating that a low PRG_score may be closely related to immune activation-related features, while a high PRG_score may be related to stromal activation-related features (Fig. 6C). More importantly, except for the A-C group, there were significant differences between every other two groups. The distributions of risk scores in the three subtypes are shown in Fig. 6C. And the distributions of risk scores in the four pyroptosis subtype-related subtypes are shown in Fig. 6D. Figure 6E showed differentially expressed of PRGs between the high-risk and low-risk groups. Figure 6F illustrates the distribution of patients' survival in the three pyroptosis subtypes, four gene subtypes, and two PRG_score groups. Each patient’s risk status, survival outcome and prognostic genes expression levels were demonstrated in Fig. 7A.D.G. The Kaplan–Meier survival curves revealed that patients with low scores had a significantly favorable overall survival compared to that in patients with high scores (log-rank test, P = 0.023, Fig. 7G). As Shown in Fig. 7M, the area under ROC curve (AUC) was 0.790 at 1 year,0.649 at 3 years and 0.626 at 5 years.

To validate the prognostic performance of the PRG_score, we used the “caret package” in R to randomly classify the patients into training (n = 175) and testing (n = 175) groups at a ratio of 1:1.We calculated the risk score of the testing cohort according to the formula obtained in the training cohort, and used the median value of the risk score of the training cohort as the cutoff value (Table S10). The training cohort was segregated into high (N = 87) and low(N = 88) risk groups. Each patient’s risk status, survival outcome and prognostic genes expression levels were demonstrated in Fig. 7B,E,H, patients had a higher probability of mortality than patients in the high-risk group vs. the low-risk group. The OS for the high-risk group was evidently worse than the low-risk group as demonstrated by the Kaplan-Meier analysis (Fig. 7K,p < .001). As Shown in Fig. 7N, the area under ROC curve (AUC) was 0.790 at 1 year, 0.649 at 2 years and 0.626 at 3 years. Similarly, the validation set was segregated into high (N = 96) and low(N = 79) risk groups. Each patient’s survival outcome, risk status and prognostic genes expression levels were demonstrated in Fig. 7C,F,I, patients had a higher probability of mortality than patients in the high-risk group vs. the low-risk group. The OS for the high-risk group was evidently worse than the low-risk group as demonstrated by the Kaplan-Meier analysis (Fig. 7L,p < .001). As Shown in Fig. 7O, the area under ROC curve (AUC) was 0.817 at 1 year, 0.707 at 2 years and 0.699 at 3 years.

Evaluation of TME and Immune Checkpoints between the high- and low-risk groups

A low PRG_score was not only closely associated with a high immune score but also associated with a high stromal score (Fig. 8A).This suggested that the low risk group has high stromal and immune cell content and low tumor cell purity, which explained why the OS of the low risk group is higher than that of the high risk group, We performed the CIBERSORT algorithm to assess the association between PRG_score and the abundance of immune cells. As shown in the scatter diagrams, the PRG_score was positively correlated with M0 macrophages, Mast cells activated, Neutrophils, T cells CD4 naive and negatively correlated with B cells memory, B cells naive, Mast cells resting, and Monocytes, T cells CD8, T cells follicular helper and T cells regulatory (Tregs)(Fig. 8C-M). We also assessed the relationship between the five genes in the proposed model and the abundance of immune cells. We observed that most immune cells were significantly correlated with the seven genes (Fig. 8B). In addition, we investigated the associations between immune checkpoints and our risk model. Figure 8N shows that 37 immune checkpoints were differentially expressed in the two groups, including PD-1, PD-L1, and CTLA-4.

Mutation and drug susceptibility analysis

More and more evidence shows that patients with elevated TMB may benefit from immunotherapy because of the large number of neoantigens[29].We then analyzed the distribution variations of the somatic mutations between two PRG_score groups among al CESC patients. The top twenty mutated genes in the high- and low-risk groups were TPGS2,ID1,KLHL41,MUM1,C17orf62,ALCAM,PCDH20,PPIL2,NEO1,F2RL3,THOC3,COL18A1,ZNF513,HIPK1,AMPH,DSCR3,DUOX2,PVT1,TREML1 and PRAMEF8(Fig. 9B-C). Figure 9D-E showed the relationship between TMB and the PRG scores. In addition, we synthesized the PRG_score and CSC index values to assess the potential correlation between the PRG_score and CSC in CESC. Figure 9F shows the results of the linear correlation between PRG_score and CSC index. We next selected chemotherapy drugs currently used for the treatment of CESC to evaluate the sensitivities of patients in the low- and high-risk groups to these drugs. Interestingly, we found that the patients in the high PRG_score group had lower IC50 value for Docetaxel and Imatinib, while IC50 values of chemotherapeutics such as Veliparib, Lenalidomide, Methotrexate, Nilotinib, Rapamycin, Sunitinib and Temsirolimus were significantly lower in the patients with low PRG_score. Together, these results showed that PRGs were related to drug sensitivity (Fig. 9A).

Construction and Validation of a Prognostic Nomogram

Considered the inconvenience clinical utility of PRG_score in predicting RFS inpatients with CESC, a nomogram incorporating the PRG_score and clinicopathological parameters was established to predict the 1-, 3- and 5-year OS rates(Fig. 10A). The predictors included PRG_score ,patient grade and age. Individual patients were given total point values, generated by combining the points corresponding to each prognostic variable. High overall point score indicated poor clinical outcome. The calibration curve for the nomogram indicated good consistency between predicted survival and actual observations (Fig. 10B). All sensitive drugs are in supplementary documents.

Cervical cancer is the second most common gynecologic malignancy in the world and the fourth leading cause of cancer death in women[30]. Most cases are diagnosed when the patient is already in an advanced stage. The standard of care for patients with cervical cancer is concurrent radiotherapy plus brachytherapy, but the prognosis is poor[31]. Therefore, there is an urgent need for potential biomarkers to assess the risk of patients with metastatic cervical cancer and the corresponding targeted agents to improve clinical outcomes.

Pyroptosis is a type of programmed cell death that plays a dual role in the occurrence and development of tumors in recent years. On the one hand, normal cells are stimulated by a large number of inflame matory factors released by pyroptosis, leading to their transformation into tumor cells[19]. On the other hand,the promotion of tumor cell pyroptosis could be a new therapeutic target[18]. It can regulate cell morphology, proliferation, infiltration, migration, and chemotherapy resistance through a variety of cell signaling pathways, thereby influencing tumor progression, and may be associated with patient prognosis[32]. However, the regulatory mechanisms and networks that link pyroptosis and CESC are not fully understood. Our study aims to elucidate this aspect and paths that can be used for developing interventions.

The results of the present study revealed global alterations in 52 PRGs at the transcriptional and genetic levels in CESC. In this study, we first studied the mRNA levels of 52 currently known pyroptosis-related genes in CESC and normal tissues and found that most of them were differentially expressed. In terms of gene expression, we first found that genes BAK1.BAX.CASP3.CASP4.CASP5.CHMP4C.CYCS.GZMB.IL18.IL1A.TP53.TP63.AIM2.CASP6.CASP8.GSDMB.GSDMC.NLRP2.NLRP7.NOD2.PYCARD and TNF were chiefly up-regulated in CESC tissues, while the genes CHMP6,ELANE,NLRP1 and NOD1 were down-regulated. Among these, those CESC patients who had low AIM2,BAX,CHMP4A,GZMA,GPX4,CHMP7,CHMP4B,GZMB,IRF2,PRKACA,PYCARD,TP53 expression and high TNF,IL1A,IL1B,CHMP4C,BAK1,CHMP2B expression had poorer survival rates. This is consistent with earlier findings that indicate that these PRGs are closely associated with tumor development. We identified three distinct molecular subtypes based on 34 PRGs by TCGA database and GSE30759 and there is a significant difference in survival between these three groups. The characteristics of the TME also differed significantly between the three subtypes. The CESC subtypes were also characterized by a significant immune activation, including antigen processing and presentation, T and B cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, JAK-STAT signaling pathways, NOD-like, and the Toll-like receptor signaling pathways. Furthermore, differences in mRNA transcriptomes between different pyroptosis subtypes were significantly related to PRG and immune-related biological pathways. We identified four gene subtypes based on the DEGs between the three pyroptosis subtypes. Thus, our findings indicate that PRGs might serve as a predictor for evaluating the clinical outcome and immunotherapy response of CESC. Therefore, we constructed the robust and effective prognostic PRG_score(include 5 genes, CD7,FAM3B,CCL19,GCNT2,CXCL3) and demonstrated its predictive ability. Using univariate Cox proportional risk regression and Lasso regression, patients in the low-risk group had a better survival prognosis. This feature was effective in predicting prognosis in patients with CESC (AUC > 0.7), and several independent prognostic analyses suggested that this feature could be considered as an independent prognostic factor. The associated Nomogram have been shown to be stable and accurate and can therefore be used satisfactorily to predict survival at 1, 3, and 5 years in patients with CESC, providing an additional option for prognostic prediction in CESC. The pyroptosis patterns characterized by immune activation and inhibition showed lower and higher PRG_score, respectively. Patients with low- and high-risk PRG_score showed significantly different prognosis, mutation, TME, immune check-points, and drug susceptibility, but no more meaningful differences were shown in TBM and CSC. Finally, by integrating the PRG_score, tumor Grade and age, we established a quantitative nomogram, which further improved the performance and facilitated the use of the PRG_score. The prognostic model can be used for prognosis stratification of patients with CESC, will assist in better understanding the molecular mechanism of CESC, and will provide new ideas for targeted therapies.

In previous studies, researchers have identified PRGs in ovarian and gastric cancers and constructed and validated several predictive models for PRGs[14, 15]. Therefore, we identified and characterized mRNAs closely associated with PRGs and formulated PRG-related features, including 5 mRNAs(CD7,FAM3B,CCL19,GCNT2,CXCL3) with prognostic value.

In our study, AIM2,BAX, GZMB, PYCARD and TP53 were significantly elevated in CESC, and the survival rate of CESC with high expression was high. While TNF,IL1A and BAK1 were on the contrary. Absent in melanoma 2 (AIM2) was initially identified in melanoma, in which it showed decreased expression. AIM2 consists of a HIN structural domain at the C-terminus and a PYD domain at the N-terminus and can id entify double-strand DNA (dsDNA) of microbes or the host. AIM2 activates CASP-1 through ASC-ediated junctional proteins to promote the maturation and release of IL-1β and IL-18 and to promote pyroptosis[14].Studies have shown that the role of AIM2 in tumorigenesis is bidirectional, with results varying depending on the type of cancer. AIM2 has a tumor suppressor function in colon cancer and liver cancer[19], but a tumor promoting function in skin squamous cell carcinoma (SCC)[33]. Interestingly, in our study, AIM2 seemed to be a cancer-promoting gene, as it was upregulated threefold in tumor tissues; however, it also contributed to prolonged patient survival. In our study, GZMB was significantly elevated in CESC, and the survival rate of CESC with high expression was high. GZMB is a serine protease and a pro-inflammatory molecule that promotes the progression of inflammatory diseases and cancers[34]. It is expressed in uroepithelial carcinoma, pancreatic cancer, and melanoma cells, and is known to promote cancer cell invasion[35–37]. For a long time, GZMB has been established as being influential in protein hydrolysis-mediated apoptosis, while its role in pyroptosis has not been given due attention. GZMB was found to be able to cleave GSDME and kill lymphocytes to activate the process of pyroptosis, thus changing apoptosis to pyroptosis[13]. BAX can form pores on the mitochondrial outer membrane (MOM) after activation, resulting in the release of mitochondrial inter-membrane components, which will cause the activation of caspase cascade[38–40]. Many studies have demonstrated that high concentrations of TNF-α (> 20 ng/mL) can induce the onset of apoptosis. For example, treatment of endothelial progenitor cells with TNF-α 20 ng/mL for 24 h resulted in a significant increase in the number of apoptotic cells[41]. Many studies have shown that PRKACA is crucial to the development and progression of a variety of cancers[42], and also that it can be used as a biomarker for hepatocellular carcinoma[43]. Therefore, these genes have the potential to be new molecular targets for CESC.

Another important finding shows that this prognostic model has a significant relationship with TME, which confirms that focal prolapse plays a significant role in TME. More and more evidences show that cell death caused by blepharoptosis is particularly important to the formation of tumor and TME, and has a more significant impact on tumor immune microenvironment[44]. Studies have confirmed that about 95% of cervical cancer cases may be caused by persistent HPV infection. Therefore, the human immune system plays an important role in the process of human infection response, which lays a foundation for immunotherapy of cancer[45]. This was confirmed in our study, which identified significant differences in TME and immune function between low-risk and high-risk groups. It was found that patients with high TME scores had a better prognosis. This finding is consistent with previous studies that concluded that patients with high TME scores exhibited a stronger anti-tumor immune response, stood to benefit more from immunotherapy, and survive longer[46, 47]. In addition, correlation analysis showed that B cells, mast cells, monocytes and CD8 + T cells were negatively correlated with risk score, while it was positively correlated with neutrophils, CD4 + T cells and macrophages M0. Wang et al. found that when the biological orthogonal shearing system was applied to pyroptosis activated by gasdermin protein, it reshaped the immune microenvironment of tumor, activated a strong anti-tumor immune response mediated by T cells, and played a strong anti-tumor role[48]. It has also been observed that tumor-associated neutrophils are involved in regulating tumor development, while IL-1 β and IL-18 direct neutrophils to tumor sites[49]. At the same time, our study confirmed that the expression of multiple immune checkpoints shows considerable differences between low-risk groups and high-risk groups. Later studies found that the abnormal expression of immune checkpoints will affect the immune microenvironment of tumors and help tumor cells escape the immune response of the body[50, 51]. These studies clearly point out that the immune response and microenvironment of the host are closely related to the progress of CESC.

There is little current research on pyroptosis, especially on its mechanism in CESC. We preliminarily studied the prognostic value of these pyroptosis-related genes and provided theoretical support for future research. However, due to a lack of data, we could not confirm whether these regulators (which have been reported in prior studies) also play corresponding roles in pyroptosis pathways in CESC, and this question deserves further in-depth studies. This study had several limitations. First, all analyses were conducted solely on data from public databases, and all samples used in our study were obtained retrospectively. Therefore, an inherent case selection bias may have influenced the results. Large-scale prospective studies and additional in vivo and in vitro experimental studies are needed to confirm our findings. Furthermore, data on some important clinical variables such as surgery, neoadjuvant chemotherapy, and chemoradiotherapy were unavailable for analysis in most datasets, which may have affected the prognosis of the immune response and pyroptosis state.

Our comprehensive analysis of PRGs revealed an extensive regulatory mechanism by which they affect the tumor- immune-stromal microenvironment, clinicopathological features, prognosis and may be helpful for exploring the functions of PRGs. We also determined the therapeutic liability of PRGs in targeted therapy and immunotherapy in CESC. These findings highlight the crucial clinical implications of PRGs and provide new ideas for guiding personalized immunotherapy strategies for patients with CESC.

Acknowledgements

All authors would like to thank the specimen donors and research groups for the TCGA and GSE30759, which provided data for this article.

Funding

This manuscript is supported by Scientific and Technological Research Foundation of Shaan'xi Province (Juan Ren, 2020JM-368);

Competing Interests

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Author Contributions

All authors contributed to the study conception and design. The first draft of the manuscript was written by Haojing Kang and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Data Availability

Informed consent All the analyzed data used in our research were collected from a public database, such as TCGA and GEO; therefore, informed consent was not required for this analysis.

Ethical approval

Since this was a retrospective study and all the data used in our study were obtained from public databases TCGA and GEO, therefore, ethical approval was not required.

Cohen, P.A., et al., Cervical cancer. Lancet, 2019. 393(10167): p. 169–182.
Ward, Z.J., et al., Global costs, health benefits, and economic benefits of scaling up treatment and imaging modalities for survival of 11 cancers: a simulation-based analysis. Lancet Oncol, 2021. 22(3): p. 341–350.
Sung, H., et al., Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin, 2021. 71(3): p. 209–249.
Crosbie, E.J., et al., Human papillomavirus and cervical cancer. Lancet, 2013. 382(9895): p. 889–99.
Kovacs, S.B. and E.A. Miao, Gasdermins: Effectors of Pyroptosis. Trends Cell Biol, 2017. 27(9): p. 673–684.
Miao, E.A., J.V. Rajan, and A. Aderem, Caspase-1-induced pyroptotic cell death. Immunol Rev, 2011. 243(1): p. 206–14.
Feng, S., D. Fox, and S.M. Man, Mechanisms of Gasdermin Family Members in Inflammasome Signaling and Cell Death. J Mol Biol, 2018. 430(18 Pt B): p. 3068–3080.
Zhang, Y., et al., Plasma membrane changes during programmed cell deaths. Cell Res, 2018. 28(1): p. 9–21.
Frank, D. and J.E. Vince, Pyroptosis versus necroptosis: similarities, differences, and crosstalk. Cell Death Differ, 2019. 26(1): p. 99–114.
Kolb, R., et al., Inflammasomes in cancer: a double-edged sword. Protein Cell, 2014. 5(1): p. 12–20.
Dupaul-Chicoine, J., et al., Control of intestinal homeostasis, colitis, and colitis-associated colorectal cancer by the inflammatory caspases. Immunity, 2010. 32(3): p. 367–78.
Xi, G., et al., GSDMD is required for effector CD8(+) T cell responses to lung cancer cells. Int Immunopharmacol, 2019. 74: p. 105713.
Zhang, Z., et al., Gasdermin E suppresses tumor growth by activating anti-tumor immunity. Nature, 2020. 579(7799): p. 415–420.
Ye, Y., Q. Dai, and H. Qi, A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer. Cell Death Discov, 2021. 7(1): p. 71.
Shao, W., et al., The Pyroptosis-Related Signature Predicts Prognosis and Indicates Immune Microenvironment Infiltration in Gastric Cancer. Front Cell Dev Biol, 2021. 9: p. 676485.
Dong, Z., et al., Identification of a Pyroptosis-Related Gene Signature for Prediction of Overall Survival in Lung Adenocarcinoma. J Oncol, 2021. 2021: p. 6365459.
Zhu, W., et al., A pyroptosis-related lncRNA signature predicts prognosis and immune microenvironment in head and neck squamous cell carcinoma. Int Immunopharmacol, 2021. 101(Pt B): p. 108268.
Xia, X., et al., The role of pyroptosis in cancer: pro-cancer or pro-"host"? Cell Death Dis, 2019. 10(9): p. 650.
Karki, R. and T.D. Kanneganti, Diverging inflammasome signals in tumorigenesis and potential targeting. Nat Rev Cancer, 2019. 19(4): p. 197–214.
Man, S.M. and T.D. Kanneganti, Regulation of inflammasome activation. Immunol Rev, 2015. 265(1): p. 6–21.
Wang, B. and Q. Yin, AIM2 inflammasome activation and regulation: A structural perspective. J Struct Biol, 2017. 200(3): p. 279–282.
Mayakonda, A., et al., Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res, 2018. 28(11): p. 1747–1756.
Wilkerson, M.D. and D.N. Hayes, ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics, 2010. 26(12): p. 1572–3.
Wang, H., et al., Precision Lasso: accounting for correlations and linear dependencies in high-dimensional genomic data. Bioinformatics, 2019. 35(7): p. 1181–1187.
Newman, A.M., et al., Robust enumeration of cell subsets from tissue expression profiles. Nat Methods, 2015. 12(5): p. 453–7.
Yoshihara, K., et al., Inferring tumor purity and stromal and immune cell admixture from expression data. Nat Commun, 2013. 4: p. 2612.
Balachandran, V.P., et al., Nomograms in oncology: more than meets the eye. Lancet Oncol, 2015. 16(4): p. e173-80.
Rooney, M.S., et al., Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell, 2015. 160(1–2): p. 48–61.
Snyder, A., et al., Genetic basis for clinical response to CTLA-4 blockade in melanoma. N Engl J Med, 2014. 371(23): p. 2189–2199.
Flanagan, M.B., Primary High-Risk Human Papillomavirus Testing for Cervical Cancer Screening in the United States: Is It Time? Arch Pathol Lab Med, 2018. 142(6): p. 688–692.
Naga Ch, P., et al., The management of locally advanced cervical cancer. Curr Opin Oncol, 2018. 30(5): p. 323–329.
Zheng, Z. and G. Li, Mechanisms and Therapeutic Regulation of Pyroptosis in Inflammatory Diseases and Cancer. Int J Mol Sci, 2020. 21(4).
Sharma, B.R., R. Karki, and T.D. Kanneganti, Role of AIM2 inflammasome in inflammatory diseases, cancer and infection. Eur J Immunol, 2019. 49(11): p. 1998–2011.
Velotti, F., et al., Granzyme B in Inflammatory Diseases: Apoptosis, Inflammation, Extracellular Matrix Remodeling, Epithelial-to-Mesenchymal Transition and Fibrosis. Front Immunol, 2020. 11: p. 587581.
D'Eliseo, D., et al., Granzyme B is expressed in urothelial carcinoma and promotes cancer cell invasion. Int J Cancer, 2010. 127(6): p. 1283–94.
D'Eliseo, D., et al., Epitelial-to-mesenchimal transition and invasion are upmodulated by tumor-expressed granzyme B and inhibited by docosahexaenoic acid in human colorectal cancer cells. J Exp Clin Cancer Res, 2016. 35: p. 24.
Dufait, I., et al., Perforin and Granzyme B Expressed by Murine Myeloid-Derived Suppressor Cells: A Study on Their Role in Outgrowth of Cancer Cells. Cancers (Basel), 2019. 11(6).
Degterev, A. and J. Yuan, Expansion and evolution of cell death programmes. Nat Rev Mol Cell Biol, 2008. 9(5): p. 378–90.
Delbridge, A.R., et al., Thirty years of BCL-2: translating cell death discoveries into novel cancer therapies. Nat Rev Cancer, 2016. 16(2): p. 99–109.
Czabotar, P.E., et al., Control of apoptosis by the BCL-2 protein family: implications for physiology and therapy. Nat Rev Mol Cell Biol, 2014. 15(1): p. 49–63.
Chen, I.Y. and T. Ichinohe, Response of host inflammasomes to viral infection. Trends Microbiol, 2015. 23(1): p. 55–63.
Stratakis, C.A., Cyclic AMP-dependent protein kinase catalytic subunit A (PRKACA): the expected, the unexpected, and what might be next. J Pathol, 2018. 244(3): p. 257–259.
Kastenhuber, E.R., et al., DNAJB1-PRKACA fusion kinase interacts with β-catenin and the liver regenerative response to drive fibrolamellar hepatocellular carcinoma. Proc Natl Acad Sci U S A, 2017. 114(50): p. 13076–13084.
Yu, P., et al., Pyroptosis: mechanisms and diseases. Signal Transduct Target Ther, 2021. 6(1): p. 128.
Kwon, M.J., et al., Clinical implication of programmed cell death-1 ligand-1 expression in tonsillar squamous cell carcinoma in association with intratumoral heterogeneity, human papillomavirus, and epithelial-to-mesenchymal transition. Hum Pathol, 2018. 80: p. 28–39.
Pasini, F.S., et al., A gene expression profile related to immune dampening in the tumor microenvironment is associated with poor prognosis in gastric adenocarcinoma. J Gastroenterol, 2014. 49(11): p. 1453–66.
Zeng, D., et al., Gene expression profiles for a prognostic immunoscore in gastric cancer. Br J Surg, 2018. 105(10): p. 1338–1348.
Wang, Q., et al., A bioorthogonal system reveals antitumor immune function of pyroptosis. Nature, 2020. 579(7799): p. 421–426.
Jorgensen, I., et al., IL-1β, IL-18, and eicosanoids promote neutrophil recruitment to pore-induced intracellular traps following pyroptosis. Eur J Immunol, 2016. 46(12): p. 2761–2766.
Saleh, R., et al., Role of Epigenetic Modifications in Inhibitory Immune Checkpoints in Cancer Development and Progression. Front Immunol, 2020. 11: p. 1469.
Rotte, A., J.Y. Jin, and V. Lemaire, Mechanistic overview of immune checkpoints to support the rational design of their combinations in cancer immunotherapy. Ann Oncol, 2018. 29(1): p. 71–83.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of Pyroptosis Genes in cervical cancer and construction of a prognostic model for pyroptosis-associated mRNA

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Sources of Cervix Cancer Datasets and Preprocessing

Mutation Analysis of PRGs and PRG_score

Identification of Differentially Expressed Pyroptosis Genes and Functional Analysis

Analysis of Pyroptosis Genes in Prognosis

Tumor Classification Based on the PRGs

Identification of Differentially Expressed Genes Among Clusters

Consensus Clustering of DEGs

Development and Validation of the Gene Prognostic Model Based on PRGs

Immune Cell Infiltration, Tumor Microenvironment

Drug Susceptibility Analysis

The predictive nomogram

Statistical analysis

Results

Genetic variation profile of PRGs in CESC

Functional enrichment analysis of the Differentially Expressed Pyroptosis Genes of PRGs and the prognostic value of PRGs

Molecular Subtypes Mediated by 34 Pyroptosis Regulators

GSVA, ssGSEA, PCA Analyse in Distinct Subtypes

Identification of gene subtypes based on DEGs

Construction the prognostic PRG_score

Evaluation of TME and Immune Checkpoints between the high- and low-risk groups

Mutation and drug susceptibility analysis

Construction and Validation of a Prognostic Nomogram

Discussion

Conclusions

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1