Identification of target genes and prognostic evaluation for colorectal cancer using integrated bioinformatics analysis

Abstract

The underlying molecular mechanisms of colorectal cancer (CRC) has attracted great attention from the scholarly community. The aim of our study is to identify pivotal genes related to the pathogenesis and prognosis of CRC. We integrated five microarray datasets from Gene Expression Omnibus (GEO) database. The differentially expressed genes (DEGs) were analyzed with the limma package. DAVID and OmicShare tools were used for Gene Ontology (GO) and KEGG enrichment analysis. The protein–protein interaction (PPI) network of DEGs was constructed. The prognostic analysis of hub genes was performed through Kaplan Meier-plotter. Finally, potential drugs were predicted in the CMap database. Through five microarray datasets, a total of 90 DEGs were detected including 54 up-regulated and 36 down-regulated genes. Biological process analysis showed DEGs were mainly enriched in positive regulation of neutrophil chemotaxis, chemokine-mediated signaling pathway, and bicarbonate transport. Signaling pathway analysis indicated that DEGs played a vital in proximal tubule bicarbonate reclamation, cell cycle and progesterone-mediated oocyte maturation. The GEPIA database confirmed that overexpression levels of hub genes were significantly associated with better survival of patients. Finally, the 20 most significant small molecules were obtained based on the CMap database. Our study has identified novel candidate biomarkers, pathways, and kinases associated with CRC prognosis.

KEYWORDS:

• Bioinformatics analysis
• colorectal cancer
• candidate small molecules
• differentially expressed genes
• hub genes

Introduction

In the world, Colorectal Cancer (CRC) is common cancer and a major public health problem. Overall, the incidence of colorectal cancer (CRC) ranks third and the cancer-related mortality ranks second, with approximately 1.88 million new cases and 915,880 deaths reported worldwide in 2021.

Colonoscopy – today’s gold standard due to its capacity to screen CRC and remove precancerous lesions – is invasive, inconvenient for the patients, and costly. The disadvantages inherent to colonoscopy makes it inappropriate for use as a primary screening tool (Benard et al. 2018). However, it has proven to be effective in reducing CRC incidence and mortality. Colorectal cancer often manifests as early lymph node invasion or distant metastases, with patients often diagnosed only when they possess advanced disease, a phenomenon that inevitably leads to dismal clinical outcomes. CRC patients with relapse respond inadequately to the second or third line of treatment making their prognosis generally poor. The 5-year survival rate for metastatic CRC remains only 12% (Siegel et al. 2019).

In addition, environmental, socioeconomic, and lifestyle in individual may be regarded as risk factors of developing CRC (Keum and Giovannucci 2019). As other malignancies, CRC is considered a polyphase disease in different stages of tumor initiation, progression, and metastasis, which are accompanied by many change, including genetic mutations, epigenetic changes, cellular contexts, and environmental influences (Boussios et al. 2019). It has shown in several studies that the occurrence and development of CRC involves multiple genes and cellular signal transduction pathways, which are used for early tumor screening. Therefore, the demand for identifying patients at a higher or lower risk of CRC death through sensitive and specific biological markers is essential, not only to better understand the molecular and cellular processes of tumorigenesis but also for more effective early diagnosis, appropriate therapeutic strategies, and improved prognosis.

Over the recent years, bioinformatics serve as an important tool for a range of life science applications, which could be used for early diagnosis of the tumor, tumor grade, and prognostic judgment (Canzoneri et al. 2019). Various bioinformatics databases, such as the GEO database, TCGA database, provide opportunities for data mining for gene expression profiles of cancer. The public database Gene Expression Omnibus (GEO) is a huge repository of high-throughput gene expression data and hybridization arrays, chips, microarrays (Jiang and Liu 2015). GEO database has provided a powerful platform for bioinformatics to explore novel biomarkers for cancer diagnosis, treatment, and prognosis analysis. For this study, we performed an integrated bioinformatics analysis to identify specific genes differentially expressed, discover the hub genes; and build protein to protein interaction (PPI) networks of hub genes by analyzing various bioinformatic online datasets to reveal the potential pathogenesis of colorectal cancer.

Materials and methods

Microarray data

The five gene expression profile matrix files and platform TXT files (GSE110223, GSE1110224, GSE113513, GSE84984, and GSE73360) were downloaded from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/). These RNA profiles were based on GPL96 platform (Affymetrix Human Genome U133A Array), GPL570 platform(Affymetrix Human Genome U133 Plus 2.0 Array), GPL15207 platform (Affymetrix Human Gene Expression Array), and GPL17586 platform (Affymetrix Human Transcriptome Array 2.0), containing a total of 100 colorectal cancer samples and 79 normal samples. The GSE110223 profile included 13 tumor samples and 13 normal samples, GSE110224 profile contained 17 tumor samples and 17 normal samples, GSE113513 profile provided 14 tumor samples and 14 normal samples, GSE84984 profile comprised from 9 tumor samples and 6 normal samples, and GSE73360 profile consisted of 37 tumor samples and 19 normal samples. The characteristics of these datasets are shown in Table 1.

