Analysis of core genes for colorectal cancer prognosis based on immune and stromal scores

Introduction

Colorectal cancer (CRC) is a common malignant tumor (Long, Lundsmith & Hamilton, 2017). Despite advances in understanding the molecular mechanisms of CRC, it still has a high mortality rate. Guided treatment based on outcome prediction is an effective method to reduce the mortality of CRC (Yu et al., 2017). In recent years, studies on the prognostication of CRC have become increasingly important, and a large number of potential outcome predictors have been identified (Sanz-Garcia et al., 2017; Rena et al., 2002). The pathogenesis of CRC is comprehensive. Tumor cell intrinsic genes, especially master transcription factors, dictate the initiation, progression, and evolution of CRC (Wang et al., 2017; Suvà et al., 2014). The tumor microenvironment has also been reported to critically influence the expression of genes in tumor tissues and the prognosis (Curry et al., 2014; Cooper et al., 2012).

The tumor microenvironment is the cellular environment of tumor cells, consisting of extracellular matrix, soluble molecules and tumor stromal cells. In the tumor microenvironment, immune cells and stromal cells are the two main types of nontumor components and have been proposed to be valuable for the diagnosis and prognosis evaluation of tumors (Fukumura et al., 2010). In 2011, scientists began to examine various stages (TNMI-IV) of CRC for immunity to authenticate grading and staging systems, and there was considerable research on immune scoring stages and TNM staging. The results showed that an immune scoring system could predict the outcomes of CRC more accurately. It has obvious advantages for predicting survival in patients with CRC (Mlecnik et al., 2011; Fridman et al., 2012; Galon, 2006). According to the study conducted by Yoshihara et al. (2013), immune and stromal scores could be used to predict the infiltration of nontumor cells by analyzing specific gene expression data from The Cancer Genome Atlas (TCGA) database of immune and stromal cells.

However, no research has been conducted on predicting the prognosis-related genes of CRC by using the tumor microenvironment score. For the first time, in this current work, by taking advantage of both the TCGA database of CRC cohorts and the ESTIMATE algorithm-derived immune scores, we extracted a list of microenvironment-associated genes that predict the outcomes of CRC patients.

Materials & Methods

Results

The matrix score was significantly correlated with stage T colorectal cancer

CRC expression profile data were downloaded from the TCGA database. Clinical information (OS>1 month) and immune and stromal scores calculated by the ESTIMATE algorithm were integrated. In addition, 362 CRC samples were finally retained for subsequent analysis. The distribution of the immune score and stromal score is shown in Fig. 1A. To determine the potential correlation between the overall survival rate and the immune and stromal scores, Maxstat software was used to find the optimal cut-offs for the immune and stromal scores. Cancer samples were divided into high and low scores according to the score of the optimal cut-off (stromal score = −207.8886/immune score = −708.0627), and the corresponding clinical information was used for survival analysis. In the survival curve of the immune score and matrix score, we found that the immune score and matrix score were significantly correlated with the survival time of the patients (p = 0.0035/ p = 0.0034, Figs. 1B–1C). At the same time, the results showed that the stromal score had no statistical significance with all stages, but it was significantly correlated with T staging in the correlation analysis (p = 0.28/ p = 0.029, Figs. 1D–1E). There was no statistical significance in the correlation between the immune score and staging. However, there was a certain difference in the scores for each stage (p = 0.34/ p = 0.93, Figs. 1F–1G).

Comparison of colorectal cancer gene expression with immune score and stromal score

