Bo Liu

1Tianjin University of Traditional Chinese Medicine,Tianjin,China.

Abstract Colorectal cancer is a common clinical tumor with a poor prognosis.In recent years,it has been demonstrated that pyroptosis is a natural immune process and plays a significant part in tumor advancement,but the correlation between the expression of pyroptosis-related genes and prognosis in the CRC is unclear.Our study identified 40 DEGs associated with pyroptosis in CRC,and based on these DEGs,all CRC patients could be classified into two subtypes.Using Cox regression analysis,we analyzed the prognostic value of all pyroptosis-associated genes in the TCGA database and constructed a 15-gene model and divided all CRC patients in TCGA into high- and low-risk groups,with the low-risk group having a considerably higher survival time than the high-risk group(P ＜0.001).Using the median risk score of the TCGA cohort,we also divided the CRC patients in the(GEO)cohort into two risk subgroups and confirmed that the overall survival(OS)time was significantly longer in the low-risk group than in the high-risk group(p=0.003).Using univariate COX analysis and multifactor COX analysis,we found that risk score was an independent predictor of OS in CRC patients.By GO and KEGG analyses,we found that DEGs associated with pyroptosis functioned mainly through regulation of extracellular functions.In a word,pyroptosis-related genes play a significant part in colorectal cancer and can be utilized to anticipate the prognosis of CRC.

Key Words:colorectal cancer;pyroptosis;TCGA;GEO;gene;prognosis

Introduction

Colorectal cancer(CRC)is a regular malignancy of the digestive system,approximately 41%percent of all colorectal cancers happen in the proximal colon,with approximately 22% percent containing distal colon and 28% percent involving rectum[1].In the United States,CRC is the third most regular cancer diagnosed,which has been an rise in amount of young patients over the last few years[2].The current treatments for CRC are surgery and chemotherapy in the early stages.In spite of progressive treatments,a large number of advanced tumors have bad prognosis[3].Given the limitations of CRC treatment,new therapeutic targets are desired to enhance treatment results in CRC;therefore,there is an extremely important need for credible novel prognostic models to make targeted therapy more feasible.

Pyroptosis is a form of programmed cell death,which is first discovered in myeloid cells infected by pathogens or bacteria in 1992[4].Pyroptosis cells are characterized by cell membrane pore formation,cytoplasmic swelling,membrane rupture and the release of cytosolic contents[5].Gasdermins are a family of proteins which play a key role in the activation of pyroptosis[6].In the presence of different microbial and endogenous stimuli,gasdermins is cleaved by pyroptotic caspases,the N-terminal domain of certain gasdermins diffuse into the lipid components,form pores in the cell membrane and conduct the pyroptosis induction role[7,8].GSDMD and GSDME are two important molecules in the Gasdermins family,which are necessary in the regulation of pyroptosis[9,10].GSDMD is a 53-kDa protein situated downstream of the pyroptotic caspases,which mainly expressed in the gastrointestinal tract and skin.Recent studies have shown that certain drugs or molecules can trigger GSDMD-mediated pyroptosis in various types of cancer,suggesting that pyroptosis is implicated in the pathogenesis of cancer and may be a new target for cancer treatment in the future[11].GSDME-mediated pyroptosis is usually triggered by chemotherapeutic agents or targeted therapeutic agents.Recent studies have shown that certain chemotherapeutic agents can modulate the expression of GSDME and increase the sensitivity of chemotherapeutic agents in cancer treatment[12,13].

Based on the current studys,we learn that pyroptosis plays a critical role in cancer development and treatment.Nevertheless,its particular function in colorectal cancer is little researched,so we managed a comprehensive study to figure out the expression levels of pyroptosis-related genes between normal and cancer samples,analyze the prognostic value of pyroptosis-related genes,evaluate the functions of the genes.

Results

Study of differentially expressed genes in normal and cancer samples

In the Genotype-Tissue Expression(GTEx)and The Cancer Genome Atlas(TCGA)database,the 53 pyroptosis-related gene expression levels were compared from 515 normal and 452 cancer samples,resulting in 40 differentially expressed genes(DEGs)linked to pyroptosis(allP＜ 0.01).As shown in Figure 1A(blue:low expression level;red:high expression level)),23 genes(CHMP2A,CHMP6,NLRP2,NLRP7,GZMA,TP63,NLRP1,NLRC4,NLRP3,TIRAP,ELANE,PRKACA,IL18,CASP5,GSDMB,BAK1,CASP9,CHMP3,CHMP2B,CYCS,CASP3,IRF2,CHMP7)of the 40 DEGs were lowly expressed in the cancer samples,while 17 genes(BAX,GPX4,PJVK,NOD2,NOD1,PLCG1,CHMP4C,HMGB1,CASP8,GZMB,IL1A,IL1B,GSDMC,IL6,CASP4,GSDMA,TP53)were upregulated.

