In the present study, gene expression profile data were downloaded from The Cancer Genome Atlas database (TCGA; http://cancergenome.nih.gov/) (18). The database contained eight colon adenocarcinoma (COAD) (19,20) samples and eight normal samples (https://portal.gdc.cancer.gov), generated in multiple studies.

The downloaded sample files were merged into a gene expression matrix. Genes with missing expression values were removed, and the expression values of the remaining genes were log₂ transformed. Following preprocessing, the expression matrix with rows and columns contained 14,662 genes and 16 samples.

The downloaded sample data contained data for samples from different batches. The batch differences were removed by batch normalization using the ComBat procedure implemented in the SVA R package (version 1.28.0) (21). Subsequently, the normalize.quantiles.robust function in the preprocess Core package (http://bioconductor.org/packages/release/bioc/html/preprocessCore) was applied to perform unified normalization. Normalized data were used for subsequent analyses, including screening for differences in gene expression and network construction.

Levels of differential gene expression were calculated using the limma R package in Bioconductor (version 3.22.7) (22) via calculation of the log₂ fold change (FC) value and the P-value of each gene. Greater |log₂FC| values of genes indicated greater differences in the expression of these genes compared with the normal group and the COAD group. In general, genes with log₂FC >1 and P<0.05 were considered to be upregulated, and genes with log₂FC <-1 and P<0.05 were considered to be downregulated. Furthermore, genes were separated into three categories, according to the log₂FC value calculated with the limma R package (22). Firstly, genes for which the difference in the expression levels between the disease group and the control group was not significant (−0.5 <log₂FC <0.5) were excluded from the present study. The remaining genes were separated into two categories for subsequent analysis: i) Genes with log₂FC ≥0.5; and ii) genes with log₂FC ≤-0.5.

The WGCNA package (version 1.64–1) (http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGNA) provides a comprehensive collection of functions for conducting weighted correlation network analysis (23). Instead of describing the correlation structure between thousands of genes and a sample trait, WGCNA analysis focuses on the association between the sample trait and a few, usually <10, modules (24). In the WGCNA algorithm, the elements in the co-expression matrix of the genes are no longer the correlation coefficients of the genes, but rather the weighted value of the correlation coefficients. The criteria for the weighted value are such that the connections between the genes contained in each gene network can follow the scale-free law, in which p(i) is inversely proportional to iⁿ, where I is the node degree (connectivity) and p(i) is the probability that a node has n links (degree i). In practical applications, the network is made to an approximate scale-free distribution by selecting the weighting coefficients such that log(i) and log[p(i)] are negatively correlated, and the correlation coefficient should be at least 0.8.

The specific construction process of WGCNA networks includes three steps. In step one, the co-expression matrix of genes is defined. The gene correlation matrix S=[Smn] is constructed based on the correlation coefficient Smn=|cor (m, n)| between the gene m and the gene n. In the second step, adjacency functions are defined. In the WGCNA algorithm, for any gene pair, the adjacency coefficient a_mn is used as a measure of inter-gene correlation: a_mn=power (S_mn, β)=|S_mn|^β. In step three, the parameter β of the adjacency function is determined according to the scale-free network principle.

After satisfying the above conditions, the network can be constructed and divided into modules. Linking modules to known features involves two aspects. The first is the calculation of the characteristic values of the module, followed by the calculation of the correlation coefficient between the feature vector of the module and the feature of interest. The second involves grouped phenotypic data (such as disease status), for which a P-value for each gene for differential expression between each group (e.g. disease and normal groups) is calculated using the t-test, and the log₁₀(P-value) is defined as gene significance (GS). The module significance (MS) of each module is defined as the mean value of the GS of the genes contained in the module. The MS values are compared. In general, if the MS value of a module is significantly higher than that of other modules, this module may be related to the existence of the disease.

In the present study, the WGCNA package (24) was used to construct the network and hierarchical clustering tree.

The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.niaid.nih.gov) can help investigators in the functional interpretation of large lists of genes (25). In the present study, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for important genes in the identified module were performed using the DAVID Bioinformatics Resources (25). P<0.05 was set as the criterion for identifying overrepresented GO terms and pathways. In addition, the connection of the genes in the selected module was visualized using Cytoscape software (version 3.1.0) (26).