To reveal the correlations between gene expression and the immune and matrix scores, we evaluated the expression profile data of CRC patients in the TCGA database. A total of 228 DEGs were identified according to the matrix score (100 cases with a high score and 262 cases with a low score), among which 163 genes were upregulated and 65 genes were downregulated. The results showed that some of the DEGs were significantly different between the two groups (Fig. 2A). A total of 579 DEGs were identified according to the immune score groups (38 cases with high scores and 324 cases with low scores), among which 301 genes were upregulated and 278 genes were downregulated. These 789 DEGs were used as key genes for subsequent analysis (Fig. 2B), and their details are listed in Table S1. To clarify the potential function of the 789 DEGs, GO and KEGG enrichment analyses were performed on the DEGs. The results are shown in Figs. 3A and 3B (the top 8 were selected for display). GO analyses revealed that these DEGs were mostly involved in tissue-specific immune responses, positive regulation of adaptive immune responses and so on. In addition, according to the KEGG pathway analyses, the DEGs were involved in the age-race signaling pathway in diabetic complications, cytokine-cytokine receptor interaction, IL17 signaling pathway and so on. The results of the enrichment analysis of KEGG pathways in GSEA showed that the genes were mostly enriched in the pathways of asthma, graft-versus-host disease, malaria, rheumatoid arthritis, and viral protein interaction with cytokines and cytokine receptors based on the ImmuneScore (Fig. 4A). In addition, the results showed that the genes were mostly enriched in the ECM-receptor interaction, hypertrophic cardiomyopathy, protein digestion and absorption, Staphylococcus aureus infection, and systemic lupus erythematosus pathways based on the Stromalscore (Fig. 4B).

Survival analysis of DEGs

To screen out the genes associated with the prognosis of CRC, we divided the 789 DEGs into high and low expression groups according to their median expression and performed survival analysis. A total of 121 survival-related prognostic genes were extracted (p < 0.05). The survival curves of some of the prognostic genes are shown in Fig. 5A.

Construction of PPI networks between genes

To better understand the interactions between the identified prognostic genes, we used the STRING database to obtain a PPI network consisting of 35 nodes and 36 edges (Fig. 6). We identified tightly linked modules (green parts) in the network and defined the tightly coupled module genes as the core prognostic genes (N = 22). GO analysis and KEGG enrichment analysis were performed on 22 core prognostic genes. The results of the enrichment analysis are shown in Figs. 7A–7B (the top 10 genes were selected for display).

Regulation of the core prognostic genes by ncRNA/TF

Based on the 51913 ncRNA-mRNA interaction relationship in the RAID2.0 database, we searched for the pivot node (ncRNA) of the regulatory function module. When the p value < 0.05, we obtained four ncRNAs that regulated the core prognosis gene (Table 3), including metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) (associated with tumor cell proliferation and metastasis), colorectal neoplasia differentially expressed (CRNDE) (promotes colorectal cancer cell proliferation), FOXF1 adjacent noncoding developmental regulatory RNA (FENDRR) and taurine upregulated gene 1 (TUG1) (regulation of resistance to colorectal cancer methotrexate).

Table 3:

Pivot (ncRNA).

ncRNA	Connection	P value
MALAT1	7	0.000346134
CANDE	6	0.000878654
FENDRR	8	0.000451542
TUG1	7	0.000286648

DOI: 10.7717/peerj.12452/table-3

Notes:

Abbreviations

MALAT1

Metastasis-associated lung adenocarcinoma transcript 1

CRNDE

colorectal neoplasia differentially Expressed

FENDRR

FOXF1 adjacent non-coding developmental regulatory RNA

TUG1

taurine upregulated gene1

Based on the interaction relationships of 9396 human TF-mRNAs contained in the TRRUST v2 database, the pivot node (TF) of the regulatory function module was searched. When p value < 0.05, 11 TFs with regulatory effects on the core prognostic genes were obtained (Table 4), including a number of cancer-related transcription factors, such as activating transcription activator 1 (ATF1), CCAAT/enhancer binding proteins C/EBP beta (CEBPB), E2F transcription Factor 3 (E2F3), ETS translocation variant 4 (ETV4), and Y-box protein 1 (YBX1). The visualization of the core prognostic genes and ncRNA/TF interactions is shown in Fig. 8A.

Table 4:

Pivot (TF).

TF	Connection	P value
ATF1	3	0.000124081
CEBPB	4	0.000733944
E2F3	2	0.00026685
ETV4	2	0.001742111
PAX3	3	6.68E−07
SOX10	3	1.63E−05
SP11	2	0.035007145
DSF1	2	0.40192337
USF2	2	0.18026951
YBX1	2	0.004397422
YY1	3	0.021654705

DOI: 10.7717/peerj.12452/table-4

Notes:

Abbreviations

