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

2.1 Data 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.

OPEN IN VIEWER

The differentially expressed genes (DEGs) were filtered using a threshold of adjusted P values 0.05 and an absolute log₂FC (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.3 GO 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.4 Protein-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.5 Hub 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. 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.6 Protein 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.7 Gene 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.8 Immune 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.9 Co-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.10 Statistical 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. 0.05 was considered statistically significant.

3.1 Identification 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).

OPEN IN VIEWER

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).

OPEN IN VIEWER

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 ( 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 ( 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 ( 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 expres-sion 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).

OPEN IN VIEWER

OPEN IN VIEWER

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).

OPEN IN VIEWER

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).

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.5 Constructing 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.

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

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.7 GSEA 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).