Total RNA was extracted from six pairs (female to male ratio, 1:2; mean age, 71; age range, 52–83) of colon cancer tissues and corresponding noncancerous colon tissues, which were collected from the Department of Gastrointestinal Colorectal and Anal Surgery, China-Japan Union Hospital of Jilin University (Changchun, China) using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA), according to the manufacturer's instructions. All participants underwent no other treatment before resection. All tissues were collected between January and April 2018. Harvested tissues were immediately frozen in liquid nitrogen and stored at −80°C prior to RNA extraction. The study was approved by the institutional ethical committee of China-Japan Union Hospital of Jilin University and informed consent was obtained from every patient. The quality and quantity of RNA samples were evaluated using an Infinite M100 PRO microplate reader (Tecan Group, Ltd., Mannedorf, Switzerland). RNA was reverse transcribed to cDNA using a PrimeScript™ RT Master Mix (Takara Biotechnology Co., Ltd., Dalian, China). All cDNA was amplified using the following primer sets: GAPDH forward, 5′-TGACAACTTTGGTATCGTGGAAGG-3′ and reverse, 5′-AGGCAGGGATGATGTTCTGGAGAG-3′); carbamoyl-phosphate synthetase 2, aspartate transcarbamylase, and dihydroorotase (CAD) forward, 5′-CCATGCACTAGACAGCCAAGA-3′ and reverse, 5′-CGGCTCAGTGTGGATACGAC-3′; transmembrane protein 147 (TMEM147) forward, 5′-ACACGCTATGATCTGTACCACA-3′ and reverse, 5′-CAGAGGTGGACGAAGGTCTC-3′; and σ-non-opioid intracellular receptor 1 (SIGMAR1) forward, 5′-CGAAGAGATAGCGCAGTTGG-3′ and reverse, 5′-TCCACGATCAGACGAGAGAAG-3′. GAPDH was used as a reference gene for normalization. Power SYBR-Green PCR Master (Thermo Fisher Scientific, Inc.) was used for qPCR, according to the manufacturer's instructions. Each reaction was performed in a final volume of 20 µl, containing 8 µl of cDNA, 1 µl of each primer and 10 µl 2X SYBR Premix EX Taq (Thermo Fisher Scientific, Inc.). RT-qPCR was performed on a Viia7 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) using the following thermocycling conditions: Denaturation at 50°C for 3 min and 95°C for 3 min; followed by 40 cycles of 95°C for 10 sec and 60°C for 30 sec. Gene expression levels were quantified using the 2^−ΔΔCq method.

All data are presented as the mean ± standard error of the mean and were analyzed using SPSS 22.0 software (SPSS, Inc., Chicago, IL, USA). Differences between colon cancer samples and control samples were determined using the Student's t-test. All experiments were repeated three times. P<0.05 was considered to indicate a statistically significant difference.

Following limma analysis, 6,134 genes that were not significantly different between the COAD samples and controls were removed, and 8,528 genes were retained. The 8,528 genes included 4,388 upregulated genes with log₂FC ≥0.5 in group A and 4,140 downregulated genes with log₂FC ≤-0.5 in group B. These genes were used for the subsequent analyses.

The WGCNA package was used to perform gene cluster analysis on the genes (log₂FC ≥0.5) (Fig. 1A). In a cluster dendrogram, height is a measure of dissimilarity according to the topological overlap matrix (23). By selecting different height cutoff values, the gene outliers were screened out and the number of genes was controlled. The genes within the first branch of the hierarchical cluster tree (Fig. 1A) were considered a research target for the follow-up analysis.

WGCNA requires the network to follow the scale-free distribution. As shown in the left panel of Fig. 1B, when β=13, the network satisfied the scale-free characteristic for the first time, and the vertical axis value exceeded 0.8 (the location of the green line in the figure), which is a prerequisite for building a WGCNA network. The figure on the right of Fig. 1B depicts the average connectivity of the network.

After determining whether the network obeyed the scale-free distribution, the hierarchical clustering tree was constructed and the gene modules were identified. As shown in Fig. 1C, branches of the hierarchical cluster tree defined nine modules with assigned colors. A total of two methods were used to examine the association between each module and colon cancer. The first was the MS value. The gene significance of the genes in each module was calculated. The MS was defined as the mean value of GS. A higher MS value for a module indicated that module had a stronger correlation with the disease. The second method was an MS correlation analysis (27). As shown in Fig. 1D and Table I, the turquoise module displayed the strongest correlation with the disease.

KEGG pathway enrichment analysis was conducted for the genes in the turquoise module (Table II). Genes in this module were mainly related to ‘RNA polymerase’ and ‘purine metabolism’. The top 30 genes with high connectivity were selected from the turquoise module for GO functional annotation. These genes were mainly enriched in the biological processes ‘translation’ and ‘gene expression’, in the cellular component ‘ribonucleoprotein complex’ and in the molecular function ‘structural constituent of ribosome’ (Table III). The associations of these 30 genes are presented in Fig. 2. The top 30 genes were all upregulated with log₂FC >1, including POTE ankyrin domain family member E, SIGMAR1, protein arginine methyltransferase 1, galactokinase 1, TMEM147 and CAD.

WGCNA analysis was also performed for genes in group B with log₂FC ≤-0.5. However, the analysis did not meet the condition of constructing a scale-free network and it was not possible to identify significant gene modules for group B.

To validate the gene expression changes identified by the aforementioned bioinformatics analysis, RT-qPCR was used to evaluate a number of potentially critical genes in color cancer, including CAD, TMEM147 and SIGMAR1, which exhibited a high degree in the constructed network. As shown in Fig. 3, the expression of TMEM147 was significantly higher in the colon cancer tissues than in the control tissues. However, no significant differences in the expression of CAD and SIGMAR1 were evident between colon cancer tissues and control tissues, perhaps reflecting the small sample size.