The raw expression data of GSE39582 (11) and GSE14333 (12) were retrieved from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/), both based on the platform of GPL570 Affymetrix Human Genome U133 Plus 2.0 Array. We used the Affy package in R to transform the CEL files of the tumor samples into an expression matrix (13). To improve the data quality, we used the k-nearest neighbors algorithm (k-NN) from the impute package in R to impute the missing expression data (14). Meanwhile, the robust multiarray average algorithm (RMA) was utilized to adjust the data for potential batch effects and for background correcting (15). Prior to WGCNA analysis, we filtered out the probes that were absent in all samples. The probe information was then transformed into the official gene symbols using Bioconductor in R. If multiple probes were applied to detect the same mRNA, the average value of the probes was used. The genes that were not differentially expressed between samples had to be excluded from WGCNA, as two genes without notable variance in expression between patients will be highly correlated. We chose the 75% most varying genes to construct the weighted gene coexpression networks. Specifically, the median absolute deviation (MAD) was used as a robust measure of variability.

In addition, the level three RNA-sequencing data of both colon adenocarcinoma (COAD) and rectum adenocarcinoma (READ) patients were downloaded from The Cancer Genome Atlas data portal (TCGA; http://cancergenome.nih.gov/). In contrast to ChIP-sequencing data, we used the voom function in package limma to normalize the TCGA data and create an expression matrix for samples for which the detailed clinical data are available (16,17). The voom method estimates the mean variance of the log counts and generates a precision weight for each observation. Thus, the WGCNA workflows originally developed for microarray analysis can be used on the RNA-seq data. Further preprocessing steps included the removal of control samples and the genes with zero counts in more than 80% of samples. As mentioned before, genes that are not differentially expressed between samples must be excluded; thus, we chose the top 12,000 genes with the highest MAD for the network building. Fig. 1 depicts a flow chart for the bioinformatic analysis.

The R package ‘WGCNA’ was used in our study to construct a gene coexpression network (10). After data collection and normalization, it is crucial that outliers be excluded. However, it was difficult to distinguish outlying samples in a dendrogram when the number of samples was large. To solve this problem, we used the standardized connectivity (Z. K) method recommended by WGCNA authors with the default threshold, Z. K score £2. After filtering out the outlying samples, expression data were tested to determine whether the samples and genes were good using the integrated function in the WGCNA package.

After filtering out the outliers and bad samples in the dataset, the next step of WGCNA is to build a scale-free network. In a scale-free network, several nodes, which are called hub nodes, are highly connected to other nodes in the network (18). In our study, we use the unsigned coexpression measure, which means that the positive correlation and negative correlation are equal. We constructed the gene coexpression network using the following steps.

First, we need a soft thresholding power β to which coexpression similarity is raised to calculate adjacency. By raising the absolute value of the correlation to a power β≥1 (soft thresholding), the weighted gene coexpression network construction emphasizes high correlations at the expense of low correlations. To determine the best soft threshold power, scale independence and average connectivity degree of modules with different power values were calculated by the gradient method. We selected the power β to ensure that the coexpression network was a ‘scale-free’ network, which was biologically close to reality. Moreover, to minimize the effects of noise and spurious associations, we subsequently constructed the Topology Overlap Matrix (TOM) from the adjacency matrix and calculated the corresponding dissimilarity (1-TOM), as well (19).

In the same way, the second coexpression network was built from TCGA data.

The traditional static tree cut method exhibits suboptimal performance on complicated dendrograms. In WGCNA, we tend to use the dynamic tree cut method by hierarchically clustering genes using the dissimilarity matrix (1-TOM) (20). The minimal size of a module was set as 30, and modules with high similarity were identified by clustering and then merged together with a height cut-off of 0.25. To determine whether the modules are reproducible, we tested the preservation of all modules with an independent gene expression dataset, GSE14333. We used the module preservation function (number of permutations set to 100) integrated in the WGCNA package to calculate the Z summary score of each module (21). In this method, a Z summary <2 indicates that the modules have no preservation, a Z summary of 2–10 indicates low to moderate preservation, and a Z summary >10 means that the module is strongly preserved.

The module eigengenes (MEs), which were measured by principal component analysis (PCA), were generated for each coexpressed module along with the module identification procedure.

We used two methods to identify the module of interest. First, we performed a module-trait relationship (MTR) analysis by calculating the correlation between module eigengenes and external clinical parameters, especially the anatomical site of the tumor. Having the module-trait relationships heatmap drawn, it was easy for us to identify which module related to the tumor location most.

Second, we measured gene significance based on the correlation of a gene expression profile with a sample trait and following module significance as an average absolute gene significance measure for all genes in a given module. Then, we plotted the barplot of the module significance for all modules detected. The highest module means it had the strongest correlation with the clinical trait.

In the key module, the hub genes were those that showed the most connections in the network. We called this property module membership, also known as eigengene-based connectivity kME, and in this instance, we used the default threshold value of 0.8. In addition to the module membership, the hub genes we need should also have a relatively higher gene significance; in this instance, we used the cut-off value as 0.4 (TGCA data set to 0.3). Combing both characteristics, we easily filtered out our hub gene in the module.

We applied Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) to detect the difference in expression levels of each hub gene between tumor and normal tissues in both the COAD and READ datasets from TCGA (22). To further validate our method, correlation plots between hub genes were generated by GEPIA, as well.

