Two microarray profiles of colon cancer samples including recurrence information (E-GEOD-39,582 and E-GEOD-33,113) were downloaded from ArrayExpress (http://www.ebi.ac.uk/arrayexpress/). E-GEOD-39,582 included 566 samples, based on the platform of AFFY HG-U133_Plus_2, which were divided into the training dataset and the validation dataset for the weighted correlation network analysis (WGCNA) network construction. There were a total of 90 colon cancer samples in E-GEOD-33,113, with the recurrence status and clinical information of the samples, which was used as the validation dataset for mining vital module associated with recurrence.

Samples in the profiles were classified using Consensus Clustering in R (http://www.bioconductor.org/packages/release/bioc/html/ConsensusClusterPlus.html) to identify their sub-types. The parameters set in this analysis are presented in Table I. Following the classification of samples, the association between colon cancer samples and corresponding clinical performance [survival, sex, cghdata, chemotherapy, mismatch repair (Mmr) status, tumor node metastasis (TNM) stage, and Tumor location] were determiend using the χ² test. P<0.05 was considered to indicate a statistically significant difference.

Data in the profiles were initially normalized using the robust multi-array analysis (RMA) method in Affy package (http://www.bioconductor.org/packages/release/bioc/html/affy.html) and were transformed using a log₂ transformation. Probes were converted into gene symbols and the average expression value was used as the only value of the gene with multiple corresponding probes. Following that, the cv method in the genefilter package was applied to filter the genes with significant variation. Genes with coefficient of variation >0.4 were recognized as the candidate genes.

The limma package in R (http://www.bioconductor.org/packages/release/bioc/html/limma.html) was used to identify the DEGs from the candidate genes. P<0.01 and log fold change (FC) >1.5 were set as the cut-off criteria.

The WGCNA method (12) was used to construct the co-expression network of the genes in the test samples of colon cancer. The interaction coefficient between genes was calculated using the following formula:

Where x_i and x_j are the expression vectors of gene i and j, respectively. The Pearson coefficient of these two vectors was cor, which was transited into the interaction coefficient using S_ij. This transition intended to give more weight to the strong connections and reduce the importance of the weak connections in the predicted co-expression network in order to improve the reliability of the co-expression network.

The connection coefficient will be transformed into a weighted coefficient W_ij using the following formula:

Subsequently, the co-expression network would be constructed based on the W matrix, followed by module mining. The reliability of the minded modules was verified using the verify sample set E-GEOD-33113. The topological properties were also confirmed. Modules with module significance <0.05 were identified to be recurrence-associated modules.

The average expression value of genes in the significant modules in each sample was calculated. Subsequently, the samples were ranked based on the expression level of modules. According to 1/4 and 3/4 value of expression level, the modules were divided into high expression level, median expression level, and low expression level. Finally, the Kaplan-Meier (KM) curves of recurrence-associated modules expression level and recurrence status were drawn. The significance of the different expression level modules was compared, and P<0.05 was considered to indicate a statistically significant difference. The identified vital module was verified using the independent dataset E-GEOD-33113.

In order to determine the functions of the genes in the recurrence-associated modules, genes were subjected to Database for Annotation, Visualization and Integrated Discovery (13) for function and pathway enrichment. P<0.05 was the threshold used for the significant terms.

Survival analysis was performed on the genes in the recurrence-associated modules. The degree (k) of a gene and the significance (P) between each gene and sample survival time were also calculated using a Cox regression. Subsequently, the interaction coefficient (coef) between k and -log₁₀(P) was computed to identify the hub genes associated with survival. Finally, selected genes were determined in the verifed samples for an association with survival. Recurrence associated genes were also subject to sample classification.

From the 19,846 irredundant genes in E-GEOD-39582, 6,600 were the variation genes. Consensus Clustering analysis revealed that cumulative distribution function (CDF) was at a high level when there were 3 sub-types, accompanied with a satisfied classification effect (Fig. 1). Therefore, the 556 samples in E-GEOD-39582 were divided into 3 subgroups: G1, G2, and G3 (Fig. 2). There were significant survival differences among these 3 types of samples, and as depicted in Fig. 3, G1 exhibited the highest survival status, whereas G3 exhibited the lowest survival status. Additionally, the c² test determined that all clinical data, with the exception of sex, were significantly different among these 3 groups (Table II).

From the 434 DEGs, 76 were the DEGs between G1 and G2, 390 were the DEGs between G1 and G3 and 63 were DEGs between G2 and G3. A total of 2 DEGs were identified in all 3 groups (Fig. 4). These two common DEGs were stress-associated endoplasmic reticulum protein family member 2 (SERP2) and long non-coding RNA-0219 (LINC0219).

The co-expression network of the training dataset was divided into 4 modules (Fig. 5), which were verified with the validation dataset. By calculating the correlation coefficient, the connections between the genes in each module and survival status of the samples were identified. It was determined that the Brown module had the highest module significance and there were various survival-associated genes (hub genes) in this module (Fig. 6). These connections were also verified in the validation dataset.

All 431 genes in the Brown module were subject to a clustering analysis and it was revealed that there were 3 types of samples (Fig. 7). As presented in Fig. 7, certain G2 type samples were also observed in with G1 and G2 type samples, suggesting it may have a medium role for the connection between the G1 and G3 status. Additionally, the recurrence interval of G3 samples was short, accompanied with downregulated genes.

The top 10 genes were collagen type VI α 3 chain (COL6A3), EGF containing fibulin like extracellular matrix protein 2 (EFEMP2), fibrillin 1 (FBN1), follistatin-like 1 (FSTL1), glycosyltransferase 8 domain containing 2 (GLT8D2), heart development protein with EGF like domains 1 (HEG1), RAB31, member RAS oncogene family (RAB31), secreted protein acidic and cysteine rich (SPARC), SPARC/osteonectin, cwcv and kazal like domains proteoglycan 1 (SPOCK1) and TIMP metallopeptidase inhibitor 2 (TIMP2) in the Brown module with higher degrees were associated with survival (Table III). Genes in the Brown module were primarily enriched in tumor recurrence-associated functions and pathways, including cell adhesion, biological adhesion, ECM organism, the ECM-receptor interaction and the focal adhesion pathways (Table IV).

As presented in the KM curves (Fig. 8), there was a significant difference in terms of recurrence status (P=1.5×10⁻⁶) among the samples with high, median, and low gene expression level of the Brown module. High expression levels inidcated higher incidence of recurrence, whereas low expression levels show lower recurrence incidence. This association between gene expression level and recurrence incidence was validated in the dataset of E-GEOD-33113 (P=2.8×10⁻²).