Using the TCGA Database to Predict and Analyze Tumor Microenvironment Genes Related to Poor Prognosis of Colon Cancer

Background

COAD is a malignant cancer of the digestive tract. It can be divided into right colon cancer and left colon cancer according to the location of the cancer, but they both are lower digestive tract cancers. According to a WHO study in 2018, there were approximately 1.8 million COAD patients worldwide, accounting for 10.2% of all tumors. The number of deaths due to COAD in 2018 reached 881 000 [1–3].

There are many types of treatments for COAD, such as surgical removal of cancerous lesions, target gene chemotherapy, and radiation therapy. However, there is currently no therapeutic means to completely cure colon cancer in all cases. There is therefore the need to further explore colon cancer treatment options to obtain a higher colon cancer cure rate. Because of the establishment of a tumor genome databases such as TCGA, we can better understand the impact of genetic composition on clinical prognosis, and can mine and analyze data from large cohorts around the world. According to information in the TCGA database, the historical type of colon cancer can be divided into 2 categories: Colon Mucinous Adenocarcinoma and Colon Adenocarcinoma. According to the NCCN2019 version of the guidelines, AFP, CEA, and CA199 mutations can be used as references for disease diagnosis [4]. Today, gene expression profiles are increasingly being recognized as part of the standard for clinical diagnosis.

According to research by Di Jia et al., the ESTIMATE algorithm-derived immune score and stromal score can be used for glioblastoma (GBM) patients. We used the same method to collect and extract a series of microenvironment-related genetic data in the TCGA database and further explored the effect on the colon genes for prognosis in colon cancer patients [5].

The purpose of this study was to use bioinformatics-related databases and tools to analyze and obtain relevant genes that constitute the tumor microenvironment and affect the survival prognosis of patients with colon cancer.

Material and Methods

RELEVANT DATA ACQUISITION:

We obtained data on 191 colon cancer patients from the TCGA database and data on RNA expression for colon cancer using AgilentG4502A_07_3 [6]. Clinically relevant data, such as age, sex, tumor pathology classification, survival rate, and histological type were obtained from the TCGA database. Relevant immune scores and stromal scores were obtained from the ESTIMATE algorithm database [7].

DIFFERENTIALLY EXPRESSED GENES (DEGS) INCLUSION CRITERIA:

Data that fold change ≤|1.0| and adj. p value ≤0.05 were included as selected genes.

HEATMAPS PRODUCTION AND ANALYSIS:

Data analysis and heat map production were performed using the ClusterVis database [8].

VOLCANO MAPPING:

Data with fold change ≤| 1.0 | and adj. P value ≤0.05 were included as selected standard, and then plotted in the volcano map.

CONSTRUCTION AND MODIFICATION OF PPI NETWORK:

We used the STRING database to build a PPI network, and used Cytoscape software to modify it [9,10]. We used the Cytoscape-related plug-in ClueGO [11] to further explore the potential immune system functions of DEGs.

CORRELATION SURVIVAL ANALYSIS AND SURVIVAL CURVE PLOTTING:

Survival-related data were analyzed using the log-rank test, and Kaplan-Meier plots were used to assess the relationship between related DEGs and survival data.

DEGS FOR GO; KEGG; IMMUNE FUNCTION ANALYSIS:

We used the DAVID database [12] for functional enrichment analysis of DEGs, such as GO analysis and KEGG pathway analysis. The false discovery rate (FDR) was considered statistically significant.

Results

CORRELATION BETWEEN HISTOLOGICAL TYPE OF COLON CANCER AND IMMUNE SCORES OR STROMAL SCORES:

We downloaded 191 colon cancer patient data from the TCGA database, such as patient age, tumor size, tumor grade, and tumor location. Patient ages ranged from 30 to 90 years. Of the 191 patients, 105 (55%) were male and 86 (45%) were female. Histological typing of 191 patients revealed that 48 (25%) had colon mucinous adenocarcinoma and 143 (75%) had colon adenocarcinoma. Patients with different histology types had different prognoses. We analyzed the selected patient data based on the ESTIMATE algorithm and obtained immune scores and stromal scores. Stromal scores ranged from −2053.75 to 976.26, and immune score ranged from −477.96 to 1511.97. Colon mucinous adenocarcinoma had higher stromal scores and immune scores than colon adenocarcinoma. The results show that both stromal score and immune score have significance in the correlation of historical type (Figure 1A, 1B).

According to the 2019 NCCN Clinical Practice Guidelines for Colon Cancer, CEA, AFP, and CA199 are related to the diagnosis and prognosis of colon cancer. However, the data we collected only contained data on AFP gene mutations, so we used AFP mutations for verification. According to the related gene mutation chart, the results show that the stromal score and immune score of the genetic mutation cases have no statistical significance (Figure 1C, 1D).

We divided the immune score of 191 patients in half into a high immune score group and a low immune score group, and compared the survival prognosis between the high immune score group and the low immune score group. The results showed that the survival rate of the high immune score group was significantly lower than in the low immune score group, and demonstrated a significant relationship between immune score and survival prognosis (p=0.0182). Using the same method, we compared the difference in survival prognosis between high and low matrix fractions. The results show that before 1000 days, the survival rate of the high stromal score group was significantly lower than that of the low stromal scores group, but after 1000 days there, they were similar (p=0.7368) (Figure 1E, 1F).

RELATIONSHIP BETWEEN IMMUNE SCORE/STROMAL SCORE AND GENE EXPRESSION IN COLON CANCER:

We compared our downloaded 191 patient data with Affymetrix microarray data, and created heatmaps based on high and low immune score/stromal score groups for visual analysis (Figure 2A, 2B). According to the inclusion criteria fold change >1, adj. P<0.05, 2 volcano maps showed 212 upregulated genes and 7 downregulated genes in the immune score group, and 155 upregulated genes and 8 downregulated genes in the stromal score group (Figure 2C, 2D). The Venn diagram summary (Figure 2E, 2F) shows there are 78 upregulated genes and 2 downregulated genes in both the immune score group and the stromal score group. We therefore focussed on these differentially expressed genes (DEGs) [8,13].

To understand the potential biological functions of DEGs, we performed an enrichment analysis on 78 DEGs that were upregulated. The TOP10 of the gene ontology (GO) and KEGG pathway analysis showed that biological processes, cellular components, and cellular components are closely related to immune mechanisms, such as inflammatory response, immune response, and chemotaxis (Figure 3A–3D).

MAPPING OF PROTEIN-PROTEIN INTERACTION (PPI) NETWORKS FOR DEGS:

To better understand the interaction between DEGs, we used the STRING database to draw a PPI network. The PPI network was composed of 3 module groups, including 78 nodes and 415 edges (Figure 4). The circle size represents the size of the Degree value, and the color shade represents the size of the LogFC value. The thickness of the line represents the combined score between the proteins. Through the PPI network, we found that some immune response genes, such as CCL8, IL18RAP, P2RY14, FPR3, CD163, CCL2, FCGR2B, AIF1, TLR2, and TNFSF13B, are in a relatively core positions in the PPI network and play an important role in protein interactions. ClueGO immune system analysis of DEGs in the PPI network showed that negative regulation of dendritic cell differentiation occupies the main part (60.42%) (Figure 4A, 4B).

DEGS-RELATED SURVIVAL PROGNOSIS ANALYSIS AND PLOTTING SURVIVAL CURVES:

To further explore the relationship between DEGs and colon cancer survival rates, we used the TCGA database to download the data we needed for survival analysis and plotted Kaplan-Meier survival curves. The results showed that among 78 DEGs, 17 genes had poor overall survival (log-rank <0.05). Five of these genes belong to the immune response pathway (Figure 5A, 5B).

CELLULAR COMPONENT ANALYSIS OF DEGS WITH POOR SURVIVAL PROGNOSIS:

We performed related cellular component analysis of DEGs associated with poor prognosis of colon cancer. The results showed that genes were mainly concentrated in extracellular space, integral component of membrane, and plasma membrane. The specific gene distribution is shown in Table 1.

Discussion

In this study, we downloaded data on 191 patients from the TCGA database. We tried to identify genes associated with colon cancer tumor microenvironment that are correlated with colon cancer survival prognosis in the TCGA database. First, we distinguished between high-score and low-score immune score and stromal score groups, and assessed the intersection of the 2 high-score groups. A total of 78 genes were co-expressed in the immune score group and the stromal score group, and a total of 78 genes were co-expressed in the immune score group and the stromal score group (Figure 6).

Next, further biological function enrichment analysis was performed. The results showed that the functions of these genes include inflammatory response, immune response, chemotaxis, chemokine-mediated signaling pathway, and innate immune response and regulation of immune response (Figure 2). According to the results of biological enrichment analysis, the biological functions of these genes are closely related to the tumor microenvironment [14–28].

Finally, through the information downloaded from the TCGA database, we identified 78 genes related to the tumor immune microenvironment, which have a significant relationship with poor prognosis of colon cancer. The results showed that 17 genes were closely related to poor survival prognosis of colon cancer. Among them, there are 5 genes in the immune response, and the remaining 12 genes were also related to the tumor microenvironment. Analysis of 17 genes of cellular components showed that the genes of extracellular space were CCL2, TNFSF13B, LY96, CCL8, SELE, and SRGN, the genes of integral component membrane were MS4A4A, C3AR1, LAIR1, C1ORF162, IL18RAP, FCGR2B, TNFSF13B, LILRB3, CD52, and SELE, and the genes of plasma membrane types were C3AR1, LAIR1, IL18RAP, FCGR2B, TNFSF13B, LY96, SELE, and S100A12 (Table 1).

In this study, the genes we were interested in were CCL8, CCL2, FCGR2B, IL18RAP, and TNFSF13B. These genes are immune feedback function genes and belong to relatively core positions in the PPI network. Our review of the literature showed that CCL8 factor regulates colonic inflammation, and the inflammatory response is closely related to the immune microenvironment. Moreover, colitis is also linked to the worsening of colon cancer, and CCL8 is a target for immunotherapy in pancreatic cancer [19,29,30]. Related studies show that CCL2 is a tumor chemokine and a marker in the immune microenvironment of colon cancer [31]. IL18RAP factor plays an important regulatory role in colitis and Crohn’s disease, and can regulate interleukins IL-2, IL-7, and IL-15. These interleukins regulate the immune response and are associated with many factors, such as lymphocytes and macrophages [32]. Deepa Bedi et al. showed that TNFSF13B has s clinical significance in PD-L1 immunotherapy, but the specific mechanism needs further exploration. Importantly, these genes are significantly related to survival prognosis in colon carcinoma.

Colon cancer has been fully studied, especially in the area of correlation between gene expression and survival prognosis. Some of these experiments were completed in animal models, and some were performed in retrospective studies of patient data. However, research related to the tumor microenvironment needs larger and more comprehensive data to be thoroughly analyzed. With the development of genome sequencing technology, similar to the existence of TCGA database, we can better study and analyze the existing big data.

In summary, we downloaded the data of 191 patients from the TCGA database, analyzed and calculated the relevant scores through the ESTIMATE calculation method, and then gathered a list of genes related to the tumor microenvironment. These data have been independently verified in the colon cancer cohort. Some previously under-appreciated genes may play important regulatory roles in the tumor microenvironment. Finally, further research on these genes may lead to new understanding of colon cancer and the tumor microenvironment.

Conclusions

By using a series of bioinformatics tools and related algorithms, we screened out a group of genes that constitute the tumor microenvironment and affect the survival prognosis of gastric cancer patients. Among these genes, 17 genes with high expression reduce the survival time of gastric cancer patients, of which 5 genes are related to immune response function.