Twenty non-selected CRC samples were applied to perform qPCR to validate coexpression of PLAG1 like zinc finger 2 (PLAGL2) and protein O-fucosyltransferase 1 (POFUT1). These experimental samples were collected at the Sir Run Run Shaw Hospital of Zhejiang University between January 2004 and December 2006. After total RNA was isolated from tumor specimens using Trizol reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA), RNA was quantified by NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Inc.) and reverse transcribed using RNeasy Mini Kit (Takara, Kyoto, Japan) according to the manufacturer's protocols. Quantitative real-time PCR was executed with SYBR Green Master Mix (Takara). Relative expression levels were calculated with 2^−ΔΔCq formula (23). Expression of mRNA was standardized according to β-actin. The primers used were as follows: β-actin_fwd, ACTCTTCCAGCCTTCCTTCC and β-actin_rev, CGTCATACTCCTGCTTGCTG; PLAGL2_fwd, GAGTCAAGTGAAGTGCCAATGT and PLAGL2_rev, TGAGGGCAGCTATATGGTCTC; POFU-T1_fwd, AACCAGGCCGATCACTTCTTG and POFUT1_rev, GTTGGTGAAAGGAGGCTTGTG. The primers were designed on online tools (https://www.genscript.com/tools/real-time-pcr-tagman-primer-design-tool) and these were synthesized by Shanghai Generay Biotech Co. Ltd. (Shanghai, China).

We performed survival analysis for hub genes using the GSE39582 dataset because of its complete overall survival information. Kaplan-Meier analysis and log-rank test were performed to evaluate the association between hub gene expression and patient survival in left- and right-sided CRC, respectively. This procedure utilized the survival package in R (24), and the Kaplan-Meier survival curves with the at-risk table were drawn using the survminer package (25).

To identify the possible pathway through which hub genes may play a part in the development of CRC, the expression data from GSE14333 was also used to perform Gene Set Enrichment Analysis (GSEA). The expression data of 290 cases were uniformly divided into two groups according to each hub gene's expression value.

We used the GSEA-p 2.0 software to conduct the enrichment analysis (26). For configuration, ‘c2.cp.kegg.v6.2.symbols.gmt’ from the Molecular signatures database (MSigDB) 3.0 (27) was used as the gene set, and the permutation number was set to 1,000 as the default. Finally, P-values <0.05 and FDR <25% were considered to be statistically significant (28).

In this study, we used Pearson correlation coefficient to measure the strength of the relationship between the variables. The coexpression of mRNA expression level of PLAGL2 and POFUT1 was presented by linear regression model. Coefficient of determination was calculated and presented. The independent samples t-test was performed for data comparison in GEPIA validation part. All statistical analyses were performed using R program. P-values <0.05 was considered to indicate a statistically significant result.

A workflow of the study is shown in Fig. 1. The dataset GSE39582 contained 585 samples from CRC patients, including 19 normal tissue samples and 566 tumor samples, while GSE14333 had 290 primary CRC tissues. We used the GSE39582 data to build our network and GSE14333 for validation purposes. After data collection, a total of 436 tumor samples with complete clinical information from GSE38582 were obtained. The clinical information of GSE39582 is shown in the clustering dendrogram with the trait heatmap (Fig. 2).

For genes, we transformed the 50,362 probe ids into 22,880 official gene symbols and calculated the median absolute deviation (MAD) of each gene in all samples mentioned above. The three-quarters genes, which equals 17,160, that have the highest MAD were used to construct the final expression network. This step also ensured that the median absolute deviation was not 0, thereby avoiding further errors when constructing the gene coexpression network.

