Gene expression data were downloaded from The Cancer Genome Atlas (TCGA) for colon cancer. TCGA colon adenocarcinoma (TCGA-COAD) is the project name (https://portal.gdc.cancer.gov/projects/TCGA-COAD). HTSeq-Counts data is included in this database, along with HTSeq-FPKM and HT-FPKM-UQ data. The RNA-seq samples of data type HTSeq-Counts were analyzed for this analysis. The total number of samples was 329. And the total number of 20,530 protein-coding genes were selected for the next analysis.

Clinical data on colon cancer were also acquired from the TCGA dataset, which included 447 colon cancer samples. There was a variety of clinical information available in the original dataset for each sample, but for this study, only the sample number and four cancer stages information were retrieved here. After following the correlation of gene expression data with the corresponding clinical data and extracting the missing information data, got 312 overlapping samples. All these samples (overlapping sample numbers) were categorized into five groups. They are healthy human tissue or normal tissue, stage I, stage II, stage III, and stage IV. There are 41 normal tissues, 45 stage I samples, 109 stage II samples, 80 stage III samples, and 37 stage IV samples. For further analysis, all the normal samples were combined with the individual stage. After that, the final samples would be 86, 150, 121, and 78 respectively. These are considered the final four differentially expressed (DE) sets.

For detecting differentially expressed genes (DE), this study mainly used gene expression analysis. To avoid the noises of raw data, data preprocessing is crucial to minimize noise. In addition, performing high-level analysis is a must for quality assessment and preprocessing of sequencing data. Counts Per Million (CPM) value should be calculated for each group of datasets due to filtering lowly expressed genes. Four stages and one normal sample dataset were used for identifying Differentially Expressed Genes. We also considered the TMM (Trimmed Mean of the M-value) algorithm (Robinson and Oshlack, 2010) because it is more effective for comparisons between samples as it does not count gene length or library size.

In this study, the edgeR package (Robinson et al., 2010) in R language was generated from Bioconductor and used to analyze and identify the differentially expressed genes (DEGs) list of the cancerous tissues for four DE sets. These analyses of DEGs were performed between control samples and CRC samples of four different stages. The p-values were calculated for all DEGs, which were further adjusted into false discovery rate (FDR) by using the Benjamini–Hochberg method. Fold Change (FC) value for each group was also calculated and only the genes with F D R   < 0.05 and l o g F C ≥ 1 were defined as differentially expressed (DE) genes.

An online tool, Vennplex (Cai et al., 2013), was used to identify not only all DEGs but also the intersection of DEGs between four DEG sets. It also showed the similarities and differences in expression changes between different groups. For this analysis, it is necessary that both gene symbol and log ⁡ 2 F C should be inputted. All common DEGs in these datasets were selected for further study.

For filtering the remarkable genes, histogram analysis was performed on both clinical and PPI data. The DEGs with l o g F C > 1.5 and in the PPI network, connections of edge confidence >0.75 were selected as the final DEGs data list as those are mostly connected to colon cancer.

In order to identify gene groups that conduct cancer related cellular functions, the functional interaction (FI) network was built. Because one of the motivations of this study is to analyze the functions of the group of genes and the FI network is unable to represent the relations between the coding genes. Four functional interaction networks were constructed using the 4 DE sets (filtered by histogram analysis). The nodes represented the genes, and the edge represented the functional associations and interaction associations between genes. The edge contains the implications of biological function, which provided the basis for the functional analysis. ReactomeFIViz (Wu et al., 2014), a tool in the Cytoscape software (Shannon et al., 2003), was utilized to analyze the functional interactions between all these DEGs.

In this study, the Markov Cluster Algorithm (MCL) was applied to each stage and whether the cluster has five nodes or above were connected, those clusters were selected to analyze the relationship between functional modules for each stage and adjacent stages. Finally, 25 clusters, 33 clusters, 27 clusters, and 29 clusters were chosen for stage I, stage II, stage III, and stage IV respectively. According to the connections of each cluster, seven networks were built, where four networks are for individual stage (stage I, stage II, stage III, and stage IV) and another three is between adjacent stages (between stage I and stage II, between stage II and stage III, between stage III and stage IV).

The main contribution of this study is the combination of seven networks with four groups (stage I, stage II, stage III, and stage IV) and seven connections (above seven networks). The Cytoscape software was used to make this combined network, called pathway interaction network, in order to analyze the relationship between functional modules. The final graphical view of the pathway interaction network was entitled as the functional evolutionary network. Where the correlation between four DEG sets was obtained and the staged genes involved in the pathway were used for further analysis.

In the functional evolutionary network analysis, the pathway enrichment analysis has been done. It was not necessary to perform a staged biological function analysis because the pathway enrichment provided a variety of the mixed results of pathway. In order to obtain results of functional analysis in stages, this study, therefore, considered using a method that represents pathways by a graph, which are enriched in different stages. This study analyzed two relationships between pathways and between adjacent stages, the first is to determine which pathways were significantly different between stages, and the second determination is to identify which pathways were associated with adjacent stages.

In this study, DAVID (Dennis et al., 2003) was used to perform the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses in order to identify the biological features of DEGs associated with biological functions as well as elucidate the functional annotation and pathway enrichment analysis to investigate the biological pathways of DEGs in each functional interaction module (criterion: FDR <0.05 and p-value <0.05).

Each biological system within a cell is controlled by proteins. Some proteins do their work on their own, but most of them rely on their interactions with others to perform their biological functions. To understand protein functions and the biological characteristics of the proteins, Protein-Protein Interaction is a must. PPI data became more available through high-throughput technologies and made it possible to visualize PPI data as networks, i.e., PPI networks.

In this study, all filtered DEGs were used to construct the PPI network based on the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (Szklarczyk et al., 2018). Here, the species was homo and the score criterion was 0.75. By using cytoHubba, a Cytoscape plug-in, finally selected the hub genes and constructed a network among those candidate genes.

This table was made with the overlapping samples (323), where various clinical information (patient’s age, gender, different stages information, and vital status) were concerned (Table 1).

A total of 312 samples were processed and normalized, where 271 colon cancer samples and 41 normal tissues were present. In the colon cancer samples, there were 45 samples for stage I, 109 for stage II, 80 for stage III, and 37 for stage IV. Finally, the DEGs were identified using adjusted and cut-off criteria. A total of 1,574, 1,607, 1,468, and 1,574 significant DEGs were found between healthy and colon cancer stages I, II, III, and IV, respectively.

There were a total of 2048 DE genes among 4 DE sets. 1,067 common DEGs were revealed through intersect function, including 1,185 upregulated and 863 downregulated DEGs in comparison with controls. And those common DEGs are defined as intersection DEGs. All the DEGs in the four sets were compared with each other by Vennplex (Figure 2) and it is found that downregulated genes were less compared with upregulated genes in each set. There were no contra-regulated DEGs which suggested that the process (molecular function) of colon cancer at all stages is identical.

Before the filtering process, there were 1,574, 1,607, 1,468, and 1,574 genes found for stages I to stage IV respectively. After making the histogram analysis based on the PPI network to filter these genes to get significant DE genes the number decreased in Table 2.

According to the DEGs detected in each stage of CRC, we constructed four FI networks at each stage, respectively as the figure shows. Where each node represents a differentially expressed gene. Edge represents the interactions between genes. Figures 3A–D are the FI network constructed by DEG-stage1, DEG-stage2, DEG-stage3, DEG-stage4, respectively. These DEGs normally formed some modules to conduct their cellular functions. That’s why MCL clustering was applied and 17, 15, 16, and 18 modules were found after the clustering at individual stages, respectively. These modules have a higher probability to enrich some functions of colon cancer.

Then the cluster interaction networks are constructed for each stage and adjacent stage. There are 7 cluster interaction networks in total. These networks are constructed based on the FI network of each cancer stage by using the MCL cluster algorithm. Where each cluster represents the functions of biological mechanisms and for the investigation of the complex molecular mechanisms of colon cancer these clusters may play important roles. In order to uncover the associations between genes and functions, the interaction between clusters is constructed, where each cluster is treated as a vertex within the interaction and the number of connections between genes in two corresponding clusters as the edge weight. According to the FI network used here, the weights denote the degree of association between clusters. The wider the edges the more closely biology functions between clusters.

Examining the functional modules of functional interaction networks at different stages of cancer can provide insight into the functional evolution of the disease and the strength of associations between genes at different stages. This is because these functional modules may contain overlapping genes, and analyzing them can provide valuable information about the progression of the disease.

Figure 4 shows the pathway interaction network between the four stages. Red nodes or “A”: DE set I, green nodes or “B”: DE set II, purple nodes or “C”: DE set III, blue nodes or “D”: DE set IV. Thick edges were constructed by overlapping gene counts. This pathway interaction network is a major contribution to the field of pathway enrichment analysis. The network was constructed after doing many analyses to get the significant pathways. This network is the result of combining 7 cluster networks. There are four connected networks and some other pathways among the four stages. In this study, the large network is more significant because all the stages are connected in different ways and the other three are connected to one, two, or three stages. This network is the final functional evolutionary network. In terms of the connections between functional modules between adjacent stages, this study constructs this functional evolution network. Nodes indicate the modules. Edges reflect the connections between modules at colon cancer stages. The overlapping genes between functional clusters at adjacent stages are represented by the network’s edges. The more overlapping genes there are between adjacent stages of functional modules, the thicker the edges. In the following sections, this network will be analyzed to investigate the molecular mechanism of colon cancer.

In this study, by analyzing the pathway interaction network and doing KEGG pathway enrichment analysis, 15 common pathways were found. They are: Neuroactive ligand-receptor interaction, Glutamatergic synapse, Circadian entrainment, Nicotine addiction, Retrograde endocannabinoid signaling, cAMP signaling pathway, Dopaminergic synapse, Adrenergic signaling in cardiomyocytes, Amphetamine addiction, Bile secretion, Wnt signaling pathway, Dilated cardiomyopathy, Long-term potentiation, Cardiac muscle contraction, Melanogenesis. Among them, cAMP signaling pathway controls essential physiological activities such as metabolism, secretion, calcium homeostasis, muscular contraction, cell destiny, and gene transcription (Cho-Chung et al., 2002). The neuroactive ligand-receptor interaction pathway is a group of receptors and ligands on the plasma membrane that are linked to intracellular and extracellular signaling pathways and this pathway is associated with prostate cancer (He et al., 2018). The APC mutant colon cancer cells remain dependent on Wnt and suppress the production of secreted Wnt antagonists epigenetically (He et al., 2005). Ke yang et al. demonstrated that the IWP inhibitors inhibit the WNT signaling pathway in colon cancer cells by disrupting the WNT ligand (Yang et al., 2016).

For analyzing the different pathways, KEGG pathway analysis was done and about five significantly different pathways were found. Pathways in cancer, Serotonergic synapse, Bile secretion, Hypertrophic cardiomyopathy (HCM), and Dilated cardiomyopathy are significantly different at all stages. The findings pathway namely Serotonergic synapse pathway is a neurotransmitter which is widely distributed in the vertebrate central nervous system, and it serves as a target for many physiologic regulations like modulators of gene transcription, steroids and neurotrophic factors. Ionotropic GluRs and 5-HT3 receptors mediate rapid synaptic transmission via serotonergic fibers making direct synaptic connections with GABAergic neurons (Maejima et al., 2013). There are active transport systems within hepatocytes and cholangiocytes. Hepatocytes secret their bile by secreting conjugate bilirubin, bile salts, cholesterol, phospholipids, and water into their canaliculi (Hundt et al., 2022). A primary myocardial disorder associated with the autosomal dominant pattern of inheritance is hypertrophic cardiomyopathy (HCM), which can be distinguished from other cardiovascular disorders primarily by the presence of myocyte hypertrophy, fibrillar disarrays, and interstitial fibrosis as histological characteristics. Dilated cardiomyopathy (DCM) is a heart muscle illness described by dilation and impaired contraction of the left or both ventricles that is responsible for progressive heart failure and sudden cardiac death from ventricular arrhythmia.

To identify the significant pathways, heatmap analysis was done based on the p-values of each stage. Table 3 shows a total of 15 pathways with the p-values of each stage.

To investigate the significant biological functions of DEGs, heatmap analysis was done on each stage (Figure 5) and among those pathways, some are significantly different such as neuroactive ligand-receptor interaction, Glutamatergic synapse, Circadian entrainment, and Nicotine addiction have low p-value (from Table 3) and a large number of genes.

From Table 3; Figure 5, it is clear that the p-values of neuroactive ligand-receptor interaction, Glutamatergic synapse, Circadian entrainment, and Nicotine addiction change and increase from the initial stage to later stage and Figure 5 also shows that the color changes from light yellow to marron and it is said for heatmap analysis is that bright color indicates high activity and dark color is vice versa.

Fang et al. and Liu et al. has shown that the neuroactive ligand-receptor interaction signaling pathway is linked to bladder cancer and renal cell carcinoma progression (Fang et al., 2013; Liu et al., 2015). It is believed that glutamatergic synapse pathways play a crucial role in a large variety of normal physiological functions due to their links to many other neurotransmitter pathways and neurodevelopmental disorders and injuries are strongly associated with glutamate dysfunction (Glutamatergic Synapse Pathway, 2022). Circadian entrainment includes retinal sensitivity as well as circadian variations in the retina that contribute to the regulation of retinal diseases and these circadian disorders are related to entrainment deficits (Golombek and Rosenstein, 2010). It is shown in various studies that Nicotine addiction increased approximately 50% chances of colon cancer because smoking has been associated with adenomatous polyps (Slattery et al., 1997). The genes encoding these pathways may be important in the molecular development of these pathways since they are strongly associated with function.

A Protein-Protein Interaction network (PPI) is a visual framework for better understanding protein functional organization. When the iteration between genes was found by above analysis, the PPI network was also constructed to check their relationship because the PPI network is the widely used network to see whether there is a strong relationship or not. In the PPI network, all the significant DE genes were investigated and a PPI network was built. Among these DEGs, some DEGs also appeared in functional interaction network and some DEGs are also known colon cancer related genes, such as CXCL11, ADH1B, PYY, SLC17A7, and so on. The network involved 794 nodes, and 7,153 edges (Figure 6).

In this study, the top 10 hub genes (INS, SNAP25, GRIA2, SST, GCG, PVALB, SLC17A7, SLC32A1, SLC17A6, and NPY) were identified according to the degree of connectivity of the DEGs and arranged it in a descending order to find the highest degree of connectivity (Table 4). Figure 7 showed the network between the top 10 hub genes.