Identification of subtype-specific genes signature by WGCNA for prognostic prediction in diffuse type gastric cancer

Detection of gene co-expression modules correlated with diffuse type GC cohort

The data was processed and analyzed by following the workflow in Figure 1. Top 50% most variable genes (9752 genes) were used for WGCNA. An obvious outlier was removed (Supplementary Figure 1A) and a soft threshold = 4 was selected to construct a scale-free network (Supplementary Figure 1B, 1C). A total of 30 gene modules were identified after setting the minimum cluster size as 30 (Figure 2A). The grey module contained genes not attributed to any modules. The identified genes associated with a clinical trait were of great value in the exploration of the molecular characteristics of that trait. In the present study, the clinical parameters of diffuse type GC patients, including age at diagnosis, gender, race, tumor differentiation, pathological stage, pathological T category, pathological N category, M category, lymphatic metastasis, cancer status, helicobacter pylori infection, status and overall survival (OS) time was extracted for analysis. we found that blue module (R=0.26, p=3e-07) and green module (R=0.27, p=7e-08) were significantly associated with diffuse type GC (Figure 2B–2D) and the blue (R=0.22, p=3e-05) and green (R=0.26, p=4e-07) module were also significantly correlated with tumor differentiation.

Figure 1. The flowchart of identifying procedure for the multi-gene signatures in diffuse type GC.

Figure 2. Identification of modules associated with the diffuse type GC. (A) Dendrogram of 9752 genes clustered based on a dissimilarity measure (1-TOM). (B) Correlation of module eigengenes with all traits. Each unit contains the corresponding correlation coefficient and P-value. The table is color-coded by correlation according to the color legend. (C, D) Scatter diagrams between blue and green modules and diffuse type GC.

Hub genes in blue and green module

Before conducting the univariate cox analysis, we first adopted the merge function in R to integrate the expression profiles of the 1794 module genes with corresponding 71 diffuse type GC patients’ survival time and status information. 63 overall survival (OS)-related hub genes (p<0.1) were picked out by univariate cox regression analysis. According to the characteristics of variable selection and regularization, while fitting the generalized linear model, LASSO regression was performed to select hub genes for predicting the prognosis of high-performance patients (Figure 3A, 3B). This approach is popular in machine learning and is implemented through the “glmnet” package. 16 hub genes ("RALA", "DDX3Y", "ERP29", "SRSF5", "SLC9A3R1", "FBXO9", "GMDS", "CCNI", "LEF1", "RFX5", "CAST", "ELMO1", "FRZB", "TMEM92", "SELP" and "NMB") were identified.

Figure 3. Identification of hub genes using LASSO regression. (A) The trajectory of each prognosis-related candidate gene’s coefficient in diffuse type GC was observed in the LASSO coefficient profiles with the changing of the lambda in LASSO algorithm. (B) After the 10-fold cross-validation, a confidence interval was got for partial likelihood deviance as the lambda changed.

Construction of risk score model based on multivariate Cox regression

First, followed by multivariate cox regression, the optimal 10 prognostic signatures in diffuse type GC samples (Table 1), including: “RALA”, “DDX3Y”, “SRSF5”, “SLC9A3R1”, “GMDS”, “LEF1”, “RFX5”, “CAST”, “FRZB” and “SELP”, were picked out. Then, the risk score (RS) model for OS was identified. All 71 diffuse type GC samples were endowed with a RS by RS model calculation and divided into high- and low- RS groups using the median value as the cut-off point. The concordance index (C-index) of this model was 0.97, indicating that this model had quite high reliability. Figure 4A showed that patients in the low-risk group had longer OS (p < 2.8e-10) than those of the high-risk group. The relation of high- and low expression of 10 genes to OS were viewed in Supplementary Figures 2–11. To determine the predictive accuracy of this prognostic model, we performed a receiver operating characteristic (ROC) curve analysis, which demonstrated that the area under the curve (AUC) was 0.98 for 1-year survival, 1 for 3-year survival and 1 for 5-year (Figure 4B). And the calibration curve supported the predicting model of 10 genes (Figure 4C–4E). The survival, 10 genes expression and risk score for TCGA samples could be viewed in Supplementary Table 1.

Figure 4. The prognostic performance of the 10 genes model in the TCGA-STAD. (A) Survival analysis of the high-risk group and the low-risk group using Kaplan–Meier curves. (B) The prognostic efficiency of the 10 genes model for survival time. ROC curves of the 10 genes signature for predicting 1-, 3 -and 5- year survival were analyzed. (C–E) The comparison between predicted and actual outcome for 1-, 3-, and 5-year survival probabilities was showed in the calibration plots.

Validation of 10-genes RS model in external independent cohort

75 diffuse type GC patients from GSE15459 were remaining after removing the non-diffuse type and survival time as zero. Consistent with the results in the TCGA cohort, the low RS group had good performance in OS (P=0.00633473) than in the high RS group (Figure 5A). Moreover, the AUCs for 1-year, 3-year and 5-year survival in the validation cohort were 0.717, 0.727 and 0.7 (Figure 5B), respectively. There were 135 diffuse type GC patients in GSE62254 by the above same processing. The low RS group had good performance in OS (P=2E-8) than the high RS group as same as above (Figure 5C). the AUCs for 1-year, 3-year and 5-year survival in the validation cohort were 0.661, 0.752 and 0.758 (Figure 5D). Thus, the RS model had again been proved to be reliable. The survival, 10 genes expression and risk score for GSE15459 and GSE62254 samples could be viewed in Supplementary Tables 2 and 3.