In the meantime, the preprocess of TCGA RNA-seq data was different. We combined the COAD and READ data into one matrix, which has a total of 19,754 genes and 644 samples. Then, we deleted 22 repeat samples and filtered out the genes with zero expression in more than 80% of samples. After voom normalization, we chose the top 12,000 genes with the highest MAD for further analysis.

In choosing the best threshold, we calculated the network topology for soft-thresholding powers from 1 to 20. As shown in Fig. 3A, power value 5, which was the lowest power for the scale-free topology fit index on 0.9, was selected. Afterward, we checked the mean connectivity (Fig. 3B) and double-checked the scale-free topology R² with a linear regression plot (Fig. 3C). Fig. 3D contains a histogram of the frequency of connections. A highly skewed histogram is said to approximate a scale-free network.

The coexpression similarity matrix was then transformed into the adjacency matrix by choosing 5 as a soft threshold, and a topological overlap matrix (TOM) was subsequently computed. Using the dynamic tree cut method, a total of 38 modules were identified. The modules with higher correlation than 0.75 were subsequently merged, resulting in 31 modules at last (Fig. 4). The gray module includes genes that were not assigned to any gene modules.

In the network built by the TCGA dataset, the soft threshold was 7 by the calculation (Fig. 5A). Ultimately, 26 gene modules were recognized (Fig. 5C).

To analyze the relationship between gene modules and sample clinical information, we employed module eigengene (ME) as the average gene expression level of the corresponding modules. It can be considered a representative of the gene expression profiles in a module. The correlations between module eigengene and clinical phenotypes in GSE39582 were calculated and plotted as a labeled heatmap (Fig. 6). The red module and orange module were significantly associated with tumor location.

We calculated gene significance based on the correlation of a gene expression profile with the samples' location traits. Then, the module significance was defined as the average absolute value of the gene significance of all genes in one module. As shown in Fig. 7A, the red and orange modules had considerably stronger correlations with tumor location than did the rest of the modules.

To determine the module's reproducibility, module preservation analysis was performed using an independent dataset GSE14333. As we can see in Fig. 7B, modules below the green dashed line (Z summary <10) are poorly preserved, while the modules above the line are well-preserved in the CRC tissues. The red module, according to the preservation test, is highly preserved in CRC; however, the orange module showed moderate preservation. Thus, we chose the red module for further analysis.

Again, the same method was applied to TCGA data, locating a similar dark-gray module (Fig. 5D).

There were 865 genes in the GSE39582 red module. After plotting the gene significance against module membership, we observed that genes with higher module memberships tended to have higher gene significance in this module (Fig. 7C). We used a relatively high criterion to select hub genes: The absolute value of gene significance >0.4 and module membership >0.8. Six hub genes were successfully identified. The genes with the highest gene significance were found to be POFUT1 and PLAGL2, which are labeled in blue print in Fig. 7C.

Meanwhile, in the TCGA dark-gray module, we used the absolute value of gene significance >0.3 to filter out 8 hub genes (Fig. 5B). After combining two datasets, we determined that there were 12 possible hub genes, 2 of which are in common (Table I).

We concentrated on PLAGL2 and POFUT1 because of their high gene significance and their presence in both datasets. We then evaluated their expression with the online TCGA-based tool GEPIA. PLAGL2 and POFUT1 were found to be significantly differentially expressed between tumor and normal tissue in both the COAD and READ datasets (Fig. 8A and B). We also performed a correlation analysis between PLAGL2 and POFUT1. The plot shows that the Pearson correlation coefficient is tightly correlated to 0.9 in CRC (Fig. 8C).

We utilized quantitative polymerase chain reaction (qPCR) to measure the RNA expression of PLAGL2 and POFUT1 in CRC samples. PLAGL2 had a high positive correlation with POFUT1 according to the qPCR results (Fig. 8D).

For survival analysis, Kaplan-Meier curves were drawn for PLAGL2 and POFUT1 in both proximal and distal CRC (Fig. 9). Although the log-rank P-value of all the analyses was >0.05 (not statistically significant), we still compared the results from different parts of the colon. In proximal CRC samples, there was a clear trend that high PLAGL2 and POFUT1 expression is related to adverse prognosis in CRC patients. However, in distal CRC samples, the expression of POFUT1 was not related to survival, and the high expression of PLAGL2 was even associated with poor survival.

We also performed a Gene Set Enrichment Analysis based on the expression level of PLAGL2 and POFUT1. As shown in Fig. 10, these two genes share a similar enriched KEGG pathway: Glycosylphosphatidylinositol GPI anchor biosynthesis and peroxisome and selenoaminoacid metabolism.