Figure 1 reveals the process of our study. We downloaded the gene expression profiles of GSE62080, which were based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse62080). We chose GSE62080 because this expression profile data include 21 samples from patients with advanced colorectal cancer, and all patients had received chemotherapy based on the FOLFIRI regimen. We used the robust multiple-array average (RMA) algorithm in the affy package of the Bioconductor in R to preprocess the gene expression profile data. We got the gene expression profile dataset containing 23,519 genes after background correction, quantile normalization, and probe aggregation. Then, 10,825 genes which were the top 50% most variant genes by analysis of variance were screened for further analysis.

After performing the WGCNA package in R, a gene co-expression network was constructed by the gene expression data profile from 10,825 genes screened before (Langfelder and Horvath, 2008). The following is the adjacency matrix aij to calculate the connection strength between each pair of nodes: S i , j = c o r X i , X j a i j = S i j β .

In this formula, vectors of the expression value for gene i and j were represented as X_i and X_j. The Pearson correlation coefficient of gene i and gene j was shown as S_ij. The power of β = 7 (scale-free R2 = 0.95) is adopted as the soft-thresholding parameter to ensure a scale-free network. For the hierarchical clustering of the weighting coefficient matrix, the genes with a high correlation were gathered in the same module and identified. We used the adjacency matrix to calculate the topological overlap measure (TOM), which represented the overlap in shared neighbors, and the functional modules of 10,825 genes were determined. T O M i , j = ∑ K = 1 N A i , j × A k , j + A i , j min K i , K j + 1 − A i , j .

In this formula, A is the weighted adjacency matrix given by A i , j = c o r X i , X j β . β = 7 is the soft-thresholding parameter. Then, we set the minimum size of the gene group to be 100 for the gene dendrogram and conducted an average linkage hierarchical clustering based on the TOM-based dissimilarity measure. After that, the DynamicTreeCut algorithm was performed to classify genes with similar expression profiles into the same gene modules.

We used 1 − T O M i , j as the distance to perform hierarchical clustering of genes and then identified modules based on the dynamic pruning method. The most representative genes in each module are called module eigengene (ME), which represents the level of gene expression in the module. ME is the first principal component of each gene module. ME is calculated as follows: M E = p r i n c o m p X i , j q .

In this formula, q is the gene module. Then, we used the Pearson correlation between the expression profile of a certain gene in all samples and the expression profile of a certain ME to measure the module membership (MM) of this gene in the module. MM is calculated as follows: M M i q = c o r x i , M E q . Here, M E q represents the ME of the module q; M M i q represents the MM of the gene i in module q. In addition, the gene significance (GS) was defined as the mediated p-value of each gene in the linear regression between gene expression and clinical characteristics. GS was used to measure the degree of association between genes and clinical characteristics. Moreover, we define the module significance (MS) as the average GS of all genes involved in the module.