In order to examine the relationship between differentially expressed genes,we imported all DEGs into the String database[14]and constructed a protein–protein interaction(PPI)network,as shown in Figure 1B.In the PPI network,we set the confidence level to 0.9 and identified IL18,IL1B,IL6,CASP8,NLRP3,CASP3,TP53,NLRC4 as the core genes in the network.Based on the PPI network,we mapped the correlation network of all pyroptosis-related genes,which is shown in Figure 1C(red:positive correlations;blue:negative correlations).

Cancer categorization derived from 40 DEGs

In the TCGA database,we carried out a consistent clustering analysis of all 452 CRC patients to analyse the connection between the expression of 40 pyroptosis-related DEGs and CRC subtypes.We increased the clustering variable( k )from 2 to 9,and we found the highest intra-group correlations and lower inter-group correlations when k = 2,implying that 452 CRC patients could be separated into two clusters derived from 40 DEGs(Figure2A).We presented the clinical characteristics of the two clusters using heatmaps,but we found little difference between the two clusters(Figure2B).At the same time,we also compared survival between the two clusters and found a dramatic difference between the two clusters(P= 0.042 ＜0.05)(Figure2C).

Production of prognostic gene modules in the TCGA series

We conducted the survival information of 446 CRC patients using univariate Cox regression analysis,resulting in the identification of 15 genes that were associated with survival information.Of 15 genes,14 genes(TMPRSS11E,MAPK12,EGFL7,TMEM88,CALB2,RNF207,HEYL,HOXC11,P2RX5,LINGO1,UPK3B,MID2,IL20RB,IFITM10)were associated with an increased risk of HR＞1,while CCL22 was a protective gene for HRs ＜1(Figure3A).A 15-gene signature was established by executing a minimum absolute shrinkage and selection operator(LASSO)Cox regression analysis based on the optimal λ(Figure3B).The distribution of risk scores for the 15 genes is as follows:risk score =(0.064＊TMPRSS11E exp)+(0.103＊MAPK12 exp)+(0.023＊EGFL7 exp)+(0.028＊TMEM88 exp)+(0.008＊CALB2 exp)+(0.063＊RNF207 exp)+(0.172＊HEYL exp)+(0.018＊HOXC11 exp)+(0.154＊P2RX5 exp)+(0.036＊ LINGO1 exp)+(0.221＊UPK3B exp)+(0.081＊MID2 exp)+(-0.178＊CCL22 exp)+(0.254＊IL20RB exp.)+(0.159＊IFITM10 exp).The 446 CRC patients were equally divided into low and high risk subgroups according to the median score calculated by the risk score formula(Figure3C).CRC patients with different risks are well divided into two clusters derived from the principal component analysis(PCA)(Figure3D).The results showed that patients in the high risk group had more deaths and shorter survival times(Figure3E)and that there was a significant difference in survival times between the low and high risk groups(P＜0.001)(Figure3F).Applying time-dependent subject operating characteristic(ROC)analysis to assess the sensitivity and specificity of the prognostic model,we observed an area under the ROC curve(AUC)of 0.782 at 1 year,0.766 at 2 years,and a 3-year survival rate of 0.757(Figure3G).

Verification of risk models in the GEO series

We used 204 CRC patients from the GEO sequence(GSE87211)as the validation set to verify the stability and plausibility of the risk models.Derived from the median score calculated by the risk score formula,203 CRC patients were equally splited into low and high risk groups(Figure4A).PCA also displayed a clear separation between the two groups(Figure4B).The low-risk group had a longer survival time and lower mortality than the high-risk group(Figure4C),and there was also a remarkable difference in survival rates between the two groups(P=0.003)(Figure4D).The area of the ROC curve is AUC = 0.724 for 1 year,AUC = 0.726 for 2 years and 0.700 for 3 years,which shows that our model has a good predictive effect(Figure4E).

Independent prognostic value of the risk model