Figure 5. The prognostic performance of the 10 genes model in the GSE15459 and GSE62254. (A) Survival analysis of the high-risk group and the low-risk group using Kaplan-Meier curves in the GSE15459. (B) The prognostic efficiency of the 10 genes model for survival time. ROC curves of the 10 genes signature for predicting 1-, 3 -and 5- year survival were analyzed in the GSE15459. (C) Survival analysis of the high-risk group and the low-risk group using Kaplan-Meier curves in the GSE62254. (D) The prognostic efficiency of the 10 genes model for survival time. ROC curves of the 10 genes signature for predicting 1-, 3 -and 5- year survival were analyzed in the GSE62254.

GSEA

In order to explore the difference of functions and pathways of high- and low- RS groups, the gene set enrichment analysis (GSEA) was performed using the risk score as the reference phenotype. The GSEA analysis revealed that the complement and coagulation cascades, neuroactive ligand-receptor interaction, hypertrophic cardiomyopathy (HCM), steroid hormone biosynthesis and dilated cardiomyopathy were upregulated in the high-risk group (Figure 6A). Then, the spliceosome, one carbon pool by folate, nucleotide excision repair, cell cycle, RNA degradation, mismatch repair, DNA replication, ubiquitin mediated proteolysis, homologous recombination, p53 signaling pathway, basal transcription factors and base excision repair were upregulated in the low-risk group (Figure 6B). To investigate the enrichment pathways of 10 prognostic genes, we divided the samples into two groups, a group samples with the expression of upper quantile for one of 10 prognostic genes and another group samples with the expression of lower quantile for one of 10 prognostic genes. All 10 genes related pathway enrichment could be viewed in Supplementary Table 4.

Figure 6. (A, B) GSEA results revealed the significantly enriched biological processes between two RS levels.

Data collection and processing

Public gene-expression data and full clinical annotation were searched in the Gene-Expression Omnibus (GEO) and TCGA database. The procedure used for data set selection in the GEO database was as follows. The following search parameters were used: (cancer) OR tumor) OR carcinoma) OR adenocarcinoma) AND ((gastric) OR Stomach)) AND "Homo sapiens"[porgn: txid9606]. In the initial search, 755 items were recognized. The eligible criteria included that: 1) owning Lauren classification; 2) owning survival information and 3) at least 50 diffuse type gastric cancer. We removed the datasets that don’t meet the criteria by checking them one by one carefully. In total, we gathered three patient cohorts with gastric cancer for this study: TCGA-STAD, GSE62254 and GSE15459.

The TCGA-STAD RNA-seq data and clinical data of the 375 GC samples were downloaded by using the “TCGAbiolinks” package in R [20], which was used as a training set for prognostic prediction of the multi-gene signature. The RNA-seq data for 19505 genes measured as fragments per kilobase of transcript per million mapped reads (FPKM), which were converted to transcripts per million (TPM) after removing duplicated genes and zero expression genes [21]. We obtained the raw expression and clinical data of GC samples from GSE622254 and GSE15459 via GEO database which was used as a validation set. The raw expression data was processed by log10-transformed. The demographics are listed in Table 2.

Weighted Gene co-expression network construction

To find modules of highly correlated with diffuse type GC, WGCNA was performed using the WGCNA R package [22] and carried out on top 50% most variable genes (9752 genes). An unsupervised co-expression relationship was initially built on the basis of the adjacency matrix of connection strengths by using Pearson’s correlation coefficients for gene pairs. This matrix was increased to β = 4 based on the scale-free topology criterion. The adjacency matrix of gene expression data for GC patients was then clustered using topological overlap matrix analysis. Finally, the dynamic tree cut algorithm was applied to the dendrogram for module identification with the mini-size of module gene numbers set as 30 and a cut height of 0.95. The module eigengenes (MEs) as the first principal component was performed with the expression data for each co-expressed module in all GC samples. The module that had the strongest association with diffuse type GC was selected for further analysis. The WGCNA algorithm was described in detail by Zhang Bin et al. [23].

The construction of multi-gene signature risk score model

To identify the prognostic genes in the module correlation with diffuse type GC, the univariate Cox regression analysis was applied using “survival” package. The genes with P<0.1 were defined to be related to the over survival. Given the already detected prognostic genes, we further investigated the significant signature associated with survival across the diffuse type GC samples by the LASSO regression model using “glmnet” package. Then, multivariate Cox regression analysis was performed with the “survival” package to screen out independent prognostic factors from these robust markers, which were conducted the risk score model: prognostic score =∑(C×EXP_mRNA), where EXP is the TPM value of the gene and C is the regression coefficient for the corresponding gene in multivariate Cox hazard model analysis. ROC plot with AUC and calibration curve were derived to assess predictive significance of the model using “rms”, “survival” and “timeROC” package in R. Additionally, survival time difference between high- and low-risk levels were estimated via Kaplan-Meier analysis. The risk score model was validated with an independent data set GSE15459 and GSE62254.

Functional pathway analysis

To investigate the biological function of different risk score groups, we further conducted GSEA using the risk score as the phenotype [24]. With the GSEA 4.0.3 software via the Java platform, we derived the “c2.cp.kegg.v7.1.symbols.gmt gene sets” as the reference set. We divided the samples into two groups, a group samples with the expression of upper quantile for one of 10 prognostic genes and another group samples with the expression of lower quantile for one of 10 prognostic genes. The enriched signaling pathways with FDR < 0.25 or nominal p < 0.05 were defined as statistically significant.