A transcriptional co-expression network-based approach to identify prognostic biomarkers in gastric carcinoma

Available microarray-based mRNA expression datasets and preprocessing

The training dataset used for co-expression construction was composed of 300 primary GC tumor specimens obtained at the time of total or subtotal gastrectomy from Samsung Medical Center, Seoul, Korea, from 2004–2007. This dataset was also part of the Asian Cancer Research Group (ACRG) study. These data were downloaded from the Gene Expression Omnibus (GEO) database using accession number GSE62254 (Cristescu et al., 2015) and comprised the largest set of samples ever downloaded from the database.

The validation dataset was constructed from 248 primary GC samples from the Singapore patient cohort known as the Gastric Cancer Project ’08. This dataset was used to confirm the relationship of gene modules or biomarkers with survival of GC. Raw data with .CEL profiles from two studies were downloaded from GEO using accession numbers GSE34942 and GSE15459 (Ooi et al., 2009); there were 56 and 192 available samples with detailed information, respectively.

All raw expression data was produced using the Affymetrix Human Genome U133 Plus 2.0 Array™ (HG-U133_Plus_2, Affymetrix, Inc., Santa Clara, CA) and normalized with robust multi-array average (RMA) algorithms (Irizarry et al., 2003) using the affy R package (Gautier et al., 2004). The validation dataset was adjusted for potential batch effects among multiple datasets using the ComBat algorithm (Pavlou et al., 2014). Probe sets with available gene symbols were reserved for subsequent analysis and probe-level expression data were transformed into gene-level expression data by merging the probes according to the official annotation file. The average expression values for the multi-probes were calculated as the corresponding gene expression value for one gene. The primary endpoints for the training dataset were overall survival (OS) and disease-free survival (DFS); overall survival (OS) was regarded as the endpoint event for the validation dataset. Gene expression profiles (Illumina HiSeq RNA Seq), level 3 data, and phenotype data of stomach cancer from The Cancer Genome Atlas Project database (TCGA-STAD) were downloaded through the UCSC Xena portal (https://xena.ucsc.edu/) to validate the identified hub genes. All of the gene expression values were in log₂ (x + 1) transformed normalized count for subsequent analysis. Patients chosen for biomarker identification met the following criteria: (1) histologic diagnosis of primary GC; and (2) available RNA expression profiles and complete clinic-pathological and follow-up data. After sample filtering, 374 patients were enrolled for further analysis. The demographics are listed in Table 1. The hub genes were also validated in the cBio Cancer Genomics portal (https://www.cbioportal.org/) (Gao et al., 2013; Cerami et al., 2012). The Human Protein Atlas (http://www.proteinatlas.org) (HPA) was applied to the protein expression dataset to validate the immunohistochemistry of the identified hub genes. Tissues were defined by the normal tissues of stomach cancer and the pathology was defined by the tumor tissues. The selection process of the prognostic biomarkers is shown in Fig. 1.

Gastric carcinoma molecular subtypes

GC patients were divided into four groups according to their molecular subtypes as described by Cristescu et al. (2015). The molecular subtype signatures had the following tumor biomarkers: MSI (microsatellite instability); MSS/EMT (epithelial-tomesenchymal transition); MSS/TP53⁺ (the tumor protein 53 (TP53)-active); and MSS/TP53⁻ (the tumor protein 53 (TP53)-inactive). Information on the molecular subtypes can be found in the supplement materials of the original publication (Cristescu et al., 2015). In this dataset, 68 samples were classified as MSI, 46 samples were classified as MSS/EMT, 79 samples were classified as MSS/TP53⁺, and 107 samples were classified as MSS/TP53⁻. A subsequent survival analysis according to the molecular subtype was not performed in this dataset because the validation dataset lacked information regarding those subtypes.

Weighted gene co-expression network detection

The training dataset (GSE62254) was the target for WGCNA using the “wgcna” R package (Langfelder & Horvath, 2008). The top 5000 most variable genes were selected according to the median absolute deviation (MAD), which restricted the analysis to those genes with a notable variance in expression; MAD is a robust measure of variability for the construction of co-expression networks.

Pearson’s correlation coefficient (PCC) was used to assess the relationship between each pair of the 5000 genes. These data were used to construct an unsupervised co-expression-based adjacency matrix with a soft threshold power of 4 based on the scale-free topology criterion to raise the matrix to simulate a realistic network structure (Zhang & Horvath, 2005). The intramodular connectivity (k.in) measured how connected, or co-expressed, a given gene was with respect to the genes of a particular module. The connection strengths were assessed by calculating the topology overlap (TO) (Yip & Horvath, 2007), and the modules were defined as sets of genes with a high TO (Zhang & Horvath, 2005). The topological overlap matrix (TOM) was the network distance measure for each gene pair from the adjacent matrix. A TOM-based dissimilarity measure (1-TOM) was utilized to achieve an average hierarchical linkage clustering (Ravasz et al., 2002). Gene modules were defined using a dynamic hybrid branch-cutting algorithm with a cutoff of 0.95 and a minimum module size and cutoff of 30 in the hierarchical clustering dendrogram on the basis of TOM dissimilarity (Langfelder, Zhang & Horvath, 2008). The module eigengene (ME) was calculated by a principal component analysis (PCA) by defining the first principal component of a given module. The module membership, also known as eigengene-based connectivity kME, related each gene expression profile with the ME of a specific module. The MEs of a summary profile were used to assess the underlying correlation of gene modules with the clinico-pathological variables and survival.

Survival analysis and identification of hub genes

Survival analysis was performed using the survival R package with the hazard ratio (HR) and its corresponding 95% confidence interval (CI) determined by the Cox regression module and Kaplan–Meier survival curves (http://cran.r-project.org/web/packages/survival/index.html). OS or DFS were considered to be the survival endpoints. Covariates, including the Lauren’s diffuse type and intestinal type, tumor type, stage and molecular type were corrected via multivariate analysis to estimate the prognostic effects of modules and hub genes. For modules or single gene-based associations, each ME or gene expression value was a continuous variable that was categorized as having a high or low expression according to the median expression value at the cutoff point. For a given ME/gene, the patients were split into two groups known as high expression (≥median expression of the ME/gene) and low expression (<median expression of the ME/gene).

The gene significance (GS) was defined as minus log 10 of the univariate Cox proportional hazard-regression p-values in the single gene-based analysis. Hub genes, namely, highly connected genes, were genes that tended to have high network connectivity (k.in) that determined the connection strength (co-expression) of a specific gene with other genes in a given module (Liu et al., 2017). Genes that satisfied the following criteria were classified as hub genes: (i) GS > 2; (ii) targeted module kME value >0.85.

Functional annotation of the targeted module

Gene enrichment analysis of the categorical biological processes of gene ontology (GO) was conducted on the targeted modules associated with GC patient survival to explore further insights into the genes via the DAVID annotation tool (http://david.abcc.ncifcrf.gov/) (Dennis Jr et al., 2003). The fold enrichment was calculated for all GO terms in the given ontologies to examine the enrichment degree of the specific genes for all genes on the array. Multiple tests were performed with p-values adjusted according to the Bonferroni, Benjamini and false discovery rate (FDR) methods.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for genes in the targeted module was conducted via the ClueGO (Bindea et al., 2009) and CluePedia (Bindea, Galon & Mlecnik, 2013) applications in the Cytoscape v3.3.0 software. Only terms with p-values <0.05 were retained.

Gene set enrichment analysis (GSEA)

To obtain further insight into the potential mechanisms of hub genes significantly associated with survival, GSEA (http://software.broadinstitute.org/gsea/index.jsp) (Subramanian, Tamayo & Mootha, 2005) was conducted based on the median expression level of the significant hub genes to map for the KEGG pathways database. The annotated gene set c2.cp.kegg.v6.1.symbols.gmt was chosen as the reference gene set. Differences with a false discovery rate (FDR) of less than 5% had statistical significance; FDR was calculated using the p.adjust function.

Detection of gene co-expression modules

To investigate the potential patterns of expression among the 5000 most varied genes; WGCNA was conducted on a public microarray-based GC dataset derived from 300 primary GC tumor tissues. Seven informative modules were identified with lengths of 43 to 1,059 genes (Fig. 2A). Each module was assigned a unique color (Table 2), and the gray module, with 2708 non-co-expressed genes, was assigned a gradient color. The topological overlap matrix was plotted according to the expression levels of all genes (Fig. 2B). The MEs and MM (kME) were calculated across all samples and all genes, respectively. The affiliations of all 5000 genes and the complete list of network indices (kME and k.in) for each gene are presented in Data S1.

Correlation of modules with clinico-pathological variables

To measure the association between the seven identified co-expression modules and clinico-pathological variables, we calculated the PCCs between MEs as continuous variables; gender, age, Lauren’s diffuse type and intestinal type, tumor type, tumor stage, molecular subtype, T (Tumor), N (Node), M (Metastasis), node count, positive node, DFS, OS, and status were also calculated (Fig. 2C). The yellow, brown, and turquoise modules yielded significant associations with Lauren’s type (yellow = −0.25, brown = +0.3, turquoise = +0.34), tumor type (yellow = −0.36, brown = +0.28, turquoise = +0.42), tumor stage (yellow = −0.23, brown = +0.22, turquoise = +0.26), and molecular subtype (yellow = −0.22, brown = +0.52, turquoise = +0.45).

Gene modules significantly correlated with survival

The relationship between OS/DFS and modules was assessed as a whole using Cox regression analysis with HRs and corresponding p-values (Table 2). The black, brown, green, turquoise, and yellow modules showed significant correlations with OS in the training dataset (module brown: HR = 1.586, p = 0.005, 95% CI [1.149–2.189]). However, only the brown module was confirmed in the validation dataset (HR = 1.664, p = 0.006, 95% CI [1.155–2.398]). We also found that ME_black, ME_brown, ME_green, and ME_turquoise were significantly correlated with DFS in the training dataset. The brown module was the focus in subsequent analyses as its high expression demonstrated poor prognosis in the training (Figs. 3A and 3B) and validation datasets (Fig. 3C).

Biological insights into the brown module

To illuminate the potential biological insights into the brown module, we performed GO biological process enrichment analysis with DAVID and KEGG pathway analyses using Cytoscape v3.3.0. Seven GO terms were identified that were significantly enriched with FDR < 0.05 (Fig. 3D) and five pathways with FDR < 0.05 (Fig. 3E). The most significant terms from the GO analysis were “extracellular matrix organization” (raw p-value = 6.29 × 10⁻¹³, Bonferroni-adjusted p-value = 9.68 × 10⁻¹⁰, FDR = 8.11 × 10⁻¹¹) and “ECM-receptor interaction” (raw p-value = 2.40 × 10⁻⁵, Bonferroni-adjusted p-value = 3.40 × 10⁻⁴, FDR = 1.20 × 10⁻⁴) from the KEGG analysis. A full list of biological GO terms and the KEGG analysis of all co-expression genes in the brown module is shown in Data S2 and S3, respectively.

Identification and validation of survival-related hub genes

Hub genes are highly likely to serve as key factors in a given module. Ten hub genes were identified from the 342 co-expression genes in the brown module (COL8A1, FRMD6, DDR2, LOC100505881, TIMP2, CNRIP1, CLEC11A, MRC2, BGN, and GPR124). These were associated with OS/DFS in the training dataset and were confirmed in the validation dataset with the exception of the relationship between CNRIP1, and OS (Table 3). Survival analysis of the hub genes COL8A1, FRMD6, TIMP2, CNRIP1, and GPR124 in the training dataset and validation dataset is shown in Fig. 4. The high expression of the five hub genes presented with poor overall survival and was validated in the TCGA dataset (Fig. 5). We also investigated the hub genes in other modules based on the above screening criterion. The results are presented in Data S1 with highlights. There were no hub genes in the blue, red, and gray modules. Immunohistochemistry (IHC) staining obtained from the Human Protein Atlas database also demonstrated the expression status of hub genes (Fig. 6). There were no related IHC samples of CNRIP1 in the database. The Oncomine database was used in our analysis and the mRNA levels of COL8A1, TIMP2, and GPR124 were higher in tumor tissues compared with normal tissues (Fig. 7A). The OncoPrint module in cBioPortal, an online tool, was used to analyze genetic alterations, including missense mutations, truncating mutations, amplifications and deep deletions (Fig. 7B).

Identification of gene modules and hub genes significantly associated with GC subtype-specific survival

The relationship between the gene modules or hub genes and GC molecular subtypes was explored and HRs and accompanying p-values were calculated to denote their significance. Analysis revealed that ME_green (HR = 1.765, 95% CI [1.043–2.99], p-value = 0.034) and ME_turquoise (HR = 2.485, 95% CI [1.205–5.124], p-value = 0.014) involving 113 and 1059 genes, respectively, were associated with poor OS outcomes within the MSS/TP53⁻ and MSS/EMT molecular subtypes (Table 4). The increased expression of genes in module turquoise indicated poor DFS prognosis in the MSS/EMT molecular subtype (HR = 2.467, 95% CI [1.198–5.081], p-value = 0.014) (Table S1).

Furthermore, we also investigated whether significant correlations could be detected among the 10 hub genes in the brown module and the molecular subtypes in the training dataset. Association analysis revealed that the expression levels of genes FRMD6 (HR = 2.822, 95% CI [1.156–6.888], p-value = 0.023), TIMP2 (HR = 2.942, 95% CI [1.199–7.217], p-value = 0.018), CNRIP1 (HR = 2.626, 95% CI [1.077–6.401], p-value = 0.034) and GPR124 (HR = 3.696, 95% CI [1.451–9.413], p-value = 0.006) were significantly associated with OS in the MSI molecular subtype. The increased expression of the COL8A1 gene showed poor prognosis (HR = 2.216, 95% CI [1.104–4.447], p-value = 0.025) in the MSS/EMT molecular subtype (Table 5). We also found that the FRMD6 (HR = 2.715, 95% CI [1.116–6.605], p-value = 0.028), TIMP2 (HR = 2.964, 95% CI [1.217–7.217], p-value = 0.017), CNRIP1 (HR = 2.942, 95% CI [1.207–7.168], p-value = 0.018), MRC2 (HR = 2.379, 95% CI [1.007–5.62], p-value = 0.048) and GPR124 (HR = 3.814, 95% CI [1.5–9.694], p-value = 0.005) genes were significantly associated with DFS in the MSI molecular subtype. In addition, COL8A1 gene expression was associated with DFS prognosis in the MSS/EMT molecular subtype (HR = 2.295, 95% CI [1.144–4.602], p-value = 0.019) (Table S2). There were significant differences noted in the expression of five hub genes when the tumor stages and molecular subtypes were compared (Figs. 8A and 8B).

GSEA analysis of key hub genes significantly correlated with subtype-specific survival of GC patients

To characterize the potential function of the real hub genes, 300 GC samples were divided into two groups (high vs. low) according to the median expression values of the above hub genes. GSEA was performed based on the expressions of COL8A1, FRMD6, TIMP2, CNRIP1 and GPR124. GC samples in the COL8A1 and CNRIP1 high expression group were significantly enriched for focal adhesion (Tables S3–S6). GC samples in the FRMD6 and TIMP2 high expression group were significantly enriched for hypertrophic cardiomyopathy (HCM) (Fig. 9; Tables S4–S5). GC samples in the GPR124 high expression group were significantly enriched for dilated cardiomyopathy (Fig. 9; Table S7).