Using univariate and multivariate Cox regression analyses,we assessed whether risk scores from the autogenic trait model could be used as an independent prognostic factor.Univariate Cox regression analysis showed that risk score was an independent predictor of patient survival in the TCGA(HR = 5.042,95% CI:3.648-6.970)(Figure5A)and GEO series(HR = 3.403,95% CI:1.271-9.113)(Figure5B).Similarly,the results of the multifactorial COX regression analysis revealed that the risk score was an independent predictor of low patient survival in the TCGA(HR = 3.658,95%CI:2.567-5.212)(Figure5C)and GEO series(HR = 3.136,95% CI:1.131-8.694)(Figure5D).Furthermore,we plotted a heat map through the clinical characteristics of the patients in the TCGA sequence and found that the TNM stages of the patients differed significantly among the high and low risk groups(P＜0.05)(Figure5E).

Functional analysis derived from risk models

To analyse the differences in gene function between the high and low risk groups,we used the "limma" R package to identify the genes that differed between the high and low risk groups in the TCGA sequence,setting FDR ＜0.05 and |log2FC | ≥0.585 as conditions,and we eventually obtained 186 DEGs.Based on these 186 DEGs,we performed Gene ontology(GO)enrichment analysis(Figure6A)and Kyoto Encyclopaedia of Genes and Genomes(KEGG)pathway analysis(Figure6B)to analyze the functions of the genes,and the results showed that DEGs are mainly involved in extracellular matrix structuralconstituent,proteindigestionandabsorption,ECM-receptor interaction,PI3K-Akt signaling pathway,Human papillomavirus infection and so on.

Figure 1 Expression and interactions of 40 DEGs

Figure 2 Cancer categorization derived from 40 DEGs

Figure 3 Construction of risk characteristics in the TCGA series

Figure 4.GEO validation of risk signatures

Figure 5.Univariate and multivariate Cox regression analyses for the risk score

Figure 6.Functional analysis derived from DEGs of the risk models in the TCGA sequence

Dicussion

In our study,we identified 40 differential genes associated with pyroptosis from CRC samples and normal samples in the TCGA database.Clustering analysis of the 40 DEGs revealed significant differences in survival time between the two clusters,but not in clinical characteristics.To explore the prognostic value of these DEGs,we used Cox univariate analysis and LASSO Cox regression analysis to construct a 15-gene signature by optimal λ values,which was well validated in the TCGA and GEO datasets.Finally,we performed a functional analysis of DEGs in the highand low risk groups,and the results showed that DEGs are mainly involved in extracellular matrix structuralconstituent,proteindigestionandabsorption,ECM-receptor interaction,PI3K-Akt signaling pathway,Human papillomavirus infection and so on.

Pyroptosis,a lytic,inflammatory type of regulated cell death that requires membrane-damaging gasdermin proteins,characterized by the swelling and lysis of cells,and release of many proinflammatory factors[15].Recent studies have shown that pyroptosis is closely linked to various human diseases,especially tumours.On the one hand,pyroptosis can prevent tumour development;on the other hand,as a form of pro-inflammatory death,pyroptosis can promote tumour growth by forming a microenvironment suitable for tumour cell growth.Moreover,the induction of tumour pyroptosis is also considered a promising cancer treatment strategy.In CRC,the role and prognostic value of pyroptosis -related genes are still poorly understood.Thus we constructed a 15-gene(TMPRSS11E,MAPK12,EGFL7,TMEM88,CALB2,RNF207,HEYL,HOXC11,P2RX5,LINGO1,UPK3B,MID2,IL20RB,IFITM10,CCL22)signature to predict survival time for CRC patients.Transmembrane protease serine 11E(TMPRSS11E),belongs to the type II transmembrane serine protease(TTSP)family.TMEM185A has been reported to play a modulating role in the development of ovarian,bladder and oesophageal cancers[16-18].It has been demonstrated that TMPRSS11E can downregulate the EGFR/AKT signaling pathway[17],which in turn is positively connected with the expression of GSDMD[19],thus it is conjectured that TMPRSS11E can have an cause on pyroptosis by regulating the EGFR/AKT signaling pathway,which in turn regulates the expression of GSDMD.In our study,we found that MAPK12 is an oncogenic gene that is highly expressed in high-risk groups.MAPK12 can affect the risk and survival of colorectal cancer patients by modulating the MAPK signaling pathway[20],MAPK can mediate the signalling of pyroptosis,which in turn affects cell proliferation and differentiation[21].Epidermal growth factor-like domain-containing protein 7(EGFL7)is overexpressed in colorectal cancer,and it acts as an oncogene,regulated CRC invasion and anoikis through PI3K/AKT signaling[22],our study also shows that EGFL7 is highly expressed in colon cancer samples,but further studies are needed to investigate its role in pyroptosis.TMEM88[23],CALB2[24],HEYL[25],HOXC11[26],P2RX5[27],IL20RB[28]have been demonstrated to be biomarkers for the diagnosis and prognosis of colorectal cancer,as verified in our study,but the regulation of pyroptosis by these genes has not been obviously reported.We found that five genes,RNF207,LINGO1,UPK3B,MID2,IFITM10,were highly expressed in the tumour samples,suggesting that it plays a role as a pro-oncogene in this study.Further studies may concentrate on whether/how these genes are related to pyroptosis and tumour promotion.Our study confirms that CCL22 is an oncogene in colorectal cancer and is lowly expressed in tumour samples.CCL22 is a chemokine and has been demonstrated to increase survival in CRC patients by recruiting beneficial T cells to CRC tumour tissue[29],but further study is demanded on the regulation of pyroptosis by CCL22.