The database g: Profiler (https://biit.cs.ut.ee/gprofiler/gost) was used for annotation, visualization, and integrated discovery. In this study, Gene Ontology (GO) and KEGG pathway analysis were technically supported based on the g: Profiler. GO analysis also included the biological process (BP), molecular function (MF), and cellular component (CC).The cut-off criterion of adjusted p < 0.001.

We measured the connectivity of genes by the absolute value of the Pearson correlation. Hub genes were defined as genes with high module connectivity (cor. gene module membership＞0.9). We measured the correlation between hub genes and certain clinical traits based on the absolute value of the Pearson correlation (cor. gene trait significance＞0.2). The clinical characters and RNA sequencing data were screened from The Cancer Genome Atlas project database (TCGA, https://cancergenome.nih.gov/). The edgeR package in R was used to normalize the mRNA sequencing data. Moreover, we used the Human Protein Atlas (HPA, http://www.proteinatlas.org) to identify associations of candidate hub genes and the immunohistochemistry of candidate hub genes. Survival analysis was exhibited on the basis of the Gene Expression Profiling Interactive Analysis (GEPIA) (Tang et al., 2017).

We obtained the expression profile data which contained 21 samples and 21,648 genes. The variance of each gene in each sample was calculated, and the top 50% of genes with a small variance were eliminated. We analyzed and clustered the samples of GSE62080 based on the average linkage method and the Pearson correlation method. Figure 2 portrays the clustering results with sample characteristics. The status was divided into R (change in indicator lesion size＜−50%) and NR (change in indicator lesion size ≥ −50% or progression of disease). The scores were defined as follows: 5: change in indicator lesion size＜−75%; 4: −75% ≤ change in indicator lesion size＜−50%; 3: −50% ≤ change in indicator lesion size＜−25%; 2: −25% ≤ change in indicator lesion size＜0%; 1: change in indicator lesion size = 0%; and 0: progression of disease. We adopted the power of β = 7 as the soft-thresholding parameter to ensure a scale-free network (Figure 3). Then, we set the minimum size of the gene group to be 100 for the gene dendrogram and conducted the average linkage hierarchical clustering based on the TOM-based dissimilarity measure. The ME of each module was calculated, and 19 modules were identified (Figure 4).

According to the eigenvectors of each module, we calculated the correlation between these modules and clinical phenotype and the result is shown in Figure 5. In the figure, the deeper the red color, the stronger the correlation, and the deeper the blue, the weaker the correlation. The figure revealed that the blue module had the strongest correlation with the significant response after FOLFIRI treatment. According to the expression level of each gene in each sample, we, respectively, calculated the correlation between the genes in these modules and status after FOLFIRI treatment. We separately counted the distribution of GS in each module (Figure 6). According to the ME of each module, we calculated the correlation between 19 modules (Figure 7). Through the aforementioned analysis results, the blue module was verified as having the highest correlation with the significant response after FOLFIRI treatment. Therefore, we take the blue module as a key module for further analysis.

All genes in the blue module were divided into the biological process (BP) group, cellular component (CC) group, and molecular function (MF) group. We used Cytoscape (vision 3.8.2) analysis and screened out the top 20 functions of enriched genes in the BP, CC, and MF groups, respectively, and used the EnrichmentMap in Cytoscape for visualization of enrichment results. The genes in the BP group were mainly enriched in extracellular matrix organization, cell adhesion, growth factor response, and anatomical development morphogenesis. The genes in the CC group were mainly enriched in functions associated with collagen extracellular space, banded collagen trimer, and cell substrate junction. The genes in the MF group were mainly enriched in sulfur compound heparin, signaling receptor binding, and constituent conferring elasticity. Moreover, through Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, we found that significant genes were mainly related to ECM–receptor interaction, glycosaminoglycan biosynthesis—chondroitin sulfate/dermatan sulfate, proteoglycans in cancer, protein digestion, and absorption. The results are shown in Figure 8.

We set the cut-off criteria as |MM| > 0.9 and |GS| > 0.2 and selected 50 hub genes with high connectivity in the blue module. Using STRING, we further performed protein–protein interaction (PPI) analysis on these 50 genes and then used MCODE in Cytoscope for further data screening and visualization. Finally, 13 key node genes (COL18A1, TIMP2, DCN, LTBP1, FBN1, TAGLN, EFEMP2, LUM, CDH11, ACTA2, BGN, COL6A2, and AEBP1) were screened out (Figure 9). We further explored the association of the expression levels of the aforementioned genes on overall survival (OS) and disease-free survival (DFS) in colon cancer patients in the GEPIA database. The results showed that the high levels of AEBP1, BGN, EFEMP2, and TAGLN expressions were associated with a shorter OS, and over-expression of AEBP1, BGN, CDH11, COL18A1, LUM, TAGLN, and TIMP2 were negatively associated with DFS (Figure 10). We selected three genes (AEBP1, BGN, and TAGLN) that were negatively related to both OS and DFS for further exploration. A significant correlation between AEBP1 and BGN expressions was found. In addition, AEBP1 and TAGLN expressions also showed a significant correlation. However, the correlation between BGN and TAGLN expressions was not statistically significant (Figure 11). According to the Human Protein Atlas database, we found that the protein expression levels of AEBP1, BGN, and TAGLN in tumor tissues were higher than those in normal tissues (Figure 12).