The p53 Pathway Related Genes Predict the Prognosis of Colon Cancer

Background

Colon cancer is one of the most common malignancies around the world. In 2020, more than 570,000 deaths and 1.1 million new cases were confirmed worldwide (http://gco.iarc.fr/). At present, the mainstream treatments for colon cancer include surgical resection, chemotherapy and radiotherapy, immunotherapy and targeted therapy. Because patients with early colon cancer have no symptoms and the risk groups are not taken seriously, most colon cancer patients are already at an advanced stage when they are diagnosed.¹ Therefore, the discovery of more new therapeutic targets is essential for the treatment of colon cancer.

p53 is a tumor suppressor gene. Mutations of p53 gene have been found in more than 50% of human tumor tissues.² It has been proven to play an important role in the process of colon adenoma canceration.³ The protein encoded by the p53 gene is a transcription factor, which is mainly distributed in the cell nucleus and can specifically bind to DNA to control the initiation of the cell cycle.⁴ The normal biological function of p53 is like a “genome guardian”. In the G1 phase, DNA damage is checked and the integrity of the genome is monitored, if there is DNA damage, the p53 protein prevents DNA replication to provide enough time for DNA repair, if the repair fails, p53 Protein triggers apoptosis.⁵ After the mutation of p53 gene, the tumor suppressor gene is transformed into an oncogene. Due to the changes in the spatial conformation of the encoded protein, the regulation of cell growth, apoptosis and DNA repair is lost, leading to the occurrence of cancer. The overexpression of p53 gene is associated with a higher recurrence rate and survival rate of colon cancer patients after surgery has been widely confirmed.^6–8 In addition, the overexpression of p53 gene is also related to the clinicopathological features and chemotherapy efficacy of colon cancer.^9–13 In recent years, the cell signal transduction pathway mediated by p53 gene has been gradually recognized.

Given the important position of p53 gene in the occurrence and prognosis of colon cancer, it is necessary to further clarify the relationship between p53 pathway-related genes and the prognosis of colon cancer. However, no relevant literature has been reported yet. This study aims to clarify the expression levels of p53 pathway-related genes in normal colon tissues and colon cancer tissues, analyze the prognostic value of these genes, and explore the correlation between p53 pathway-related genes and the tumor immune microenvironment.

Methods

Data Acquisition and Preparation

This study used samples clinicopathological and prognostic information from The Cancer Genome Atlas (TCGA) dataset, which was downloaded from an open-access resource website (https://www.cbioportal.org/).¹⁴ After removing invalid values and duplicate values from the data, obtain clinicopathological and prognostic information of 540 samples, including sex, age, stage, overall survival and status, disease-specific survival and status, progression-free survival and status. The mRNA data is downloaded in the TCGA database (https://www.cancer.gov/). After data sorting, 41 normal colon tissue mRNA data and 480 colon cancer tissue mRNA data are obtained. 200 p53 pathway-related genes were acquired from an open-access resource website (http://www.gsea-msigdb.org/gsea/msigdb/search.jsp), contributed by Arthur Liberzon (MSigDB team).

Statistical Analysis

The statistical analysis of this study was completed by R version 3.6.3 (http://www.r-project.org). P <0.05 on both sides is considered statistically significant. Through the “limma” package, the differentially expressed p53 pathway-related genes in normal colon tissues and colon cancer tissues that satisfy |log2FC| ≥ 1 and FDR<0.05 are screened out. A consistent cluster analysis was performed on 540 patients to explore the relationship between differentially expressed genes (DEGs) and colon cancer subtypes. The clustering variable (K) is set to 2~10, and the variable with the lowest correlation between the groups and the highest correlation within the group is selected and the corresponding cluster is obtained. The Kaplan-Meier method is used to compare overall survival (OS), disease-specific survival (DSS) and progression-free survival (PFS) between clusters and draw the corresponding survival curve. Carry out univariate and multivariate Cox regression analysis on DEGs and draw forest plots to obtain DEGs related to survival. The risk score of each patient was obtained by multiplying the expression of DEGs related to the survival of each patient by the corresponding regression coefficient and adding them. According to the median risk score, patients were divided into a low-risk score group and a high-risk score group. Based on survival-related DEGs, use the “prcomp” function and the “Rtsne” package to evaluate whether patients with different risks are well divided into two clusters through principal component analysis (PCA) and t-distributed Stochastic Neighbor Embedding (t-SNE), and used the Kaplan-Meier method to compare the survival differences between the two clusters. Match each patient’s risk score with clinicopathological information, perform univariate and multivariate Cox regression analysis and display it through a forest plot. Through the “limma” package, select DEGs in the high/low risk scoring group, and then uses the “clusterProfiler” package to analyze DEGs’ gene ontology(GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genome(KEGG) pathway analysis.¹⁵ The “gsva” package is used to perform single-sample gene set enrichment analysis (ssGSEA) to calculate the proportion of infiltrating immune cells and evaluate the activity of immune-related pathways.

Results

We identified a total of 28 DEGs, including 6 up-regulated genes (SERPINB5, TNFSF9, DDIT4, CDKN2A, KLK8, TRIB3) and 22 down-regulated genes (TPRKB, FAM162A, HMOX1, SERTAD3, HDAC3, ZBTB16, KLF4, RXRA, TXNIP, BMP2, FAS, FUCA1, MXD1, BAK1, BTG1, NHLH2, CTSD, RETSAT, PTPRE, TGFA, PMM1, EPHX1). We found that when K=2, 480 samples can be well divided into two clusters according to the above DEGs (Figure 1A). There was a significant difference in PFS between the two clusters (P=0.011) (Figure 1B). There was no significant difference in OS and DSS between the two clusters (P=0.180; P=0.120, respectively) (Figure 1C and D), but it was different from cluster 1. In comparison, cluster 2 has better OS and DSS. The correlation between clusters and clinical characteristics including gender (Male or Female), age (<65 or ≥65 years), stage (Stage I–IV) and survival status (Censored or Progression) is shown in the heatmap, except for survival status, we found almost no difference in clinical characteristics between the two clusters (Figure 1E).

Figure 1 Tumor classification based on p53 pathway related genes. (A) According to consistent cluster analysis (k=2), 480 colon patients were divided into two clusters; (B) Kaplan-Meier PFS curve of two clusters; (C) Kaplan-Meier OS curve of two clusters; (D) two clusters The Kaplan-Meier DSS curve; (E) Is a heat map based on two clusters, clinicopathological characteristics and PFS status. Abbreviations: PFS, progression free survival; OS, overall survival; DSS, disease free survival.

By matching the DEGs with the patient’s PFS and status, and through univariate COX regression analysis, we obtained three survival-related DEGs (CDKN2A (HR = 1.0005, 95% CI: 1.0002–1.0009, P=0.0040), BAK1 (HR = 1.0000, 95% CI: 1.0000–1.0000, P=0.0143), BTG1 (HR = 1.0001, 95% CI: 1.0000–1.0001, P=0.0284), Figure 2A), the above survival-related DEGs were subjected to multivariate COX regression analysis, We found that the differential expression of CDKN2A gene has a significant effect on PFS (P=0.0281, Figure 2B). By calculating the risk score of each patient, the patients were divided into low-risk groups and high-risk groups. Both PCA and t-SNE showed that patients with different risks were well divided into two clusters (Figure 2D and E). As the risk score increases, the survival time of patients is relatively reduced, and the number of patients progressing increases (Figure 2C). There is a significant difference in PFS between the low-risk group and the high-risk group (P=0.015, Figure 2F). The univariate Cox regression analysis after matching the clinicopathological information with the risk score showed that the risk score is a significant indicator for predicting PFS of colon cancer ((HR = 1.6195, 95% CI: 1.0940–2.3975, P=0.016), Figure 3A) The multivariate COX regression analysis after adjusting for other confounding factors also showed that risk score is an independent prognostic factor for predicting PFS in colon cancer ((HR = 1.7008, 95% CI: 1.1479–2.5201, P=0.008), Figure 3B). In addition, based on the heatmap of the high/low risk score group and clinicopathological information, we only found that the patient’s survival status was significantly unevenly distributed between the high risk score group and the low risk score group (Figure 3C). Through GO analysis and KEGG enrichment analysis, we believe that DEGs in the high/low risk score group are mainly related to the neutrophil-mediated immunity pathway, chemokine signaling pathway, Rap1 signaling pathway and MAPK signaling pathway (Figure 4A and B).

Figure 2 Construction of risk characteristics based on differentially expressed p53 pathway-related genes. (A and B) Univariate and multivariate Cox regression analysis based on differentially expressed p53 pathway-related genes; (C) PFS status of each patient (low-risk population: the left side of the dotted line; high-risk population: the right side of the dotted line); (D and E) Based on each PCA chart and tSNE chart of patient risk score; (F) Kaplan-Meier curve of PFS of high-risk group and low-risk group. Abbreviations: PFS, progression free survival; PCA, principal component analysis; tSNE, t-distributed Stochastic Neighbor Embedding.

Figure 3 Univariate and multivariate Cox regression analysis of risk score. (A and B) Univariate and multivariate Cox regression analysis of risk score; (C) Heatmap based on risk score group, clinicopathological characteristics and PFS status.

Figure 4 Functional analysis of DEGs between two risk groups. (A) GO enriched bubble chart (the larger the bubble, the more enriched genes, the darker the red, the more obvious the difference); (B) The bar chart of the KEGG pathway (the longer the bar indicates the more enriched genes, and the darker the red Indicates that the difference is more obvious). Abbreviations: DEGs, differentially expressed genes; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genome.

Based on functional analysis, we further compared the enrichment scores of 16 immune cells and the activity analysis of 13 immune-related pathways between the high and low risk score groups through ssGSEA. Immune cells with lower levels of infiltration in the high-risk group include activate dendritic cells (aDCs), B-cells, dendritic cells (DCs), immature dendritic cells (iDCs), Macrophages, Mast-cells, plasmacytoid dendritic cells (pDCs), and T helper (Th) cells (Tfh, Th1), tumor-infiltrating lymphocytes (TIL), and higher levels of infiltrating immune cells include Neutrophils, natural killer (NK) cells, T helper cells, and regulatory T (Treg) cells (Figure 5A). Except for human leukocyte antigen (HLA) pathway, major histocompatibility complex class 1 (MHC-class-1) pathway, parainflammation pathway, type-1 Interferon lambdas (IFN) response pathway, the activity of other immune pathways in the high-risk group is lower than that in the low-risk group (Figure 5B).

Figure 5 Comparison of the level of immune cell infiltration and the enrichment score of immune pathways between the two risk groups. (A) Comparison of 16 types of immune cells between the low-risk (green box) and high-risk (red box) groups; (B) Enrichment score of 13 immune-related pathways between the low-risk (blue box) and high-risk (red box) groups. *P<0.05, **P<0.01, ***P<0.001.

Discussion

In this study, we first screened the differentially expressed p53 pathway-related genes in normal colon tissue and colon cancer tissue and obtained two clusters through consistent cluster analysis of DEGs. We found that there are significant differences in PFS between the two clusters. Unfortunately, The difference in OS and DSS between the two clusters is not significant, and there is no significant difference in clinicopathological information between the two clusters. To further explore the impact of these DEGs on PFS, three prognostic genes (CDKN2A, BAK1, BTG1) were obtained through univariate COX regression analysis, and one independent prognostic gene (CDKN2A) was obtained by multivariate COX regression analysis of the above three genes. The risk score of each patient was calculated based on 3 prognostic-related DEGs, and the patients were divided into high/low risk score groups based on the median risk score. There were significant differences in PFS between the high/low risk score groups. In addition, after matching the patient’s risk score group with gender, age and stage, univariate and multivariate COX regression analysis, we found that the risk score group can be used as an independent prognostic factor for PFS in colon cancer patients. The distribution of gender, age, and tumor stage between the high/low risk score groups was uniform. Functional analysis showed that DEGs between high/low risk groups were related to neutrophil-mediated immune-related pathways, chemokine signaling pathways, Rap1 signaling pathways, and MAPK signaling pathways. Comparing the level of infiltrating immune cells and the activity of immune-related pathways between the high/low risk groups, we did not find a general rule.

The p53-mediated cell signal transduction pathway plays an important role in mediating the normal life activities of cells. At the same time, the p53 gene is the gene with the highest correlation with human tumors so far, and its correlation with colon cancer has been widely studied and reported.^2,3,16 However, in colon cancer, the relationship between p53 pathway-related genes and patient prognosis is still unclear. In our study, we identified 3 (CDKN2A, BAK1, BTG1) genes that are significantly related to PFS, and we believe that CDKN2A gene is an independent prognostic factor for PFS. The CDKN2A gene is also known as the multiple tumor suppressor gene. In 1993, Serrano et al discovered the Cyclin-dependent kinase inhibitor 2A (CDKN2A) gene when they used the yeast two-hybrid protein correlation screening method to study proteins that interact with cyclin dependent kinase 4 (CDK4).¹⁷ Early studies have shown that point mutations, deletions and methylation of CDKN2A gene may be related to cell senescence and cancer progression.¹⁸ In recent years, the close relationship between the expression level of CDKN2A gene and the occurrence and prognosis of gastrointestinal tumors has been reported successively,^19,20 which are similar to our findings. CDKN2A gene β transcription product p14ARF can activate the p53 response, which is manifested by increased levels of MDM2 and p21CIP1 and cell cycle arrest in G1 and G2/M phases. In contrast, p53 gene negatively regulates p14ARF expression.²¹ BRI1-associated receptor kinase 1 (BAK1) is considered to be a gene that promotes apoptosis. Research by Shi et al showed that in bladder cancer, the down-regulation of BAK1 mediated by miR-125b may contribute to the progression of bladder cancer and resistance to treatment.²² In addition, in breast cancer, miR-125b-mediated down-regulation of BAK1 inhibits paclitaxel-induced apoptosis and leads to drug resistance. Restoring the expression of BAK1 can overcome miR-125-mediated paclitaxel resistance.²³ The down-regulation of BAK1 is also associated with the occurrence of cervical cancer and non-small cell lung cancer.^24,25 Regarding colon cancer, the latest research clarifies that miR-410 inhibits the expression of BAK1 and inhibits the pathway of cell apoptosis, leading to the occurrence of colorectal cancer.²⁶ Our research further validates the previous research. p53 gene can induce the expression of BAK1, but the specific regulatory mechanism is still unclear.²⁷ B-cell translocation gene 1 (BTG1), as an anti-proliferative gene, was first reported in B-cell chronic lymphocytic leukemia.²⁸ However, Kamalakaran et al found that the high expression of BTG1 can promote the angiogenesis of breast cancer tumors.²⁹ On the other hand, BTG1 is believed to inhibit breast cancer cell growth by inducing cell cycle arrest and apoptosis, and can reduce the possibility of metastasis.^30,31 Therefore, the BTG1 gene may have different effects on different tissues. With further research, the relationship between down-regulation of BTG1 and the occurrence or prognosis of gastric cancer, bladder cancer, and kidney cancer has been gradually recognized.^32–34 The latest meta-analysis proves that the BTG1 gene can be used as an important biomarker for solid tumors.³⁵ At present, many studies have confirmed that the low expression of the BTG1 gene promotes the occurrence and invasion of colon cancer through a variety of ways and affects the prognosis of patients.^36–40 It is a pity that there is no research on the regulatory relationship between the BTG1 gene and the p53 gene.

Until now, many mechanisms affect tumor progression and prognosis, and there may be interactions and coexistence among various mechanisms. The deletion or mutation of p53 gene not only promotes the malignant transformation of cells by changing the control of the cell cycle, but also promotes tumor growth by inhibiting the components of the immune system in the tumor microenvironment.⁴¹ GO analysis of DEGs in the high/low risk score group suggests that the immune response mediated by neutrophils may play a key role in the p53 gene-related pathway. In the immune cell enrichment analysis, we found that the neutrophil infiltration level of the high-risk group was significantly higher than that of the low-risk group. However, early studies have shown that in the early stages of tumor development, neutrophils can produce reactive oxygen species (ROS), and ROS-induced apoptosis signals can kill tumor cells.⁴² But Uribe-Querol et al believe that when ROS produced by neutrophils is not enough to kill tumor cells, it will promote tumor proliferation through genotoxicity and DNA damage.⁴³ Therefore, neutrophils may play a double-edged sword in the process of tumorigenesis. The current mainstream view is that the deletion or mutation of p53 gene can lead to the decrease of MHC-I expression, resulting in a decrease in the activity of NK cells and a decrease in the level of Treg cell infiltration in the tumor microenvironment, leading to a poor prognosis.⁴⁴ In addition, type I interferon can promote the tumor suppressor function of p53, and promote the infiltration of NK cells, macrophages and cytotoxic T cells in the tumor microenvironment, and enhance the killing of tumor cells.^45,46 In this study, the opposite result was obtained. We believe that there are three possible reasons. First, there are differences in the immune infiltration components and the activity of immune pathways between the recurring tumor and the primary tumor. Second, the level of immune cells and the activation of immune pathways in the tumor microenvironment are the result of a comprehensive influence of multiple factors, and the immune system has complex internal controls. Third, although the level of immune cell infiltration is high, the activity of immune cells decreases or the tumor produces immune escape. This intricate immune microenvironment is also one of the important reasons for tumor proliferation, metastasis and recurrence, and further research is imminent.

Our research shows that three p53 pathway-related genes are significantly related to PFS in patients with colon cancer, and CDKN2A has been identified as an independent prognostic factor. In addition, we found that DEGs in the high/low risk score group were significantly related to the infiltrating cells of the tumor immune microenvironment and the activity of immune pathways. However, this research still has shortcomings. First, our data comes from public open-source databases, and the data may have a certain offset. Second, limited by the availability of data, all factors related to the prognosis of colon cancer have not been included in the study, such as tumor grade, nerve/vascular invasion, tumor location, etc. Third, there is no verification of the results.

Conclusion

p53 pathway-related genes play an important role in regulating the immune microenvironment of colon cancer and are significantly related to PFS in colon cancer patients.