Table 1. Details of the GEO colorectal cancer data.

References	Sample	GEO	Platform	Tumor	Normal
Vlachavas EI et al. (2019)	colorectal adenocarcinomas	GSE110223	GPL96	13	13
Vlachavas EI et al.(2019)	colorectal adenocarcinomas	GSE110224	GPL570	17	17
Peng J et al.	colorectal cancer	GSE113513	GPL15207	14	14
Condorelli DF et al.(2018)	colorectal adenocarcinoma	GSE84984	GPL17586	9	6
Condorelli DF et al. (2018)	colorectal adenocarcinoma	GSE73360	GPL17586	37	19

GEO, gene expression omnibus.

Screening of DEGs

The analysis was carried out by using GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/), a convenient online analysis tool that is based on R programming language and can be used to perform comparisons on GEO datasets (Yao and Liu 2018). Adjusted P < 0.05 and |log fold change (FC)| > 1 were chosen as the cutoff value for statistical analysis of each dataset. The adjusted P-value (adj. P) was applied to help correct false positives. The heat map and volcano plot were drawn with the ‘ggplots’ package in R 3.6.2. The upregulated and downregulated gene lists were saved as Excel files for subsequent integration analysis.

Gene ontology and Kegg analysis of DEGs

To depict the potential biological functions of the overlapping DEGs, we performed Gene Ontology (GO, http://www.geneontology.org) enrichment analysis covering three domains: biological process (BP), molecular function (MF), and cellular component (CC) (Hinderer et al. 2019). OmicShare website, a free online platform for Go enrichment analysis (http://omicshare.com/tools), Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.ad.jp/kegg/), an extensive database of genomes, biological pathways, diseases, drugs, and chemical substances, are often used to explore functions and metabolic pathways involved in the overlapping DEGs (Du et al. 2014). DAVID (http://david.ncifcrf.gov), an online bioinformatics resource, provides highly integrated gene-annotation and comprehensive data coverage (Jiao et al. 2012). It has covered a wealth of biological data and relevant analysis tools for the biological function annotation information of genes or proteins. P < 0.05 was considered as the cut-off criterion with a significant difference. Downloading the TXT file of the GO and KEGG pathway enrichment analysis for subsequent analysis. A visual network analysis of the KEGG analysis was performed using Cytoscape 3.7.2 software.

Pathway enrichment analysis of DEGs

The Kyoto Encyclopedia of Genes and Genomes (KEGG) database mapped genes and their respective pathways together (Du et al. 2014). ClueGO (Mlecnik et al. 2018), an app that selects Gene Ontology terms and pathways for a list of genes or proteins and visualizes them into functional networks, runs on Cytoscape (Otasek et al. 2019). The plug-in CluePedia (Mlecnik et al. 2018) is another Cytoscape tool for finding markers that may be related to pathways. Additionally, better results were achieved by combing ClueGO with CluePedia. In this study, DEG functions were studied using ClueGO and CluePedia toolkits through KEGG pathway enrichment analysis. A P value of <.05 and kappa coefficient of 0.4 were considered as threshold values.

Protein–protein interaction (PPI) network construction and module analysis

The STRING, a database of a customizable protein–protein network, and functional characterization of user-uploaded gene/measurement sets (STRING, https://string-db.org/) (Szklarczyk et al. 2021) aims to predict functional associations between proteins for a large number of organisms. The DEGs were imported to STRING to assess the interactive relationships, and a combined score > 0.4 was considered significant. Subsequently, the PPI network was visualized and further analyzed using Cytoscape software (www.cytoscape.org/) (Otasek et al. 2019). The MCODE (Molecular Complex Detection) plugin in was utilized to explore significant modules in the PPI network. The degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and max depth = 100 were set as cutoff criteria for advanced options. Using the Cytoscape plugin cytoHubba (Chin et al. 2014), we determined the hub proteins in the PPI network.

Comparison of the hub genes expression level and validation of hub genes

GEPIA (http://gepia.cancer-pku.cn/detail.php) web server organizes cancer and normal gene expression data for analysis and provides online interactive analysis following standard methods (Tang et al. 2017). It provides customizable functions such as tumor and normal differential expression analysis. Then the relationship for hub genes was indicated in boxplots. We also analyzed the protein expression of the hub genes between CRC and normal tissues using the human protein atlas (HPA, www.proteinatlas.org) database, which is an online tool widely used for determining the protein level of genes in clinical samples.

Survival analysis of hub genes

Kaplan Meier-plotter (KM plotter, http://kmplot.com/analysis/) could assess the effect of 54675 genes on survival using 10,461 cancer samples, including 5143 breast, 1816 ovarian, 2437 lung and 1,065 gastric cancer patients with the average follow-up time was 69, 40, 49, and 33 months. The relapse-free and overall survival information was based on GEO (Affymetrix microarrays only), EGA, and TCGA database. The hazard ratio (HR) with 95% confidence intervals and log rank P value were calculated and displayed on the plot.

Identification of small molecules

We queried the Connectivity Map (CMap, http: //www. broadinstitute.org/cmap/) to detect the candidate small molecule drugs based on the gene signature of CRC. CMap is a gene transcription-expression profiling database that links disease, genes, and drugs by similar or opposite gene expression profiles. First, these overlapping DEGs were divided into upregulated and down-regulated groups. Then the probesets from each group were utilized to query the CMap database (Subramanian et al. 2017). Finally, the enrichment scores ranging from −1 to +1 were calculated, which represented the similarity. A positive connectivity value (closer to +1) demonstrated that a small molecule can induce the gene expression of CRC cells, whereas a negative connectivity value (closer to −1) demonstrated that a small molecule can imitate the status of normal cells.

Results

Identification of DEGs in CRC

In the aggregate, 90 colorectal adenocarcinoma samples and 69 matched normal mucosa samples from 5 datasets 2019 were analyzed. The major characters of the qualified studies were shown in Table 1. After the colorectal cancer chip expression datasets GSE110223, GSE110224, GSE113513, GSE84984, and GSE73360 were normalized, we then applied the GEO2R online analysis tool to screen the DEGs, using adjust P value < 0.05 and logFC ≤ −1 or logFC ≥ 1 as the cut-off criteria. We captured 90 DEGs including 54 upregulated DEGs and 36 downregulated DEGs. The GSE110223 dataset contained 727 differential genes, including 285 upregulated genes and 442 downregulated genes. The GSE110224 dataset contained 1086 differential genes, including 469 upregulated genes and 617 downregulated genes. The GSE113513 dataset contained 2866 differential genes, including 1155 upregulated expression genes and 1711 downregulated expression genes. The GSE73360 dataset contained 4621 differential genes, including 3000 upregulated expression genes and 1621 downregulated expression genes. In addition, the GSE84984 dataset contained 8723 differential genes, including 4184 upregulated genes and 4539 downregulated genes. The volcano plot showed the up-regulated and down-regulated genes in each dataset (Figure 1A). The Venn diagram from the five datasets summarizng results was shown in Figure 1(B and C) shows the heat map of the top 100 genes grouped by clusters.

Figure 1. (A)Volcano plot of gene expression profile data between clcorectal cancer and normal tissues in each dataset. Red dots: significantly up-regulated genes in CRC; Blue dots: significantly down-regulated genes in CRC; Black dots: non-differentially expressed genes. P < 0.05 and |log2 FC| >1 were considered as significant.(B) The Venn diagram of 172 overlapping DEGs among five datasets.(C) Cluster heat map of the top 100 DEGs in five GEO databases. Red indicates relative upregulation of gene expression; blue indicates the relative down-regulation of gene expression; gray indicates no significant change in gene expression; and white indicates that the signal intensity is not high enough to detect.

GO functionenrichment analysis of DEGs

Based on GO analysis, the uDEGs exhibited distinct differences in regulation of neutrophil chemotaxis, chemokine-mediated signaling pathway, and extracellular matrix disassembly in biological processes (BPs) (Figure 2A, Table 2); the dDEGs are mainly in bicarbonate transport, regulation of intracellular pH, and one-carbon metabolic process (Figure 2B, Table 2). For molecular function (MF), the uDEGs are enriched in CXCR chemokine receptor binding, growth factor activity, and chemokine activity (Figure 2A, Table 2); and the dDEGs are enriched in carbonate dehydratase activity, chloride channel activity, and zinc ion binding (Figure 2B, Table 2). Cellular component (CC) analysis revealed that the uDEGs are concentrated in the extracellular region, extracellular space, and extracellular exosome (Figure 2A, Table 2); and the dDEGs are concentrated in the basolateral plasma membrane, plasma membrane, and extracellular exosome (Figure 2B, Table 2).

Figure 2. (A) GO analysis of upregulated DEGs. (B) GO analysis of downregulated DEGs.(C) KEGG pathway of DEGs. (D) Function annotation of DEGs. Each node is a Kyoto Encyclopedia of Genes and Genomes pathway item. The node size reflects pathway significance (fdr): the smaller the fdr value, the larger the node size is. Edge between nodes reflects shared or common genes: the wider the edge, the larger the overlap is. Different node colors represent different functional groups.

Table 2. Gene ontology analysis of DEGs related to colorectal cancer.

Expressed	Category	Term	Count	%	PValue	FDR
Upregulated	GOTERM_BP	positive regulation of neutrophil chemotaxis	3	15	2.48E-04	0.299118335
	GOTERM_BP	chemokine-mediated signaling pathway	3	15	0.002581357	3.079205442
	GOTERM_BP	extracellular matrix disassembly	3	15	0.002951142	3.51311397
	GOTERM_BP	cell proliferation	4	20	0.006571075	7.667655489
	GOTERM_BP	inflammatory response	4	20	0.007235163	8.411761691
	GOTERM_BP	regulation of epithelial cell proliferation	2	10	0.012791867	14.42581271
	GOTERM_BP	response to lipopolysaccharide	3	15	0.01309069	14.73872445
	GOTERM_BP	positive regulation of protein oligomerization	2	10	0.014908761	16.62000317
	GOTERM_BP	extracellular matrix organization	3	15	0.018342044	20.06917335
	GOTERM_BP	visual perception	3	15	0.019231499	20.94113799
	GOTERM_BP	cell development	2	10	0.042041175	40.53142925
	GOTERM_BP	regulation of MAPK cascade	2	10	0.043070501	41.30004314
	GOTERM_BP	activation of protein kinase activity	2	10	0.047177366	44.2768663
	GOTERM_BP	positive regulation of pathway-restricted SMAD protein phosphorylation	2	10	0.050246585	46.41042667
	GOTERM_BP	SMAD protein signal transduction	2	10	0.064446555	55.3403133
	GOTERM_BP	collagen catabolic process	2	10	0.066458678	56.48880317
	GOTERM_BP	cell chemotaxis	2	10	0.067463206	57.0519747
	GOTERM_BP	immune response	3	15	0.073675265	60.38871792
	GOTERM_BP	cell adhesion	3	15	0.085538938	66.10935876
	GOTERM_BP	proteolysis	3	15	0.098959723	71.6611433
	GOTERM_MF	CXCR chemokine receptor binding	3	15	4.30E-05	0.041411974
	GOTERM_MF	growth factor activity	4	20	7.51E-04	0.720843063
	GOTERM_MF	chemokine activity	3	15	0.00136762	1.309464185
	GOTERM_MF	cytokine activity	3	15	0.01645129	14.76576966
	GOTERM_MF	extracellular matrix binding	2	10	0.028876728	24.58911097
	GOTERM_MF	serine-type endopeptidase activity	3	15	0.032823596	27.48977165
	GOTERM_MF	calcium ion binding	4	20	0.044533029	35.51644964
	GOTERM_MF	transforming growth factor beta receptor binding	2	10	0.046252454	36.62546176
	GOTERM_MF	endopeptidase activity	2	10	0.05909052	44.38035175
	GOTERM_CC	extracellular region	10	50	1.31E-05	0.012150714
	GOTERM_CC	extracellular space	9	45	3.15E-05	0.029319651
	GOTERM_CC	extracellular exosome	8	40	0.01883732	16.21966573
	GOTERM_CC	proteinaceous extracellular matrix	3	15	0.031251675	25.58174481
Downregulated	GOTERM_BP	bicarbonate transport	6	30	4.18E-10	4.36E-07
	GOTERM_BP	regulation of intracellular pH	4	20	4.97E-06	0.005189023
	GOTERM_BP	one-carbon metabolic process	3	15	3.65E-04	0.379865389
	GOTERM_BP	chloride transmembrane transport	3	15	0.003461968	3.555115106
	GOTERM_BP	transport	4	20	0.004043789	4.141244566
	GOTERM_BP	oxalate transport	2	10	0.010434482	10.37058701
	GOTERM_BP	ethanol oxidation	2	10	0.011377993	11.25858918
	GOTERM_BP	sulfate transmembrane transport	2	10	0.011377993	11.25858918
	GOTERM_BP	excretion	2	10	0.034693462	30.82695793
	GOTERM_BP	cellular response to cAMP	2	10	0.048434249	40.44123883
	GOTERM_BP	regulation of membrane potential	2	10	0.069147783	52.66534589
	GOTERM_MF	carbonate dehydratase activity	3	15	9.70E-05	0.088695576
	GOTERM_MF	chloride channel activity	3	15	0.001487056	1.352297764
	GOTERM_MF	zinc ion binding	6	30	0.006323352	5.638425427
	GOTERM_MF	arylesterase activity	2	10	0.006381638	5.689052262
	GOTERM_MF	alcohol dehydrogenase activity, zinc-dependent	2	10	0.006381638	5.689052262
	GOTERM_MF	secondary active sulfate transmembrane transporter activity	2	10	0.011670268	10.18332919
	GOTERM_MF	sulfate transmembrane transporter activity	2	10	0.011670268	10.18332919
	GOTERM_MF	oxalate transmembrane transporter activity	2	10	0.011670268	10.18332919
	GOTERM_MF	bicarbonate transmembrane transporter activity	2	10	0.015882019	13.62495221
	GOTERM_MF	Goxidoreductase activity	3	15	0.018864211	15.99030148
	GOTERM_MF	transporter activity	3	15	0.019220202	16.26876758
	GOTERM_MF	inorganic anion exchanger activity	2	10	0.022167875	18.54312013
	GOTERM_MF	anion:anion antiporter activity	2	10	0.02425472	20.11984241
	GOTERM_CC	basolateral plasma membrane	4	20	8.17E-04	0.736183624
	GOTERM_CC	plasma membrane	11	55	0.004140857	3.680563952
	GOTERM_CC	extracellular exosome	8	40	0.01883732	15.79056824
	GOTERM_CC	apical plasma membrane	3	15	0.036338915	28.43206115
	GOTERM_CC	brush border membrane	2	10	0.051878416	38.21084849
	GOTERM_CC	integral component of plasma membrane	5	25	0.054847417	39.93764033
	GOTERM_CC	apical part of cell	2	10	0.075399694	50.76037166

KEGG signal pathway analysis of DEGs

A KEGG pathway analysis of the integrated DEGs is shown in Table 3 and Figure 2(C). The significantly enriched pathways, by DEGs, included the Proximal tubule bicarbonate reclamation pathway, Cell cycle pathway, Progesterone-mediated oocyte maturation pathway, Nitrogen metabolism pathway, Bile secretion pathway, Pancreatic secretion pathway, and Oocyte meiosis pathway. Additionally, cyclin-dependent kinase 1 (CDK1) was found to participate in the cell cycle pathway, Progesterone-mediated oocyte maturation pathway, Oocyte meiosis pathway and p53 signaling pathway, while in the cell cycle pathway, cell division cycle 6 (CDC6) was enriched and in the Pancreatic secretion pathway, mitotic arrest deficient-like 1 (yeast) (MAD2L1) was enriched, respectively. The network diagram was drawn using Cytoscape software (Figure 2D).

Table 3. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of integrated DEGs.

Pathway	ID	Count	PValue	Genes
Proximal tubule bicarbonate reclamation	hsa04964	4	6.12E-04	CA4, CA2, SLC4A4, PCK1
Cell cycle	hsa04110	6	0.002125	CCNB1, CDK1, CDC6, MAD2L1, CCNA2, CDC25B
Progesterone-mediated oocyte maturation	hsa04914	5	0.003794	CCNB1, CDK1, MAD2L1, CCNA2, CDC25B
Nitrogen metabolism	hsa00910	3	0.006826	CA4, CA2, CA1
Bile secretion	hsa04976	4	0.014269	AQP8, CA2, SLC4A4, ABCG2
Pancreatic secretion	hsa04972	4	0.031235	SLC26A3, CLCA4, CA2, SLC4A4
Oocyte meiosis	hsa04114	4	0.048697	CCNB1, CDK1, MAD2L1, AURKA
Biosynthesis of antibiotics	hsa01130	5	0.071226	SHMT2, SQLE, PSAT1, PPAT, PCK1
p53 signaling pathway	hsa04115	3	0.087884	CCNB1, CDK1, RRM2
PPAR signaling pathway	hsa03320	3	0.087884	SCD, MMP1, PCK1

PPI network construction and module selection.

Considering the critical role of protein interactions in protein function, we used the STRING database and Cytoscape software to generate PPI network once we had identified the 90 DEGs. The results showed that there were dense regions in PPI, that is, genes closely related to colorectal cancer (HUB genes) modules.

A total of 77 nodes and 300 edges were selected to plot the PPI network, which consisted of 54 up-regulated genes and 36 down-regulated genes (Figure 3A). Subsequently, a pivotal module of 19 genes (CDK1,TPX2, DLGAP5, UBE2C and so on) was identified with the degree≧2 as the cut-off value by using MCODE (Figure 3C). Another important module of 9 genes including only up-regulated genes was also identified (Figure 3D). The top 10 HUB genes were identified by cytoHubba (Top 10 genes ranked in Degree) and the PPI net work was established as well (Figure 3B), the details of which were showed in Table 4.

Figure 3. PPI and MCODE analyses of DEGs. (A) Protein–protein interaction network of 90 DEGs. (B) The protein-protein interaction(PPI) network of top 10 hub genes.(C) A signifcant module, containing 19 up-regulated proteins, was selected from PPI network. (C) Another module selected from PPI network. For (A, C, D), red nodes are up-regulated proteins, and blue nodes are down-regulated proteins. Circles represent genes, lines represent interactions between gene-encoded proteins and line colors represent evidence of interactions between proteins.

Table 4. The degree values of the top 10 hub genes.

Gene symbol	Gene description	Log FC	Degree
AURKA	aurora kinase A	1.21	23
CDK1	cyclin dependent kinase 1	1.36	23
MKI67	marker of proliferation Ki-67	1.15	22
TPX2	microtubule nucleation factor	1.96	22
CCNB1	cyclin B1	1.58	22
UBE2C	ubiquitin conjugating enzyme E2 C	1.48	21
CCNA2	cyclin A2	1.34	21
CDC6	cell division cycle 6	1.14	21
MAD2L1	mitotic arrest deficient-like 1 (yeast)	1.32	20
DLGAP5	DLG associated protein 5	1.16	19

FC, fold change.

KEGG analysis of the top 2 most important modules were conducted. The results showed that an apparent enrichment of genes in module 1 was observed in the Cell cycle, Progesterone-mediated oocyte maturation, p53 signaling pathway, and Oocyte meiosis, while genes in Module 2 are mainly enriched in ECM-receptor interaction (Table 5).

Table 5. KEGG enrichment of genes in the top 2 modules.

Module	Pathway	Count	PValue	Genes
Module 1	Cell cycle	5	1.23E-05	CCNB1, CDC6, CDK1, MAD2L1, CCNA2
	Progesterone-mediated oocyte maturation	4	1.59E-04	CCNB1, CDK1, MAD2L1, CCNA2
	p53 signaling pathway	3	0.0031471	CCNB1, CDK1, RRM2
	Oocyte meiosis	3	0.008485555	CDK1, MAD2L1, AURKA
Module 2	ECM-receptor interaction	3	0.003187	CD44, COL1A2, COL1A1
	Rheumatoid arthritis	2	0.08622	MMP3, MMP1
	Protein digestion and absorption	2	0.08622	COL1A2, COL1A1

Expression levels and prognostic analysis of hub genes

GEPIA is an interactive online server for exploring large data sets from the TCGA and GTEx projects. To confirm the reliability of the top 10 identified hub genes, we used GEPIA to verify the correlation between them, the results obtained from the GEPIA demonstrated that the levels of CDK1, MKI67, CDC6, and MAD2L1 were significantly differentially expressed between tumor (include colonic adenocarcinoma, COAD, and rectum adenocarcinoma, READ) and normal tissues (Figure 5A), the overexpression of which in CRC tissues were also confirmed in the above 5 GEO datasets (Figure 4). Since gene expression does not always correlate with the level of a protein encoded by the corresponding gene, we further measured the protein level of CDK1, MKI67, CDC6, and MAD2L1 in clinical CRC tissues from the HPA database. The results of immunohistochemical indicated that the protein expression level of CDK1, MKI67, CDC6, and MAD2L1 exhibited a significant difference between CRC tissues and normal tissues (Figure 5B).

Figure 4. The expression level of CDK1, MKI67, CDC6 and MAD2L1 between CRC and normal tissues in five datasets.

Figure 5. (A) The expression level of CDK1, MKI67, CDC6 and MAD2L1 hub genes between CRC and normal in COAD and READ in GEPIA database. (B) Representative immunohistochemistry staining results reveal the protein level expression of CDK1, MKI67, CDC6 and MAD2L1 in colorectal cancer. (C) The top 20 candidate small molecules targeting the gene expression of CRC.(D) The prognostic value of hub genes according to the Kaplan Meier-plotter online database. Abbreviations: COAD:colon adenocarcinoma; READ:rectal adenocarcinoma; CRC: colorectal cancer.

To assess the prognostic value of 10 hub genes, we used Kaplan Meier-plotter online database for prognostic analysis. The results showed that high expression levels of CDK1, MKI67, CDC6, and MAD2L1 may correlate with a better prognosis in patients with CRC (P < 0.05) (Figure 5D).

Identification of related active small molecules

To screen out candidate small molecule drugs, the CMap database was utilized to analyze consistent differently expressed probesets between CRC tissues and normal tissues. Table 6 and Figure 5(C) showed the predicted small molecules that could inhibit CRC-associated gene expression. Among these small molecules, DL-thiorphan (enrichment score = −0.977) and blebbistatin (enrichment score = −0.828) showed a higher negative correlation between the gene expression of CRC and the small molecules, which could be the most promising small molecules to reverse the tumoral status of CRC. The prediction of small molecule drugs aims to exploit the potential drugs targeting CRC and make existing drugs fully utilized. However, further studies were urgently needed to validate the above results.

Table 6. List of the 20 most significant small molecule drugs that can reverse the tumoral status of colorectal cancer.

Rank	cmap name	n	Enrichment	p
11	DL-thiorphan	2	−0.977	0.00107
53	sanguinarine	2	−0.926	0.01131
6	repaglinide	4	−0.901	0.00016
8	ellipticine	4	−0.884	0.00042
113	piperlongumine	2	−0.866	0.03603
44	milrinone	3	−0.841	0.00809
46	ursodeoxycholic acid	3	−0.838	0.00843
162	blebbistatin	2	−0.828	0.05799
50	trazodone	3	−0.825	0.01076
193	5248896	2	−0.802	0.077
196	butein	2	−0.801	0.07746
218	menadione	2	−0.789	0.08722
72	ronidazole	3	−0.785	0.02029
30	acepromazine	4	−0.78	0.00487
77	alsterpaullone	3	−0.78	0.02177
34	Prestwick-1084	4	−0.771	0.00559
13	meticrane	5	−0.767	0.0013
37	thioguanosine	4	−0.766	0.00607
40	daunorubicin	4	−0.764	0.00635
90	triflusal	3	−0.761	0.0278

Discussion

The failure to early screening and diagnosis in CRC patients results in poor overall survival rates and prognosis. Therefore, sensitive and specific biomarkers for CRC are urgently needed to be detected.

In the present study, a total of 90 DEGs, including 54 upregulated and 36 downregulated genes, were identified in CRC tissues compared to in-paired adjacent normal tissues. Our results showed that AURKA, CDK1, MKI67, TPX2, CCNB1 and so on were the hub proteins with higher degrees in the PPI network. Furthermore, the proximal tubule bicarbonate reclamation, the cell cycle, Progesterone-mediated oocyte maturation, and Nitrogen metabolism signaling pathways were significantly enriched by DEGs. In addition, CDK1 was enriched in the cell cycle pathway, progesterone-mediated oocyte maturation pathway, Oocyte meiosis pathway, and p53 signaling pathway, CDC6 was enriched in the cell cycle pathway and MAD2L1 was enriched in the pancreatic secretion signaling pathway. We performed a prognostic analysis of these 10 hub genes using the Kaplan Meier-plotter online database. The results demonstrated that the prognosis of colorectal cancer patients with high expression levels of CDK1, MKI67, CDC6, and MAD2L1 is better than that with low expression of theirs’. Of the identified differentially expressed hub genes, we found that they were generally involved in the regulation of the cell cycle and were closely associated with the prognosis of colorectal cancer. AURKA is a cell cycle-regulated kinase and is involved in colon carcinogenesis through promoting genomic instability. In addition, AURKA could connect deregulated genes in the Wnt and Ras-MAPK signal pathways, which promots the proliferation, invasion, and migration of CRC cells (Jacobsen et al. 2018). One previous study indicated that CDK1 upregulation was linked to BRAFV600E mutation in a subset of colorectal cancer patients (Barras et al. 2017). Zhang et al (2018b) suggested that BRAFV600E colorectal cancer cell lines expressing CDK1 were sensitized to apoptosis upon siRNA knockdown or small-molecule inhibition with RO-3306 (CDK1 inhibitor) or dinaciclib (CDK1, 2, 5, 9 inhibitors). Combination of RO-3306 or dinaciclib with cobimetinib (MEK inhibitor) cooperatively enhanced apoptosis and reduced clonogenic survival versus monotherapy. MKI67 is a nuclear non-histone protein required for cell proliferation and is widely used as a biomarker to quantify the proportion of cleavage cells to grade tumors (Miller et al. 2018). Ki67 expression is a marker for cells that have enter the cell cycle and is also proven to be correlated with many associated genes in CRC. Deletion of SNRPA1 can down-regulate the expression of MKI67 as an upstream molecule in CRC cells (Zeng et al. 2019). Through its cooperation with MYC, TPX2 is a driver gene which promotes colorectal adenoma to carcinoma (Takahashi et al. 2015). Cancer cells overexpressing CCNB1 grow and proliferate more rapidly in CRC (Fang et al. 2014). UBE2C is an E2 ubiquitin-conjugating enzyme that catalyzes the degradation of APC/C substrates through polyubiquitination in the regulation of the cell cycle. A growing body of evidence suggests that UBE2C could be connected to colorectal cancer progression and to drug resistance in breast cancer (Zhang et al. 2018). A specific cell-cycle-associated gene associated with poor prognosis for CRC patients is CCNA2 (Gan et al. 2018). CDC6 functions as a regulator at the early steps of DNA replication and is reported as a promising target for overcoming oxaliplatin resistance for CRC (Cai et al. 2019). MAD2L1 participates in triggering the spindle assembly checkpoint (SAC) and its hetero-tetramer complex, mitotic arrest deficient-like 1 (MAD1L1), to promote the cell cycle checkpoint signal (Zhong et al. 2015). DLGAP5, DLG associated protein 5, encoded DLG7 (disc large homolog 7), which is a microtubule-associated protein (MAP), and has a crucial role during spindle assembly, kinetochore fibers (K-fibers) stabilization, and chromosomal segregation during mitosis. DLGAP5 overexpression has been found in many malignancies such as CRC, breast cancer, hepatocellular carcinoma, urinary bladder cancer, meningioma, and adrenocortical cancer. These previous studies have already described an association between DLGAP5 deregulation and genomic instability (Branchi et al. 2019; Sinha et al. 2019). Thus, DLGAP5 overexpression may contribute to carcinogenesis. Meanwhile, we performed a meta-analysis to pool the 4 hub genes expression levels in different datasets. Our meta-analysis also revealed CDK1, MKI67, CDC6, and MAD2L1 were significantly elevated in CRC tissues as compared to normal tissue (Figure 6).

Figure 6. Meta-analysis of the expression amount of CDK1, MKI67, CDC6 and MAD2L1 among the five datasets. Abbreviations: CRC: colorectal cancer.

In addition, based on the overlapping DEGs, we obtained several potential small molecule drugs for the treatment of CRC from the CMap database. These small molecules had the potential to alter the gene expression of CRC, thereby controlling the progression of tumors. For example, repaglinide (enrichment score = − 0.901) has cytotoxic effects against hepatic, breast, and cervical carcinoma cells in previous researches and has been clinically used for the treatment of glioblastoma multiforme. In the meantime, it is still unknown whether it is effective or safe on CRC (El Sharkawi et al. 2014; Xiao et al. 2017). Therefore, further studies are urgently demanded to confirm the huge potential of these candidate small molecules in treating CRC. Summary, the top 10 hub genes identified in the PPI network are closely associated with tumorigenesis and tumor progression, which suggests that they can be prognostic factors or therapeutic targets in CRC.

Ontologies of gene products are an important way to describe and integrate biological processes, cellular components, and molecular functions across species. In our study, Based on GO and KEGG enrichment analysis, DEGs were especially enriched in processes related to extracellular matrix disassembly, positive regulation of neutrophil chemotaxis, cell proliferation, and growth factor expression. Genes which regulate growth factors secretion cause tumorigenesis, migration, invasion, and metastasis of cancer cells by stimulating the release of large amounts of TGF-β, HGF, and FGF (Erdogan and Webb 2017). Extracellular matrix (ECM) regulates the development and maintains tissue homeostasis. The loss of ECM's homeostasis is generally regarded as defining the transitional phases of cancer progression and metastasis (Rigoglio et al. 2020). Alterations in tumor ECM actively promote tumor growth, survival, invasion, and metastasis, as well as angiogenesis and immune function (Varol and Sagi 2018). As well, DEGs are mapped on major colorectal cancer pathways, including cell cycle signaling pathways, chemokine-mediated signaling pathways, and MAPK cascades. Signal transduction pathways are involved in the regulation of cell growth, survival, adhesion, movement, and invasion. Our study showed that cell cycle signaling pathways were significantly enriched by DEGs. Cell cycle control is crucial for the maintenance of genomic integrity. Many tumors rely on bypassing physiological checkpoints in cell cycle control as a central component of uncontrolled cellular proliferation. Targeting these cell-cycle dysregulation mechanisms is an attractive method to inhibit tumorigenesis. It is well established that cell cycle control has remained essentially constant throughout evolution, from simple unicellular eukaryotes to complex mammals, including humans. A crucial transition in the cell cycle is driven by the dynamic interplay between cyclins and their associated cyclin-dependent kinases (CDKs) (Xu and McArthur 2016). CDKs are serine/threonine kinases that regulate the cell cycle by interacting with specific cell-cycle–regulatory cyclins. Among the CDKs, CDK1 is necessary for G2–M transition as well as the progression of G1 and the transition from G1 to S (Haneke et al. 2020; Kalous et al. 2020). The non-selective CDK inhibitor dinaciclib inhibits tumor growth in multiple types of cancer in preclinical studies by arresting cell-cycle progression and triggering caspase-8, −9, and related apoptosis (Buzzetti et al. 2021; Howard et al. 2021), cyclins are abnormally expressed in colorectal cancer cells, resulting in increased proliferation, greater viability, and longer survival times (Shen et al. 2019). The results of this study are largely consistent with findings reported for prostate cancer, urinary bladder, and hepatocellular carcinoma (Chang et al. 2015; Wang et al. 2018; Ji et al. 2021). Chemokines and their receptors also function in the development and function of multiple tissues. Chemokine expression has prognostic significance in several cancers apparently due to chemokine-mediated growth and metastatic spread. Endocytosis is required for CXC chemokine receptor type 4 (CXCR4)-mediated Akt activation and antiapoptotic signaling. This mechanism whereby compartmentalized chemokine-mediated Akt signaling from endosomes suppresses the cancer-related process known as anoikis. CXCR4 receptors may be targets of this signaling pathway for the treatment of metastatic cancer (English et al. 2018). In other words, studying these pathways can help elucidate the mechanisms of proliferation and invasion of colorectal cancer and predict its progression.

In this study, we also performed a module analysis on the constructed PPI network and selected the two most important modules, and we then performed a KEGG pathway analysis of the genes in these modules. All the 2 identified modules with hub nodes were involved in the important biological function and tumor pathways. The results showed that Module 1 genes are primarily associated with the cell cycle, progesterone-mediated oocyte maturation, p53 signaling pathway and Oocyte meiosis, and Module 2 genes are mainly involved in the ECM-receptor interaction. This indicates that cancer development and metastasis are closely related to these signal transduction pathways, and that blocking the cell cycle and blocking these pathways can provide a beneficial treatment for cancer.

However, there were a few limitations in our study. First, the analysis of the 5 datasets from GEO will result in partial bias, and too few samples in each dataset will not lead to new findings. Although we perform normalization of each dataset individually, the heterogeneity of scientific dissimilarity within every single study cannot be disentangled, result in the meta-analysis is still unclear. Second, most of the included gene expression studies come from Western countries and may pose geographical bias. Therefore further studies based on the information from ongoing trials are required to verify potential biomarkers of CRC.

Conclusion

In conclusion, this study provides some perspectives on potential biomarkers linked to the prognosis of CRC patients. In addition, this study highlights the importance of PPI and modules analysis as a potent framework to gain insight into the key hub nodes underlying the prognosis of CRC and to recognize potential biomarkers of CRC. We also found DL-thiorphan was the most promising small molecule to reverse the gene expression of CRC. These findings could provide new sights about the future genomic individualized treatment of CRC and survival prediction.

Acknowledgements

We express our cordial thanks to all those who participated in this study. Rui Chen, Shoucheng Ma and Fei Su: Conception, design, data analysis, and manuscript writing, Hui Qiao and Lina Wang: Data collection, Quanlin Guan: Guide, supervise, and review.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability

The following information was supplied regarding data availability: The data is available at NCBI GEO: GSE110223, GSE1110224, GSE113513, GSE84984 and GSE73360. The GEO database is available at the following websites: https://www.ncbi.nlm.nih.gov/geoprofiles/

Additional information

Funding

This research was supported in part by the Hospital foundation of the First Hospital of Lanzhou University [grant number ldyyyn2019-84] and Natural Science Foundation of Gansu Province youth Fund [grant number 21JR7RA386].