Pyroptosis is a form of programmed cell death,which is characteristiced by cell membrane pore formation,cytoplasmic swelling,membrane rupture and the release of cytosolic contents such as IL-1β but into the extracellular environment,amplifying the local or systemic inflammatory effects[5,30].Our studies have also demonstrated the involvement of pyroptosis-related DEGs in extracellular matrix structural constituents,protein digestion and absorption,ECM-receptor interaction,PI3K-Akt signaling pathway,Human papillomavirus infection and other biological functions.Based on the results we obtained,we can hypothesize that DEGs in the high and low risk groups induce pyroptosis by regulating extracellular functions.

There is little current research on pyroptosis,especially on its mechanism in CRC.Our study identified a 15 pyroptosis-related gene model to predict OS for CRC patients.However,due to the lack of adequateresearch,we cannot determine the role of these 15 genes related to pyroptosis in CRC,and this question deserves further in-depth studies.

In conclusion,our study confirmed that pyroptosis is associated with CRC,as most pyroptosis-related genes are expressed differently between normal and CRC samples.Also,the scores produced based on our risk profile of 15 pyroptosis-related genes were independent risk factors for predicting OS in the TCGA and GEO cohorts.The function of DEGs between the low- and high-risk groups was related to extracellular regulation.Our research offers a novel genetic signature for predicting prognosis in CRC patients.

Materials and methods

Datasets

We downloaded RNA sequencing(RNA-seq)data and clinical information from the TCGA database for 514 normal samples and 452 CRC samples.At the same time,we downloaded the corresponding RNA sequencing(RNA-seq)data and the corresponding clinical features from the GEO database(https://www.ncbi.nlm.nih.gov/geo/,ID:GSE87211)as external validation sequences.

Confirmation of differentially expressed pyroptosis-related genes

Based on 514 normal and 452 CRC samples from the TCGA database,we obtained 40 DEGs using the "limma" package and setting a p-value of for ＜0.05,and constructed a PPI network of these DEGs using the String database(https://www.string-db.org/).

Construction and validation of a prognostic model for pyroptosis-related genes

We used COX regression analysis to assess the correlation between genes and survival time in the TCGA database and identified 15 genes that were associated with survival.We then developed a prognostic model using LASSO Cox regression and determined the risk scores for the 15 genes using the 'scale' function in R software.All CRC patients in the TCGA database were categorised into low and high risk subgroups based on median risk scores.Using "survival","survminer"and "timeROC" in the R package,we compared the OS times between the two groups and 3-year ROC curves were plotted.PCA derived from the 15-gene signature was also performed by the“prcomp”function in the“stats”R package.Then we validated this using relevant cohort from the GEO database(GSE87211),applying median risk scores from the TCGA cohort,and patients in the GSE87211 cohort were also divided from low- or high-risk subgroups,which were then compared to validate the genetic model.

Standalone prognostic assessment of risk scores

We obtained clinical information on CRC patients from the TCGA database and the GEO database and analysed this clinical information using univariate and multivariate Cox regression models.

Functional analysis of DEGs in high and low risk groups

Based on specific criteria(|log 2 FC| ≥0.585 and FDR ＜0.05),we filtered out DEGs between the low- and high-risk groups and used the"clusterProfiler" package to perform GO and KEGG analyses on these DEGs.