Identification of Hub genes with prognostic values in colorectal cancer by integrated bioinformatics analysis

2.1Data source

As one of the biggest collections of gene chips in the world, the GEO database is a comprehensive gene expression library at the National Center of Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/ geo/). Expression profiles of GSE21510 and GSE110224 based on the GPL570 platform and GSE74602 based on the GPL6104 platform were obtained from the GEO database. GSE21510 contained 123 CRC samples and 25 normal tissues; GSE110224 included 17 CRC samples and 17 normal tissues, and GSE74602 contained 30 CRC samples and 30 normal tissue. Figure 1 shows the overall flowchart of the study.

Figure 1.

The flow diagram of study.

2.2Differential expressed gene analysis

The differentially expressed genes (DEGs) were filtered using a threshold of adjusted P values < 0.05 and an absolute log2FC (fold change) > 1 by GEO2R online tools (http://www.ncbi.nlm.nih.gov/geo/geo2r) [22]. Venn analysis was used to select overlapping DEGs among the three datasets mentioned above.

2.3GO and KEGG pathway analysis

To further understand the biological functions and related signaling pathways of the DEGs, we performed Gene Ontology (GO) enrichment analysis and Kyoto Encylopaedia of Genes and Genomes (KEGG) pathway enrichment analysis using the Metascape (http:// metascape.org/) [23, 24]. An adjusted P value of < 0.01 was considered to identify the enriched terms, and the results were visualized.

2.4Protein-protein interaction (PPI) comprehensive analysis

PPI comprehensive analysis was performed using the online tool Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/) [25, 26]. The selected Hub genes were imported into the STRING, and the confidence score of > 0.4 was considered significant. Then, the PPI network information was obtained. Next, the data downloading from STRING was further analyzed and visualized using Cytoscape (version 3.8.2). The core function of Cytoscape is to provide basic functional layout and query network, and build a PPI network based on the combination of basic data into a visual network. The app MCODE plugged into Cytoscape was used to identify the paramount modules in the PPI network (MCODE score > 5, degree cutoff = 2, maximum depth = 100, K-core = 2, node cutoff = 0.2). Finally, the visual network graph was presented in a circle layout according to the score.

2.5Hub gene selection and analysis of survival and differential expression

The PPI analysis of DEGs was performed using the online tool STRING, and results with a minimum interaction score of 0.4 were visualized in Cytoscape software (version 3.8.2). The obtained results were used to select the Hub genes through the MCC algorithm in the Cytoscape cytoHubba application, and the top 30 genes were selected as potential Hub genes.

Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) is a web server that analyzes cancerous and normal gene expression profiles and interactions in The Cancer Genome Atlas (TCGA, https://genomecancer.ucsc.edu/) and the Genotype Tissue Expression (GTEx, https://gtexportal.org/ home/) projects [27]. The overall survival of colorectal cancer patients and the expression level validation of Hub genes were evaluated using GEPIA. P< 0.05 was defined as the significant threshold to indicate candidate Hub genes. Genes not statistically significant were removed (P> 0.05). We also analyzed the impact of the expression levels of Hub genes on the disease free survival (DFS) of CRC patients using GEPIA. In addition, the expression levels in CRC from the GEO database were analyzed by Graphpad Prism 8 and scatter plots were drawn. The Tumor Immune Estimation Resource (TIMER) database was also used to explore the gene expression in tumors and adjacent normal tissues [28]. The differences of gene expression among different subtypes of CRC (POLE, MSI, CIN and GS) were investigated using R software.

2.6Protein expression and receiver operating characteristic (ROC) curves

The protein expression levels of Hub genes in human normal and cancer tissues were determined using the Human Protein Atlas (HPA) database (https://www.prot einatlas.org/), which respectively displays the expression of proteins in cells, normal tissues, and cancerous tissues. Receiver operating characteristic (ROC) curves were obtained from the online tool Xiantao Academy (https://www.xiantao.love/).

2.7Gene mutation status and survival analysis

cBioPortal (http://cbioportal.org) is an open-access resource for exploring, visualizing and analyzing multi-dimensional cancer genome data [29]. The cBioPortal was used to analyze the type and frequency of Hub gene mutations, and the effects of mutation on the patients’ prognosis.

2.8Immune infiltration analysis

The TIMER database contains 32 cancers and 10,897 tissue sample information from the TCGA database. It can realize systematic analysis of the correlation between immune infiltrates and other wide spectrum of factors, including related gene expression in tissues and prognosis, gene mutations and copy number of cancer patients [28]. In this study, we evaluated the infiltration of immune cells (CD8⁺ T cells, CD4⁺ T cells, B cells, dendritic cells, macrophages, and neutrophils) in CRC patients through TIMER database. The relationship between gene expression and tumor purity was also explored and visualized.

2.9Co-expression analysis

The Linked Omics database (http://www.linkes.org/) is a web-based platform for analyzing 32 TCGA cancer-related cubes [30]. The Link Finder module of Linked Omics was used to study the differentially expressed genes related to selected Hub genes in the TCGA-COAD, and the Pearson correlation coefficient was used for statistical analysis. All results are presented graphically in a volcano map, heat map or scatter plot. The Link-Interpreter module of Linked Omics performs pathway and network analysis of differentially expressed genes. Use the comprehensive functional classification database in the Web-based Web Gestalt to sign and sort the data in the Link Finder results, and use Gene set enrichment analysis (GSEA) to analyze the GO (CC, BP, MF) and KEGG channels. The GSEA program was run with 500 simulations, and the significance level was top 25. P value and false discovery rate (FDR) were both less than 0.05, the gene set was considered significantly enriched.

2.10Statistical analyses

Data are represented as the mean ± standard deviation (SD), and the t-test was used for comparisons between the two groups. GraphPad Prism 8.3.1 (GraphPad Software, Inc., San Diego, CA, United States) was utilized for statistical analysis and generating graphs. P< 0.05 was considered statistically significant.

3.1Identification and functional enrichment analysis of DEGs

We obtained 170 CRC tissues and 72 normal tissues from the GSE21510, GSE74602 and GSE110224 datasets. According to the cutoff criteria, a total of 7148 DEGs were screened from these three datasets. Among them, 4722, 1650, and 776 DEGs were obtained from the GSE21510, GSE74602, and GSE110224 datasets, respectively. Then 346 DEGs common to the three GEO datasets were detected via the Venn diagram, including 117 upregulated genes and 229 downregulated genes were identified (Fig. 2A and Table 1).

GO and KEGG pathway analysis were used to explore the biological functions and signal pathways of DEGs related to CRC. GO analysis indicated that the DEGs were mainly involved in the biological processing of the mitotic cell cycle process, regulation of hormone levels and response to nutrient levels (Fig. 2B), and associated with the cellular components including the apical part of cells, extracellular matrix and apical plasma member (Fig. 2B). The DGEs genes were linked to structural molecule activity, and signaling receptor activator activity (Fig. 2B). KEGG pathway analysis showed DEGs related to CRC were involved in bile secretion, pyruvate metabolism, viral protein interaction with cytokine and cytokine receptor, IL-17 signaling pathway, etc. (Fig. 2B).

Figure 2.

Differentially expressed of DEGs were selected from three GEO datasets. A, The DEGs with P< 0.05 and absolute log2FC > 1 were filter from three GES data sets using GEO2R online tools, and Venn analysis was used to select overlapping up-regulated and down-regulated DEGs. B, GO and KEGG pathway enrichment analysis were used to analyzed the biological process, cellular component, molecular function and signal pathways of DEGs.

3.2PPI network construction and Hub genes analysis

The PPI network diagram of DEGs was constructed by STRING, and analyzed and visualized by Cytoscape. The PPI network included 295 nodes and 1170 edges (Fig. 3A and B). Finally, we identified the most important nodes by the MCC algorithm. The new network included 30 nodes and 424 edges, and the top 30 central nodes were identified as potential Hub genes (Fig. 3C).

Figure 3.

The Hub genes were identified. A, STRING was used to constructed PPI networks of DEGs. B, Cytoscape was used to analyzed and visualized of PPI networks. Upregulated genes are shown in red and downregulated genes are shown in blue. C, The top 30 genes were selected using the MCC algorithm by Cytoscape’s plug-in cytoHubba. D–G, GEPIA was used to analyzed the correlation between overall survival and the expression of AURKA (D), CCNA2 (E), CCNB1 (F) and EXO1 (G) in colorectal cancer. H–K, The association between disease-free survival and the expression of AURKA (H), CCNA2 (I), CCNB1 (J) and EXO1 (K) in colorectal cancer by GEPIA.

Then, we further analyzed the relationship between the expression of 30 potential Hub genes and the overall survival of patients with CRC using GEPIA, including 270 CRC patients. According to the screening criteria (P< 0.05), four genes, Aurora kinase A (AURKA), Cyclin B1 (CCNB1), EXO1 and Cyclin A2 (CCNA2), were selected as candidate Hub genes. From the results of GEPIA analysis, patients with high expression of these four genes have higher overall survival than those with low expression (Fig. 3D–G). Furthermore, we next investigated the influence of these four genes on DFS in patients with CRC. The results showed that the patients with low expression of AURKA, CCNA2 and CCNB1 had lower DFS than those with high expression (Fig. 3H–K).

We verified the four Hub gene expression levels in colon adenocarcinoma (COAD) patients and normal via GEPIA, which included 275 tumors and 349 normal tissue. The result showed the expression levels of AURKA, CCNB1, EXO1 and CCNA2 significantly increased in tumors than those in normal tissue (P< 0.05, Fig. 4A). Following, the different expressions of AURKA, CCNB1, EXO1 and CCNA2 in CRC and normal tissues in GSE21510, GSE110224 and GSE74602 databases were analyzed. The results showed that the expression level of AURKA, CCNB1, EXO1 and CCNA2 was higher in CRC tissues than that in normal tissues (P< 0.001, Fig. 4B). Overall, the AURKA, CCNB1, EXO1 and CCNA2 genes were highly expressed in CRC tissues.

Subsequently, we investigated the expression of these four genes in various tumor tissues using the DiffExp module in TIMER 2.0. The statistical significance of differential expression was evaluated by the Wilcoxon test. We found that compared to normal tissues, AURKA, CCNB1, EXO1 and CCNA2 mRNA expression were significantly increased in COAD tumor tissues (Fig. 4C).

The HPA database was used to investigate the protein expression of AURKA, CCNB1, EXO1 and CCNA2 in CRC tissues. The protein expression of AURKA, CCNB1 and CCNA2 were increased in CRC tissues than that in normal tissues. The EXO1 protein expression was not found in tissue from the HPA database (Fig. 5A). ROC curves displayed that the AUC of AURKA, CCNB1, CCNA2 and EXO1 was 0.958, 0.929, 0.881, and 0.942, respectively (Fig. 5B). It indicated that AURKA, CCNB1, CCNA2 and EXO1 have good diagnostic value for CRC. The differences of gene expression among different subtypes of CRC were investigated. The results displayed that there are differences in the expression levels of Hubs among different subtypes of CRC (POLE, MSI, CIN, and GS). AURKA was a higher expression in CIN subtypes. The CCNA2 expression was highest in the POLE subtypes of four types. The CCNB1 and EOX1 expressions were higher in MSI and POLE subtypes than in CIN and GS subtypes (Fig. 5C).

Figure 4.

The Hub genes was high expressed in colorectal cancer. A, GEPIA was used to analyzed the expression level of AURKA, CCNB1, EXO1 and CCNA2 in normal and colorectal cancer tissues. Red indicates tumor tissue and blue indicates normal tissue. ^{*P <} 0.05. B, The expression level of AURKA, CCNA2, CCNB1 and EXO1 were obtained from GSE21510, GSE110224 and GSE74602 datasets. C, The AURKA, CCNB1, EXO1 and CCNA2 in different tumors were obtained from TIMER database, ^{*P <} 0.05, ^{**P <} 0.01, ^{***P <} 0. 001.

Figure 5.

The protein expression of Hub genes were up-regulated in colorectal cancer. A, HPA database was used to analyzed the protein expression of AURKA, CCNB1 and CCNA2 in normal and tumor tissue by immunohistochemistry. B, ROC curves of AURKA, CCNB1, EXO1 and CCNA2 were obtained from online tool Xiantao Academy. C, The differences of gene expression among four subtypes of colorectal cancer (POLE, MSI, CIN and GS) were investigated using R software.

3.3Genomic mutation of the Hub genes in CRC

We analyzed the gene mutations of AURKA, CCNB1, CCNA2 and EXO1 in CRC using the cBioPortal database. The results showed that AURKA had different mutation frequencies in different data sets, 13.64% of 22 cases (MSK, Cancer Discovery 202), 8.55% of 269 cases (TCGA, Nature 2012), 7.58% of 594 cases (TCGA, PanCancer Atlas), 4.26% of 47 cases (MSK, JCO Precis Oncol 2022), 3.82% of 471 cases (MSK, Gastroenterology 2020), 3.79% of 1134 cases (MSK, Cancer Cell 2018), 3.29% of 152 cases (MSK, Nat Commun 2022), 2.26% of 619 cases (DFCI, Cell Reports 2016) (Fig. 6A). Of the 3308 CRC patients in selected 8 data sets, 153 had a change in AURKA, and amplification is the most common type of AURKA mutation in CRC (Fig. 6B). There were 36 mutation sites in the AURKA gene (including 3 duplicate mutations in patients with multiple samples) (Fig. 6C). Furthermore, the results of survival analysis displayed that AURKA mutation had no effect on overall survival (P= 0.213, Fig. 6D), but impacted the progression free survival time of CRC patients (P= 0.0269, Fig. 6E).

Figure 6.

The Hub gene mutation affected the prognosis of patients with colorectal cancer. A–E, The mutation frequency (A), mutant type (B), mutation site (C) and the effect of mutation on overall survival (D) and progression free survival time (E) of AURKA gene in colorectal cancer were obtained from cBioPortal. F–J, The mutation frequency (F), mutant type (G), mutation site (H) and the effect of mutation on overall survival (I) and disease free (J) of CCNB1 gene in colorectal cancer were obtained from cBioPortal. K–O, The mutation frequency (K), mutant type (L), mutation site (M) and the effect of mutation on overall survival (N) and progression free survival time (O) of CCNA2 gene in colorectal cancer were obtained from cBioPortal. P–T, The mutation frequency (P), mutant type (Q), mutation site (R) and the effect of mutation on overall survival (S) and disease free (T) of EXO1 gene in colorectal cancer were obtained from cBioPortal.

The CCNB1 mutation frequencies in different data sets were as follows: 1.68% of 594 cases (TCGA, PanCancer Atlas), 1.12% of 269 cases (TCGA, Nature 2012) and 0.81% of 619 cases (DFCI, Cell Reports 2016) (Fig. 6F). Of the 2650 CRC patients in the selected 5 data sets, 18 had a change in CCNB1 (mutation rate was 0.7%), and deep deletion is the common type of CCNB1 mutation in CRC (Fig. 6G). There were 13 mutation sites in the CCNB1 gene (including 1 duplicate mutation in patients with multiple samples) (Fig. 6H). The results of survival analysis showed that the mutations of CCNB1 gene did not affect the overall survival (P= 0.502) and disease free time of CRC patients (P= 0.474) (Fig. 6I–J).

The CCNA2 mutation frequencies in different data sets were as follows: 1.52% of 594 cases (TCGA, PanCancer Atlas), 1.12% of 269 cases (TCGA, Nature 2012) and 0.81% of 619 cases (DFCI, Cell Reports 2016) (Fig. 6K). Of the 2650 CRC in the selected 5 databases, 17 had a change in CCNA2 (mutation rate was 0.6%) (Fig. 6L). There were 14 mutation sites in the CCNA2 gene (including 2 duplicate mutations in patients with multiple samples) (Fig. 6M). The mutations in the CCNA2 gene did not affect the overall survival (P= 0.655) and progression free survival time of CRC patients (P= 0.924, Fig. 6N–O).

The EXO1 mutation frequencies in different data sets were as follows: 2.36% of 594 cases (TCGA, PanCancer Atlas), 2.26% of 619 cases (DFCI, Cell Reports 2016), and 1.49% of 269 cases (TCGA, Nature 2012) (Fig. 6P). Of the 2650 CRC patients in the selected 5 databases, 32 had a change in EXO1, and deep deletion and amplification are the common type of EXO1 mutation in CRC (Fig. 6P and Q). There were 29 mutation sites in the EXO1 gene (including 1 duplicate mutation in patients with multiple samples) (Fig. 6R). Moreover, the EXO1 mutation had no effect on overall survival time (P= 0.150) and disease free time (P= 0.807) (Fig. 6S–T).

3.4Correlation analysis between the Hub genes expression and infiltrating immune cells

Tumor infiltrating lymphocytes affect the survival of patients with CRC. So we using the TIMER database analyzed the correlation of AURKA, CCNB1, CCNA2 and EXO1 with six kinds of infiltrating immune cells and tumor purity in CRC. The results displayed that the AURKA expression was correlation with tumor purity (r= 0.159, P= 1.33e-03), but no association with the immune cells (P> 0.05, Fig. 7A). In addition, compared with normal tissue, different copy states of AURKA have some effect on B cells, CD8+ T cells, neutrophils and dendritic cells immersion (Fig. 7B). The expression level of CCNB1 was correlation with B cells (r= 0.131, P= 8.16e-03), CD8+ T cells (r= 0.178, P= 3.17e-04), CD4+ T cells (r=-0.116, P= 2.01e-02), neutrophils (r= 0.194, P= 9.60e-5), and dendritic cells (r= 0.104, P= 3.75e-02) in CRC (Fig. 7C). Compared with normal tissue, different copy states of CCNB1 effected the level of B cells and CD8+ T cells immersion (Fig. 7D). The expression level of CCNA2 was correlation with B cells (r= 0.194, P= 8.78e-05), CD8+ T cells (r= 0.259, P= 1.19e-07), neutrophils (r= 0.237, P= 1.56e-06), and dendritic cells (r= 0.169, P= 6.76e-04) in COAD (Fig. 7E). The copy states of CCNA2 effected the level of B cells and CD8+ T cells infiltration (Fig. 7F). The expression level of EXO1 was correlation with B cells (r= 0.172, P= 5.10e-04), CD8+ T cells (r= 0.246, P= 5.36e-07), neutrophils (r= 0.343, P= 1.57e-12), and dendritic cells (r= 0.226, P= 4.83e-06) in CRC (Fig. 6G). The copy states of EXO1 have some effect on the level of B cells, CD8+ T cells and dendritic cells infiltration (Fig. 7H).

3.5Constructing PPI networks

The functional interaction between proteins is necessary for the molecular mechanism and metabolism of malignancy. Therefore, we constructed the PPI network of AURKA, CCNB1, CCNA2 and EXO1 protein using STRING. Cytoscape software was used to analyze the PPI network formed by each gene and Hub genes (Fig. 8). Their results were helpful to reveal the pathogenesis of CRC and to search for therapeutic targets and prognostic biomarkers.

Figure 7.

The expression of Hub gene affected the immune cell infiltration level in colorectal cancer. A, C, E, G, TIMER database was used to analyze the correlation between AURKA (A), CCNB1 (C), CCNA2 (E) and EXO1 (G) gene expression and immune cell infiltration level, respectively. B, D, F, H, TIMER database was used to analyze the relationship between the copy number of AURKA (B), CCNB1 (D), CCNA2 (F) and EXO1 (H) infiltration level, respectively. ^{*P <} 0.05, ^{**P <} 0.01, ^{***P <} 0. 001.

Figure 8.

The PPI network of Hub genes in colorectal cancer. A–D, The PPI network diagrams of AURKA (A), CCNB1 (B), CCNA2 (C) and EXO1 (D) were structured by STRING tool.

Figure 9.

The co-expression analysis of Hub genes. A–D, Linked Omics database was applied to analyzed the co-expressed genes of AURKA (A), CCNB1 (B), CCNA2 (C) and EXO1 (D). E, The correlations between Hub genes and their co-expressed genes were analyzed using Linked Omics database.

Figure 10.

The biological functions and signal pathways of Hub genes and co-expressed genes. A–D, The biological process (A), cellular component (B), molecular function (C), and signal pathways (D) of AURKA and co-expressed genes were analyzed by GSEA analysis. E–H, The biological process (E), cellular component (F), molecular function (G), and signal pathways (H) of CCNB1 and co-expressed genes were analyzed by GSEA analysis. I–L, The biological process (I), cellular component (J), molecular function (K), and signal pathways (L) of CCNA2 and co-expressed genes were analyzed by GSEA analysis. M–P, The biological process (M), cellular component (N), molecular function (O), and signal pathways (P) of EXO1 and co-expressed genes were analyzed by GSEA analysis.

3.6Analysis of co-expression genes associated with the Hub genes

The mRNA sequences of 379 patients with TCGA-COAD were analyzed by the functional module method. The volcanic map showed that the number of genes with positive correlation with AURKA, CCNB1, CCNA2 and EXO1 was higher than that of negative correlation (Fig. 9A–D). The 50 important genes that were positively correlated and 50 genes that were negatively correlated with AURKA, CCNB1 and EXO1 were shown in Fig. 8A–D. AURKA was strong positive with TPX2 (Pearson correlation = 0.85, P= 1.101e-105), UBE2C (Pearson correlation = 0.82, P= 3.014e-94) (Fig. 9E). CCNB1 was strong positive with DEPDC1B (Pearson correlation = 0.79, P= 1.006e-71), CCNA2 (Pearson correlation = 0.78, P= 3.162e-67) (Fig. 9E). CCNA2 showed strong positive with MAD2L1 (Pearson correlation = 0.91, P= 1.779e-147), PLK4 (Pearson correlation = 0.84, P= 2.510e-100) (Fig. 9E). EXO1 showed strong positive with DTL (Pearson correlation = 0.86, P= 7.280e-111), MCM10 (Pearson correlation = 0.81, P= 4.487e-90) (Fig. 9E).

3.7GSEA of the Hub genes and the co-expression genes

Furthermore, we analyzed the biological function of Hub genes and the co-expression genes by GSEA. AURKA co-expressed genes were mainly involved in the biological processing of chromosome segregation, DNA replication, cell cycle checkpoint and translational elongation, and associated with the cellular components including mitochondrial protein complex, mitochondrial inner membrane and ribosome (Fig. 10A and B). KEGG pathway analysis showed AURKA and co-expressed genes participated in Cell cycle, RNA transport and Ribosome biogenesis in eukaryotes (Fig. 10D). CCNB1, CCNA2 and their co-expressed genes were mainly involved in the biological processing of chromosome segregation, DNA replication and cell cycle checkpoint (Fig. 10E and I), and associated with cellular components including mitochondrial protein complex, mitochondrial inner membrane, mitochondrial matrix and ribosome (Fig. 10F and J). The molecular functions related to CCNB1 and CCNA2 included structural constituent of ribosome, catalytic activity-acting on DNA and RNA, extracellular matrix structural constituent, cyclin-dependent protein kinase activity (Fig. 10G and K). KEGG pathway analysis showed that CCNB1, CCNA2 and the co-expressed genes involved in cell cycle, ribosome, spliceosome, DNA replication and RNA transport (Fig. 10H and L). EXO1 and co-expressed genes were mainly located in chromosomal, mitochondrial and ribosome, and are mainly involved in physiological processes such as DNA replication, chromosome segregation, cell cycle checkpoint, DNA recombination, cytokinesis, chromatin assembly or disassembly, RNA localization and DNA-templated transcription, termination (Fig. 10M and N). The molecular functions of EXO1 and co-expressed genes include catalytic activity, acting on DNA and RNA, ATPase activity, nucleotidyltransferase activity and structural constituent of ribosome (Fig. 10O). KEGG path analysis shows that EXO1 and co-expressed genes participate in cell cycle pathways, DNA replication pathways, spliceosome pathways, RNA transport pathways, ribosome and lysosome pathways in CRC (Fig. 10P).