ATF1

activating transcription activator 1

CEBPB

CCAAT/enhancer binding proteins C/EBP beta

E2F3

E2F transcription factor 3

ETV4

ETS translocation variant 4

YBX1

Y-box protein 1

PAX3

Paired box 3

SOX10

SRY-related HMG-box 10

SP11

Synaptophysin 11

DSF1

Double-skin façade 1

USF2

Upstream stimulus factor 2

YY1

Yin Yang 1

Mutation analysis of the hub genes

As shown in Fig. 9, the 22 hub gene samples had mutations only in 90 samples. Among these, ERBB2 had the highest mutation rate, most of which were missense mutations, and most of the samples had mutations of C > T.

The ROC curve based on the hub genes

The area under the curve (AUC) of the 22 hub genes ranged from 0.344 to 0.669. Twelve hub genes, ENO2, GRP, HIST1H2AC, HIST1H2AK, HIST1H4H, HIST4H4, HOXB4, HSPB1, MIFT, POMC, QRFP, and TNFAIP6, had weak value in predicting the survival of colorectal cancer with >0.5 AUCs (Fig. 10).

Discussion

In the current work, we identified tumor microenvironment-related genes that contribute to CRC overall survival in the TCGA database and found that both the immune score and matrix score were significantly correlated with the patient survival time. Importantly, we were able to validate 789 DEGs. In addition, 121 prognostic genes were found to be associated with survival. Twenty-two genes were identified as core prognostic genes. Four ncRNAs that regulate the core prognostic genes, 11 TFs that had a regulatory effect on the core prognostic genes, and 2 drugs that have regulatory effects on CRC were also identified.

Both the immune score and the matrix score can predict the purity of the tumor, as well as the number of stromal and immune cells. The more stromal and immune cells there are, the lower the purity of the tumor. In a study of immune scoring by Galon et al., patients with early TNM (stage I and II) were followed up for survival and recurrence of colon cancer, and it was found that patients with high immune scores had a longer survival period. Ninety-five percent of patients with a high score had no tumor recurrence within 18 years after surgery. Fifty percent of patients with low scores had a recurrence within 2 years after surgery (Galon et al., 2012). In 2018, Pages et al. reported that immune scores could provide a reliable assessment of the risk of recurrence of colon cancer, which supports the use of immune scores as part of the new TNM immunotyping (Pagès et al., 2018). Therefore, as a perfect supplement to the TNM staging system, the immune score will be a superior staging system for predicting the survival of tumor patients.

Searching for more metastatic genes related to CRC and studying their biological characteristics and the specific mechanism of metastasis can help control tumor metastasis, prevent recurrence and improve the prognosis of patients. In terms of multigene early diagnosis and a prognostic model of CRC, Marshall et al. obtained an early diagnostic model of CRC based on the gene expression profile of the peripheral blood of CRC patients (Marshall et al., 2010). Sveen et al. (2012) proposed a prognostic model by studying gene expression in tumor tissues, which has potential prognostic value in stage IV and V CRC patients. In addition, a large number of genes, including P53 (Lüchtenborg et al., 2004), K-ras (Andreyev et al., 1990), B-raf (De Roock et al., 2010), and N-ras (Parsons et al., 2012), have been shown to be associated with the prognosis of CRC.

In the present study, 22 genes were identified as the core prognostic genes of CRC. Among them, MAPK8 (Slattery, Lundgreen & Wolff, 2012), HSPB1 (Nadin et al., 2012), IL1B (Sanabria-Salas et al., 2017), PRIM1 (Cloutier et al., 1997) and so on were shown to be associated with the susceptibility to or prognosis of CRC. However, only the PRIM1 gene was validated in the GEO database.

PRIM1 (primase) plays a key role in the process of DNA synthesis initiation by synthesizing RNA primers for Okazaki fragments (Cloutier et al., 1997). Mutations in PRIM1 cause extensive apoptosis of retinal neurons through activation of the DNA damage checkpoint and tumor suppressor p53 (Yamaguchi et al., 2008). PRIM1 plays an important role in the strict control of the DNA replication fork during tumor cell proliferation and was shown to be involved in estrogen-induced breast cancer formation through activation of the G2/M cell cycle checkpoint (Lee et al., 2018). In addition, PRIM1 (Job et al., 2018) has been reported to be associated with CRC progression, and this gene is likely to be a prognostic marker of the CRC immune microenvironment. However, the mechanism by which PRIM1 affects CRC is not clear.

A total of four core ncRNAs were obtained that could regulate the core prognostic genes in the present study. MALAT1 is a long noncoding RNA with important biological functions and is located on chromosome llql3.1 (Zhang et al., 2017). MALT1 expression is upregulated in CRC tumor cells. A-kinase anchoring protein 9 (AKAP9) is a target gene for the cancer-promoting effect of MALAT1. AKAP9 is related to the polarization of cell spinning chains, and the overexpression of MALAT1 in normal intestinal epithelial cells causes them to divide and differentiate and promotes tumorigenesis (Xu et al., 2011). In vitro experiments have demonstrated that its knockout can inhibit the migration of p-catenin from the cytoplasm to the nucleus, resulting in decreased expression of c-MYC and MMP-7, whereas overexpression of MALAT1 can activate the Wnt/p-catenin pathway and promote the invasion and metastasis of CRC (Ji et al., 2013). Zheng et al. (2014) further confirmed that MALAT1 was significantly upregulated in CRC tissues and was associated with a poor prognosis in stage II/III CRC patients.

CRNDE was first identified as a biomarker of CRC. High expression of IncRNA CRNDE was detected in more than 90% of colorectal adenomas and adenocarcinoma cells (Ning et al., 2014). Experiments showed that the regulatory effect of CRNDE was related to the IGF signaling pathway and was involved in the carcinogenic process of intestinal epithelial cells (Huan et al., 2017). CRNDE can activate the Wnt/x-catenin signaling pathway in CRC patients, thus promoting tumor cell proliferation, invasion and other processes (Clevers & Nusse, 2012). In addition, Li et al. (2018) showed that CRNDE can be used as a potential noninvasive serum marker to predict the efficacy of first-line FOLFOX regimens for metastatic colorectal cancer and it is associated with a poor prognosis. Although CRNDE can enhance the proliferation, migration and invasion of CRC cell lines, the mechanism by which CRNDE promotes the proliferation and migration of CRC cells remains unclear.

Taurine upregulated gene 1 (TUGl) is a newly discovered carcinogenic lncRNA located on chromosome 22q12 (Li et al., 2016). Numerous studies have demonstrated that TUGl expression levels are significantly elevated in CRC tissues and cell lines. The expression level of TUGl in CRC tissues was 4-6 times higher than that in adjacent noncancerous tissues (Sun et al., 2016). In vitro tests, overexpression of TUGl can induce the formation of tumor cell colonies and activate the expression of EMT-related genes to improve the invasion and metastasis of tumors and promote liver metastasis of CRC (Sun et al., 2018). TUGl promotes cancer cell proliferation, migration, invasion, and epithelial-mesenchymal transition (EMT) and inhibits tumor cell apoptosis through complex mechanisms, including competitive inhibition of miRNA function. Nevertheless, the specific oncogenic mechanism of TUGl remains to be further elucidated (Huang et al., 2016).

Transcription factors activate the transcription of downstream target genes by binding to the target gene promoter region. Studies have found that transcription factors are closely related to tumors. A total of 11 TF genes regulating the core prognosis were obtained in the present study. We are particularly interested in ATF1, CEBPB, E2F3, and ETV4. ATF1 belongs to the cAMP response element binding protein (cyclic AMP response element binding protein, after CREB) family (Huang et al., 2016). A study found that ATF1 in a wide variety of tumors plays a role as both an oncogene and a tumor suppressor gene by influencing the cell signal transduction pathways involved in tumor cells related to biological processes, such as proliferation, apoptosis, angiogenesis, migration, invasion, and immune surveillance, which affects the occurrence and development of tumors (Pu, Storr & Ahmad, 2018). In mouse colorectal cancer cells, active phosphorylated CREB is elevated (Sampurno et al., 2013), and p300 can assist the phosphorylation of CREB and activate intestinal stem cell transcription factors, including Myb, that regulate the proliferation of intestinal epithelial cells (Ramaswamy et al., 2018). It was speculated that ATF1 is involved in human colorectal cancer cell proliferation through the p300-MYB-CREB axis.

The transcription factor CEBPB is involved in a number of biological processes, including cell differentiation, metabolic balance, proliferation, tumorigenesis, apoptosis, immune and stress responses, energy metabolism, and blood production (Ramji & Foka, 2002). Some studies have shown that CEBPB affects tumorigenesis by interacting with other genes to form a regulatory network during tumorigenesis and development (Abreu & Sealy, 2010). Other studies have also shown that the LIP subtype of the C/EBP transcription factor can induce apoptosis of human breast cancer cells and induce its own phagocytosis, which may be acting as a tumor suppressor by inducing tumor autophagy (Koslowski et al., 2009).

Additionally, the E2F transcription factor is an important regulator of G1 phase entry into S phase in the cell cycle, and it is closely related to tumor occurrence and cell apoptosis (Leone et al., 1998). Hurst et al. (2008) showed that E2F3 inhibited pRb and P53 through two important pathways of cell proliferation regulation and tumor monitoring, p16(Ink4a)-cycd/cdk4-rb-e2f and Arf/mdm2-p53, respectively, leading to the disorder of cell cycle regulation and thus promoting the occurrence and development of tumors. Wang, Zhao & Yuan (2011) found that the expression of E2F3 was high in colorectal carcinoma. There is a positive correlation with the expression of E2F3 in colorectal carcinoma.

Furthermore, ETV4 is a member of the polyomavirus enhancer activator 3 (PEA3) subfamily of ETS transcription factors, in which transcription factors can recognize and bind to GGAA/T sequences to regulate the expression of multiple target genes, thus affecting the occurrence and development of diseases. Existing studies have found that ETV4’s ability to promote cancer migration is closely related to MMP (Fung et al., 2016). For example, in breast cancer, ETV4 can promote the migration of cancer cells by promoting the expression of MMP2 (Bièche et al., 2004); in esophageal cancer, ETV4 can promote the metastasis of cancer cells by promoting the expression of MMP1 (Keld et al., 2010). In addition, ETV4 was found to promote the invasion and metastasis of cancer cells by promoting the greening of COX2 (Subbaramaiah, 2002). Of course, the mechanisms of the pro-proliferative and pro-migratory effects of ETV4 in CRC are not fully understood.

Notably, we also predicted that quercetin might be a drug that regulates CRC. Quercetin is a flavonoid with a wide range of biological activities, including antioxidant activity, scavenging oxygen free radicals, antifibrosis activity, lowering blood pressure, lowering blood glucose, protecting heart muscle and antitumor activity (Nguyen et al., 2017). In recent years, many studies have found that quercetin can inhibit a variety of cancer cells, especially colon cancer. Quercetin can not only inhibit the proliferation and induce the apoptosis of colon cancer cells but also reduce the number of abnormal gland crypts in the colon (Miyamoto, Yasui & Ohigashi, 2010). However, there are few studies on the effect of quercetin on colon cancer, and the specific molecular mechanism is still not fully understood. Luo et al. (2014) reported that quercetin may inhibit cell proliferation and induce apoptosis through the bcl-2 and c-myc genes. Park et al. (2005) found that quercetin could significantly inhibit the transcriptional activity of catenin/Tcf signals in SW480 colon cancer cells by reducing the protein levels of large-catenin and Tcf in the nucleus. Thus, quercetin can be used as an adjuvant drug to inhibit the growth of colon cancer cells and has potential value in the drug treatment of colon cancer.

Conclusions

In summary, from the functional enrichment analysis of the TCGA database applied by ESTIMATE algorithm-based immune scores, we extracted a list of tumor microenvironment-related genes. These genes could be useful for determining the prognosis of CRC patients. In addition, further investigation of these genes as well as their regulators, including ncRNAs and TFs, provides a stronger predictor of survival than individual genes. Finally, the detection of drugs that regulate CRC could lead to novel insights into potential treatments for CRC.

Supplemental Information