Weighted correlation network analysis identifies FN1,COL1A1 and SERPINE1 associated with the progression and prognosis of gastric cancer

Abstract

BACKGROUND:

Gastric cancer (GC) is one of the most deadliest tumours worldwide, and its prognosis remains poor.

OBJECTIVE:

This study aims to identify and validate hub genes associated with the progression and prognosis of GC by constructing a weighted correlation network.

METHODS:

The gene co-expression network was constructed by the WGCNA package based on GC samples and clinical data from the TCGA database. The module of interest that was highly related to clinical traits, including stage, grade and overall survival (OS), was identified. GO and KEGG pathway enrichment analyses were performed using the clusterprofiler package in R. Cytoscape software was used to identify the 10 hub genes. Differential expression and survival analyses were performed on GEPIA web resources and verified by four GEO datasets and our clinical gastric specimens. The receiver operating characteristic (ROC) curves of hub genes were plotted using the pROC package in R. The potential pathogenic mechanisms of hub genes were analysed using gene set enrichment analysis (GSEA) software.

RESULTS:

A total of ten modules were detected, and the magenta module was identified as highly related to OS, stage and grade. Enrichment analysis of magenta module indicated that ECM-receptor interaction, focal adhesion, PI3K-Akt pathway, proteoglycans in cancer were significantly enriched. The PPI network identified ten hub genes, namely COL1A1, COL1A2, FN1, POSTN, THBS2, COL11A1, SPP1, MMP13, COMP, and SERPINE1. Three hub genes (FN1, COL1A1 and SERPINE1) were finally identified to be associated with carcinogenicity and poor prognosis of GC, and all were independent risk factors for GC. The area under the curve (AUC) values of FN1, COL1A1 and SERPINE1 for the prediction of GC were 0.702, 0.917 and 0.812, respectively. GSEA showed that three hub genes share 15 common upregulated biological pathways, including hypoxia, epithelial mesenchymal transition, angiogenesis, and apoptosis.

CONCLUSION:

We identified FN1, COL1A1 and SERPINE1 as being associated with the progression and poor prognosis of GC.

Keywords

WGCNA gastric cancer progression prognosis hub genes prediction diagnosis

Abbreviation

AUCs	The areas under the curves
BP	biological process
CC	cellular component
COL11A1	Collagen Type XI Alpha 1 Chain
COL1A1	Collagen Type I Alpha 1 Chain
COL1A2	Collagen Type I Alpha 2 Chain
COMP	Cartilage Oligomeric Matrix
	Protein
DEGs	differentially expressed genes
ECM	extracellular matrix
FDR	false discovery rate
FN1	Fibronectin 1
GC	gastric cancer
GEO	Gene Expression Omnibus
GEPIA	Gene Expression Profiling
	Interactive Analysis
GO	Gene Ontology
GSEA	Gene set enrichment analysis
KEGG	Kyoto Encyclopedia of Genes and
	Genomes
MAD	median absolute deviation
MCC	maximal clique centrality
MF	molecular function
MMP13	Matrix Metallopeptidase 13
POSTN	Periostin
PPI	Protein protein interaction
ROC	receiver operating characteristic
SERPINE1	plasminogen activator inhibitor 1
	(PAI-1)
SPP1	Secreted Phosphoprotein 1
TCGA	The Cancer Genome Atlas
THBS2	Thrombospondin 2
TOM	topological overlap matrix
WGCNA	Weighted gene co-expression
	network analysis

1. Introduction

According to the 2018 Global Cancer data, GC is the fifth most common and the third most deadliest tumour worldwide [1]. In 2020, it is estimated approximately 27600 new cases of GC and 11010 GC deaths in the United States [2]. The high incidence and mortality of GC cause a very large burden of disease globally and individually each year. The prognosis of GC depends on the pathological stage, tumour site, histological type, biological behaviour and therapeutic measures [3, 4]. Tumours with late stage or high-level grade usually have a poor prognosis. A large retrospective study from Japan revealed that the 5-year overall survival rate for stage IV GC patients was far lower than that for stage IA GC patients (16.4% vs 91.5%) [5]. In fact, most GC patients were diagnosed at advanced stages and have an unfavourable prognosis due to a lack of evaluation and treatment strategies. Therefore, it is of great significance to find new therapeutic and prognostic molecular biomarkers of GC.

Currently, emerging microarray and high-throughput sequencing technologies have provided a large amount of available sequencing data. Bioinformatics analysis can help researchers discover new drugs and new molecular markers related to cancers [6, 7, 8, 9]. The Cancer Genome Atlas (TCGA) analysis project [10] and Gene Expression Omnibus (GEO) [11] are the two most widely used databases for researchers. Weighted gene co-expression network analysis (WGCNA) is an R package widely used to identify modules related to different clinical traits [12]. By clustering co-expression genes based on expression data at the gene level, WGCNA can identify significant modules responsible for clinical traits. Some studies that have identified some gene markers for GC by WGCNA commonly use differentially expressed genes to construct networks [13, 14, 15, 16], which may destroy genetic diversity and not reflect the true cluster results. Instead, top variant genes are more reasonable for WGCNA. The pathological stage and histological grade are two important factors affecting prognosis. Identifying the module of interest related to stage, grade and overall survival simultaneously might help to identify true genes that could act as diagnostic and prognostic biomarkers of GC. However, related research has not been reported.

Figure 1.

Flow chart of the study design. Firstly, the RNA-seq and clinical data were downloaded from the TCGA database and top 5000 variant genes were selected to construct WGCNA network. Secondly, by WGCNA network, we identified the interest module that was related to survival status, stage and grade, in addition, the GO and KEGG enrichment analyses of interest module were also explored. Thirdly, we selected 10 hub genes from the interest module by using cytoscape software and explored their expression and prognostic value in gastric cancer (GC), in which we identified three hub genes that were differentially expressed in GC and related to poor prognosis of GC. At last, we performed verification analysis by using several sequencing data from GEO database and our experiments including immunohistochemical staining (IHC), western blot (WB) and quantitative real-time fluorescence PCR (qRT-PCR) with gastric tissue specimens. Meanwhile, univariate and multivariate cox analysis, GSEA and ROC curve were also performed.

In this study, we first performed WGCNA to find significant modules correlated with clinical traits including age, sex, stage, grade, and overall survival (OS) of GC based on GC samples from the TCGA database. Then Gene Ontology (GO) and KEGG pathway enrichment analyses were performed to identify the biological significance of the selected modules. Third, the 10 hub genes of interest module were selected by a PPI network using Cytoscape software and were analysed at mRNA expression levels and prognostic values on the GEPIA online resource. Furthermore, we selected statistically significant genes at the expression and prognosis levels as true hub genes. In addition, we performed verification experiments using GEO datasets and our clinical gastric specimens. The ROC curves were also plotted in the prediction of GC. Finally, GSEA was conducted to explore the potential pathogenic mechanism of true hub genes.

2. Materials and methods

2.1 Data collection from the TCGA and GEO database

The workflow for the current study is presented in Fig. 1. The HTSeq count matrix of 407 gastric samples (375 tumours and 32 normal tissues) was downloaded from the TCGA database using Genomic Data Commons Client (https://portal.gdc.cancer.gov/). The clinical data of the TCGA-STAD dataset were extracted from supplementary materials of integrated the TCGA Pan-Cancer Clinical Data published in the journal of Cell [17]. After removal of 32 normal samples and 40 samples with incomplete clinical data, a total of 335 GC samples with 19712 protein coding genes were included in the study. The microarray data of GSE54129 and GSE118916 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GEO expression matrix was annotated with gene symbols and log2 transformed in R (version 3.6.1) if necessary.

2.2 Data process for WGCNA

For WGCNA, we used the voom function [18] in the limma [19] package to normalize the TCGA count data. A large number of genes are not differentially expressed between samples. These genes must be excluded because they may produce unreal correlation. In this study, we selected the top 5000 most variant genes for WGCNA. The median absolute deviation (MAD) was used as a robust measure of variability [20, 21].

2.3 Gene co-expression network construction

A gene co-expression network was constructed using the R WGCNA package [12]. The samples were clustered with hierarchical clustering analysis by the hclust function [22] to check if outlier samples existed. The function pickSoftThreshold [23] was used to select appropriate soft thresholding powers $\beta$ value based on a set number from 1 to 20 when the scale-free topology fitting indices R2 reached 0.9 to ensure a scale-free network. One-step network construction method was used. The minimum module size of 30 and medium sensitivity (deepSplit $=$ 2) were used in cluster splitting to generate different modules.

2.4 Identification of module of interest

The moduleEigengenes function was used to calculate the module eigengene (ME). The correlation between MEs and clinical traits, including age, sex, stage, grade, survival status (OS), and overall survival time (OS.time), was evaluated by Pearson’s correlation tests. $P$ values less than 0.05 were considered to be significantly correlated. To explore the potential biological significance of selected modules, GO and KEGG pathway enrichment analyses were performed with the clusterprofiler package in R [24].

2.5 Identification of hub genes with a PPI network

Hub genes usually refer to those with a high degree of intra-module connectivity. To identify hub nodes, a PPI network based on the selected module was constructed in the STRING database (https://www.string-db.org/) and visualized in Cytoscape software [25]. The maximal clique centrality (MCC) algorithm of the cytoHubba plugin, which was reported to be effective predicting hub nodes in the PPI network [26], was used to calculate 10 hub genes. The correction analysis of 10 hub genes was performed using the ggcorrplot function in R.

2.6 Differential expression and survival analysis

Gene Expression Profiling Interactive Analysis (GEPIA) is a bioinformatics platform based on expression and survival information from the TCGA database, as well as expression data from the GTEx database [27]. To analyse the expression levels of 10 hub genes between GC and normal gastric tissue, data from both the TCGA and GTEx databases were used. One-way ANOVA was used for the hypothesis test. Genes with $|$ log2FC $|$ cutoff $=$ 1 and $P$ -value cut-off $=$ 0.01 were considered differentially expressed genes (DEGs). To explore the prognostic values of ten hub genes, overall survival (OS) analysis was conducted on the GEPIA online resource. The median expression levels of each hub gene were used as the cut-off value to define the high and low groups. Log-rank test was used for hypothesis test. A $P$ value less than or equal to 0.05 was considered statistically significant.

2.7 Verification of differentially expressed and prognostic genes by GEO databases

To verify the expression levels of hub genes, we selected two different types of GC datasets. GSE54129 contains 21 normal gastric tissues and 111 GC samples; GSE118916 contains 15 paired GC and adjacent normal tissues. The mRNA expression levels of the GC and normal groups were compared with $t$ -tests (GSE54129) or paired $t$ -tests (GSE118916) and plotted with GraphPad Prism 7.0. Furthermore, to test the prognostic value of selected hub genes, we selected another two GEO databases (GSE15459, and GSE62254) because they have a large sample size and complete clinical data. The survival analyses were conducted on the Kaplan-Meier online platform [28]. $P<$ 0.05 was considered statistically significant.

2.8 Gastric tissue specimens and immunohistochemistry

We collected 20 pairs of paraffin-embedded human GC and normal gastric tissue specimens from surgical samples at the First Affiliated Hospital of Nanchang University. All patients signed informed consent forms and the study was approved by the Ethics Committee of the First Affiliated Hospital of Nanchang University. Immunohistochemical staining was conducted at a concentration of 1:1200 for FN1 (66042-1-Ig), 1:2000 for COL1A1 (67288-1-Ig) and 1:1800 for SERPINE1 (66261-1-Ig). All antibodies were purchased from Proteintech ${}^{\text{TM}}$ Company (Wuhan, Hubei, China).

Table 1
The primers used in the real-time quantitative PCR analysis

Genes	Primers	Sequences	Product size (bp)
FN1	Forward	GAGAATAAGCTGTACCATCGCAA	200
FN1	Reverse	CGACCACATAGGAAGTCCCAG
COL1A1	Forward	AAGGTGTTGTGCGATGACG	116
COL1A1	Reverse	TGGTCGGTGGGTGACTCTG
SERPINE1	Forward	GACCGCAACGTGGTTTTCTC	131
SERPINE1	Reverse	GCCATGCCCTTGTCATCAAT
GAPDH	Forward	ACCCAGAAGACTGTGGATGG	124
GAPDH	Reverse	TCAGCTCAGGGATGACCTTG

2.9 Western blotting

Twenty pairs of GC and normal gastric tissues preserved at $-$ 80 ${{}^{\circ}}$ C were dissolved on ice, and a certain amount of RIPA tissue cell lysate (Cat# R0020, Solarbio, China) with Halt ${}^{\text{TM}}$ protease and phosphatase inhibitor cocktail (Cat# 78444, Thermo Fisher, USA) was added for electric grinding. A Pierce ${}^{\text{TM}}$ BCA Protein Assay Kit (Cat# 23227, Thermo Fisher, USA) was used to detecting protein concentration. The tissue lysates were boiled and denatured at 100 ${}^{\circ}$ C for 10 minutes with protein loading (Cat# DL101,TransGen, China). For western blotting (WB), a total of 30 $\mu$ g protein was separated by electrophoresis on a 10% SDS-PAGE gel and transferred to nitrocellulose membranes (NC) under wet conditions. The membranes were then blocked in 5% milk-TBST(g/ml) at room temperature for 1.5 hours. After that, the NC membranes were cut into strips according to corresponding molecular weight and incubated with primary antibodies against FN1 (1:2000, 66042-1-Ig), COL1A1 (1:5000, 67288-1-Ig), SERPINE1 (1:1000, 66261-1-Ig) and GAPDH (1:5000, 60004-1-Ig) in 5% BSA-TBST at 4 ${{}^{\circ}}$ C overnight. The next day, the strips were washed three times with TBST and incubated with anti-mouse secondary antibodies (1:5000, Cat# HS201-01, TransGen, China) in 2% milk-TBST at room temperature for 1 hour. After that, the strips were washed three times with PBS again and imaged by the iBright FL1000 smart imaging system (Thermo Fisher, USA). The proteins of interest were quantified by comparison with the greyscale value of GAPDH using ImageJ software.

2.10 RNA extraction and real-time quantitative PCR analysis

Total RNA was extracted using TRIzol reagents (Takara, Dalian, China) and reverse transcribed into cDNA using the PrimeScript RT reagent kit (Takara, Dalian, China) according to the instructions as the same with our previous study [29]. SYBR Green Master Mix (Thermo Fisher Scientific) was used to detect the mRNA expression levels of FN1, COL1A1 and SERPINE1 in the QuantStudio ${}^{\text{TM}}$ 5 Real-time PCR system (Thermo Fisher Scientific). The endogenous GAPDH mRNA level was used to normalize the relative mRNA levels of three genes of interest. The primers used in the study are shown in Table 1.

2.11 Diagnostic efficiency evaluation for GC

To explore whether the mRNA levels of the three hub genes exhibit diagnostic efficiency for distinguishing GC from normal gastric tissues, we performed receiver operating characteristic (ROC) curve analysis. Briefly, the expression profile of 407 gastric samples (375 tumour tissues and 32 normal tissues) was downloaded from the TCGA-STAD dataset and log2 transformed in R (version 3.6.1). The pROC package [30] was used to plot ROC curves and calculate the areas under the curves (AUCs) values in R (version 3.6.1).

2.12 GSEAs of three hub genes based on the TCGA RNA-seq data

Gene set enrichment analysis (GSEA) is a computational approach for identifying potential biological pathways based on RNA expression profiles [31]. To explore the possible carcinogenic process of selected hub genes, we performed GSEA using GSE software (version 4.3.0) based on the TCGA data. A total of 375 GC samples were divided into high and low expression groups according to the median expression value of each hub gene. The hallmark gene_sets (h.all.v7.0 symbols.gmt) were used for analyses. The number of permutations was 1000, and the permutation type was phenotype (gene_high vs gene low). NES (normalized enrichment score) $>$ 2 and FDR (false discovery rate) $q$ -value $<$ 0.01 were considered significantly enriched.

Figure 2.

Construction of WGCNA network. (A), Clustering dendrogram of 335 GC samples and clinical traits. Colors: white means low, red means high. (B), Visualizing the gene network with a heatmap plot. The heatmap depicts the topological overlap matrix (TOM) among all genes in the analysis. (C), Clustering dendrogram of genes together with assigned module colors. (D), Analysis of network topology for various soft-thresholding powers.

Figure 3.

Identification of the module of interest. (A), Heatmap of module-trait relationships. (B), The clustering eigengenes of network among different modules and grade. (C), The eigengene adiacency heatmap of network among different modules and grade. (D) The clustering eigengenes of network among different modules and stage. (E) The clustering eigengenes of network among different modules and OS. (F–H), Scatterplot of gene significance and Module Membership in the magenta module for OS (F), grade (G) and stage (H).

3. Results

3.1 Gene co-expression network and module analysis

The top 5000 most variant genes based on the MAD value in 335 GC samples from the TCGA database were first used for hierarchical clustering analysis. There were no outlier samples to be excluded. The relationship between clinical features and the sample dendrogram was visualized (Fig. 2A). The WGCNA network was then visualized with a topological overlap matrix (TOM) heatmap of all genes (Fig. 2B). According to the clustering dendrogram of genes correlated with different modules, ten co-expressed modules were identified (Fig. 2C). In this study, power 4 (scale free R2 $=$ 0.9) was selected as the appropriate soft-thresholding parameter to ensure a scale-free network for the whole analysis (Fig. 2D).

3.2 Identification of modules associated with clinical characteristics

To identify modules that were significantly associated with clinical traits, the module-trait relationship was visualized with a heatmap (Fig. 3A). The results showed that the yellow module was highly correlated with grade and stage. Nevertheless, it is noteworthy that the magenta module was positively associated with stage, grade, and OS and negatively associated with OS.time simultaneously. To quantify module similarity with grade, stage or OS by eigengene correlation, an eigengene network was conducted. The eigengene dendrogram of grade (Fig. 3B) and the eigengene adjacency heatmap of grade (Fig. 3C) showed the magenta module was more relevant to grade than the yellow module. Similarly, we found that magenta module was also highly clustered with stage (Fig. 3D) and OS (Fig. 3E) by eigengene dendrogram. Furthermore, scatter plots of gene significance (GS) and module membership (MM) in magenta module analysis showed that magenta module genes were highly related to OS (Fig. 3F), stage (Fig. 3G) and grade (Fig. 3H). Therefore, the magenta module was considered as the most relevant to progression and prognosis of gastric cancer and selected for the subsequent analysis.

3.3 Enrichment analyses of the selected module

To explore the biological process of the magenta module, 83 genes in the magenta module were included to perform GO and KEGG pathway enrichment analyses. GO enrichment analysis showed that extracellular matrix, collagen fibril organization in the biological process (BP), collagen-containing extracellular matrix, collagen, fibrillar collagen trimer in cellular component (CC) and extracellular matrix, collagen binding and metallopeptidase activity in molecular function (MF) were significantly enriched (Fig. S1A). KEGG pathway analysis showed that ECM-receptor interaction, focal adhesion, PI3K-Akt pathway, and proteoglycans in cancer were mainly enriched (Fig. S1B). Therefore, we speculated that magenta module genes were associated with extracellular matrix-related cancer signalling pathways.

3.4 Identification of hub genes by PPI network construction

The genes in magenta module were constructed with a PPI network in Cytoscape software. The MCC algorithm of the plugin cytoHubba plugin was used to calculate the 10 hub genes: COL1A1, COL1A2, FN1, POSTN, THBS2, COL11A1, SPP1, MMP13, COMP, SERPINE1 (Fig. 4A). The correlation heatmap revealed that COL11A1, POSTN, COL1A2, COL1A1, THBS2, and FN1 were highly relevant to each other (Fig. 4B).

Figure 4.

Identification of 10 hub genes. (A), 10 hub genes were identified by cytoscape software. Genes (also named nodes) with deeper colors indicate higher scores by MCC algorithm of cytohubba plugin. (B), the correlation of expression levels among 10 hub genes. The deeper colors indicate higher correlation between each other.

Figure 5.

Expression profiles of ten hub genes based on TCGA and GTEx database on GEPIA platform. Red bar refers to gastric cancer samples ( $N=$ 408), black bar refers to normal gastric samples ( $N=$ 211) for COL1A1 (A), COL1A2 (B), FN1 (C), POSTN (D), THBS2 (E), COL11A1 (F), SPP1 (G), MMP13 (H), COMP (I) and SERPINE1 (J). One-way ANOVA was used for hypothesis test. $|$ log2FC $|$ cutoff $=$ 1 & $p$ -value cutoff $=$ 0.01 were used to evaluate significant expression genes.

Figure 6.

Survival analysis of ten hub genes based on TCGA database on GEPIA platform. The median expression level of each hub gene was as the cut-off value of overall survival for COL1A1 (A), COL1A2 (B), FN1 (C), POSTN (D), THBS2 (E), COL11A1 (F), SPP1 (G), MMP13 (H), COMP (I) and SERPINE1 (J). The Log-rank test was used for hypothesis test. The value of cox proportional hazard ratio (HR) was also calculated. Genes with Logank $p$ value less than 0.05 were considered significantly prognostic-related genes.

3.5 Expression and prognosis analyses of hub genes

To identify the expression levels of 10 hub genes in GC tissues and normal gastric tissues, both the TCGA and GTEx expression profiles were analysed on the GEPIA platform. The results showed that COL1A1, COL1A2, FN1, THBS2, COL11A1, SPP1, COMP and SERPINE1 were significantly more highly expressed in GC tissues than in normal gastric tissues (Fig. 5A–J). The survival analysis suggested that four genes (COL1A1, FN1, POSTN, and SERPINE1) were identified as significant prognosis-related genes (Fig. 6). By intersecting the DEGs and significant prognosis-related genes, we identified three hub genes (COL1A1, FN1, and SERPINE1) associated with the progression and prognosis of GC.

Figure 7.

The relative mRNA levels of FN1, COL1A1 and SERPINE1 in GEO datasets and our experiments with surgical gastric specimens. (A–C), the relative FN1 mRNA levels in GSE54129, GSE118916 and our clinical specimens; (D–F), the relative COL1A1 mNRA levels in GSE54129, GSE118916 and our clinical specimens; (G–I), the relative SERPINE1 mNRA levels in GSE54129, GSE118916 and our clinical specimens. $P<$ 0.05 was considered statistically significant. ${}^{*}$ , $P<$ 0.05, ${}^{**}$ , $P<$ 0.01, ${}^{***}$ , $P<$ 0.001.

Figure 8.

The results of immunohistochemical staining for FN1, COL1A1 and SERPINE1 in our experiments with surgical gastric specimens. (A), the FN1 expression in normal gastric tissues and GC specimens. (B), the COL1A1 expression in normal gastric tissues and GC specimens. (C), the SERPINE1 expression in normal gastric tissues and GC specimens. The objective lens of 100 $\times$ , 200 $\times$ and 400 $\times$ were shown in the results.

3.6 Hub gene validation at the mRNA and protein expression levels

To verify the expression levels of the three hub genes in GC, we used two independent GEO databases and our clinical gastric tissues. The results showed that the relative mRNA levels of FN1, COL1A1 and SERPINE1 were significantly higher in GC in GSE54129 (Fig. 7A, D, G), GSE118916 (Fig. 7B, E, H) and our 20 pairs of clinical gastric specimens (Fig. 7C, F, I). The IHC results showed that FN1, COL1A1 and SERPINE1 were more highly expressed in GC tissues than in normal gastric tissues (Fig. 8A–C). In addition, WB results indicated that FN1, COL1A1 and SERPINE1 were significantly more highly expressed in GC (Fig. 9A–E). Therefore, the results above verified that FN1, COL1A1 and SERPINE1 were three hub genes influencing the development of GC.

3.7 Hub gene validation in prognosis values and diagnostic evaluation

To verify the prognostic value of the three hub genes, we first performed survival analysis of two independent GEO databases at the Kaplan-Meier Plotter online resource. The results showed that overexpression of the three hub genes was all significantly associated with poor prognosis of GC in both GSE15459 (Fig. 10A–C) and GSE62254 (Fig. 10D and E). In addition, univariate and multivariate Cox regression analyses revealed that high expression levels of FN1, COL1A1 and SERPINE1 were independent risk factors for poor prognosis of GC (Table 2).

Table 2
Univariate and multivariate analysis of FN1, COL1A1 and SERPINE1 mRNA expression levels and clinical parameters in TCGA gastric cancer patients

Variables	Univariate anaysis			Multivariate analysis
	HR	95% CI	$P$ value	HR	95% CI	$P$ value
Age (years)	1.63	1.10–2.42	0.015 ${}^{*}$	1.76	1.17–2.62	0.006 ${}^{**}$
$>$ 60
$\leqslant$ 60
Gender	1.49	1.01–2.2	0.043 ${}^{*}$	1.46	0.98–2.17	0.062
Male
Female
Stage	1.79	1.24–2.60	0.002 ${}^{**}$	1.23	0.70–2.19	0.471
III $+$ IV
I $+$ II
Grade	1.39	0.95–2.02	0.088	NA	NA	NA
G3
G1 $+$ G2
T stage	1.65	1.05–2.58	0.03 ${}^{*}$	1.09	0.65–1.85	0.743
T3 $+$ T4
T1 $+$ T2
N stage	1.69	1.10–2.59	0.016 ${}^{*}$	1.42	0.79–2.55	0.237
N1 $\sim$ N3
N0
M stage	1.92	1.04–3.67	0.049 ${}^{*}$	2.00	1.02–3.93	0.043 ${}^{*}$
M1
M0
FN1	1.61	1.13–2.30	0.01 ${}^{*}$	1.59	0.95–2.65	0.076
High
Low
COL1A1	1.15	1.02–1.28	0.02 ${}^{*}$	0.79	0.48–1.30	0.347
High
Low
SERPINE1	1.22	1.10–1.36	$<$ 0.001 ${}^{***}$	1.61	1.07–2.42	0.022 ${}^{*}$
High
Low

Note: ${}^{*}$ , $P<$ 0.05; ${}^{**}$ , $P<$ 0.01. These analyses were conducted using R. High and Low were defined according to the expression median value of each gene.

Figure 9.

The results of western blot for FN1, COL1A1 and SERPINE1 in our experiments with surgical gastric specimens. (A–B), the ECL exposed images of WB in 12 pairs of GC and normal gastric tissues. (D–E), the statistical histograms of WB for FN1, COL1A1 and SERPINE1. ${}^{*}$ , $P<$ 0.05, ${}^{**}$ , $P<$ 0.01, ${}^{***}$ , $P<$ 0.001.

Table 3

Co-enriched Hallmark genesets of COL1A1, FN1 and SERPINE1 by GSEA anlysis

Hallmark names	COL1A1			FN1			SERPINE1
	NES	NOM p-val	FDR q-val	NES	NOM p-val	FDR q-val	NES	NOM p-val	FDR q-val
HYPOXIA	2.507	0.000	0.000	2.363	0.000	0.000	2.387	0.000	0.000
EPITHELIAL_MESENCHYMAL TRANSITION	2.469	0.000	0.000	2.512	0.000	0.000	2.491	0.000	0.000
ANGIOGENESIS	2.444	0.000	0.000	2.245	0.000	0.001	2.393	0.000	0.000
APICAL_JUNCTION	2.417	0.000	0.000	2.330	0.000	0.000	2.210	0.000	0.000
APOPTOSIS	2.398	0.000	0.000	2.101	0.000	0.002	2.238	0.000	0.000
INFLAMMATORY_RESPONSE	2.356	0.000	0.000	2.286	0.000	0.000	2.447	0.000	0.000
MYOGENESIS	2.327	0.000	0.000	2.523	0.000	0.000	2.075	0.000	0.003
TGF_BETA_SIGNALING	2.312	0.000	0.000	2.093	0.000	0.002	2.158	0.000	0.001
COAGULATION	2.299	0.000	0.000	2.373	0.000	0.000	2.490	0.000	0.000
TNFA_SIGNALING_VIA_NFKB	2.298	0.000	0.000	2.122	0.002	0.002	2.622	0.000	0.000
KRAS_SIGNALING_UP	2.222	0.000	0.000	2.066	0.000	0.003	2.404	0.000	0.000
UV_RESPONSE_DN	2.176	0.000	0.000	2.284	0.000	0.000	2.104	0.000	0.002
IL2_STAT5_SIGNALING	2.118	0.000	0.001	2.108	0.000	0.002	2.192	0.000	0.000
COMPLEMENT	2.099	0.000	0.001	2.065	0.002	0.003	2.281	0.000	0.000
IL6_JAK_STAT3_SIGNALING	2.075	0.000	0.001	2.070	0.002	0.003	2.397	0.000	0.000

Note: NES: normalized enrichment score; NOM $p$ -value: nominal p val; FDR q-val: false discovery rate. NES $>$ 2 & FDR $<$ 0.01 as criteria of significant enrichment.

Figure 10.

Survival analysis of FN1, COL1A1 and SERPINE1 based on GSE15459 (A–C) and GSE62254 (D–F) datasets from Kaplan-Meier Plotter platform. Logrank method was used for hypothesis test. $P<$ 0.05 was considered statistically significant.

Figure 11.

ROC curve of FN1, COL1A1 and SERPINE1 mRNA expression levels for prediction of GC. (A), ROC curve of FN1; (B), ROC curve of COL1A1; (C), ROC curve of SERPINE1.

The ROC curve showed that the AUC areas of FN1, COL1A1 and SERPINE1 were 0.702, 0.917 and 0.812, respectively, suggesting that three hub genes might be good diagnostic biomarkers of GC.

3.8 GSEAs for three hub genes

To evaluate the potential carcinogenic mechanisms of FN1, COL1A1 and SERPINE1, we performed GSEA based on hallmark gene sets. The results showed that the majority of significantly enriched pathways among COL1A1, FN1 and SERPINE1 were the same. A total of 15 coenriched pathways were identified, including hypoxia, epithelial mesenchymal transition (EMT), angiogenesis, apical junction, apoptosis, inflammatory response, TGF- $\beta$ signalling, TNF- $\alpha$ signalling via NF- $\kappa$ B, IL-2/STAT5 signalling, and IL-6/JAK/STAT3 signalling (listed in Table 3). These pathways might be involved in the pathogenic process of GC for the FN1, COL1A1 and SERPINE1 genes.

4. Discussion

Identifying potential prognostic and therapeutic biomarkers is very important for cancers [32, 33]. In fact, there have been some studies using WGCNA to explore the correlation between gene modules of gastric cancer and clinical traits [13, 34]. However, as most of the studies use differentially expressed genes as input data of WGCNA, it may destroy genetic diversity and not reflect the true cluster results. In the current study, WGCNA was performed based on the TCGA database at the overall level.

We identified a magenta module that was highly correlated with stage, grade and overall survival. Enrichment analysis showed that the magenta module mainly functioned as extracellular matrix (ECM) components and was related to the ECM-receptor interactions, focal adhesion, and the PI3K-Akt pathway. Increasing evidence has shown that overexpressed ECM proteins such as collagen, fibroblasts, and laminin contribute to promoting cancer cell proliferation, migration and invasion [35, 36, 37]. A recent study reviewed the role of the ECM in different stages of gastric cancer, and indicated that ECM components, receptors and associated signalling molecules should be used as predictive biomarkers or as potential targets in gastric cancer [38].

To identify hub genes of the magenta module, we performed a PPI network [39] and selected ten hub genes for further study. The ten hub genes were highly relevant to each other, indicating that they might have similar functions in GC. Among ten hub genes, eight genes (COL1A1, COL1A2, FN1, THBS2, COL11A1, SPP1, COMP and SERPINE1) were significantly more highly expressed in GC tissues than in normal gastric tissues, among which three DEGs (FN1, COL1A1 and SERPINE1) were significantly associated with poor prognosis in GC.

The COL1A1 gene encodes type I collagen which constitutes a major component of the extracellular matrix. Both GEO and experimental verification from our clinical specimens showed that COL1A1 was more highly expressed in GC tissues than in normal gastric tissues. Univariate Cox regression analyses revealed that high expression levels of COL1A1 were an independent risk factor for poor prognosis of GC. In fact, COL1A1 has been reported to be overexpressed in various tumours such as breast cancer [40], colorectal cancer [41] and hepatocellular cancer [42]. A previous study found that the expression of COL1A1 mRNA was significantly upregulated in GC tissues compared with normal tissues and was related to poor overall survival [43]. Another integrated bioinformatics analysis identified that COL1A1 correlated with the pathogenesis of GC [44]. Together with our results of WGCNA, COL1A1 was considered a crucial oncogene in the progression and prognosis of GC. FN1 encodes fibronectin 1 protein, a glycoprotein of the extracellular matrix that acts as a ligand of the integrin family [45]. In this study, we found that FN1 was highly expressed in GC and was an independent risk factor for GC. Increased FN1 expression levels have been reported in various aggressive tumours and related to the migration and invasion of cancer cells [46, 47, 48]. Although some studies based on bioinformatics or experiments have confirmed that FN1 is overexpressed in GC tissues [49, 50, 51], the potential pathogenic mechanism of FN1 in GC remains unclear. SERPINE1 encodes plasminogen activator inhibitor 1 (PAI-1),which regulates the adhesion balance of cells to ECM [52]. It has been reported that high expression of SERPINE1 may promote invasiveness and aggressiveness and predict poor prognosis in various tumours such as colorectal cancer [53], oral squamous cell carcinoma [54] and breast cancer [55]. Here, we identified SERPINE1 as upregulated in GC tissues and as a prognostic biomarker of GC, which is consistent with some studies based on DNA microarray data [44, 56].

ROC curves are usually used to assess the diagnostic efficiency between two methods. The AUC values of the ROC curves of FN1, COL1A1 and SERPINE1 were all greater than 0.7, suggesting that these genes might be good diagnostic biomarkers of GC.

GSEA revealed that overexpression of three hub genes shared 15 common upregulated biological pathways including hypoxia, epithelial mesenchymal transition, angiogenesis, apical junction, apoptosis, inflammatory response, myogenesis, TGF- $\beta$ signalling, coagulation, TNFA signalling via NF- $\kappa$ B, KRAS signalling, IL-2/STAT5 signalling, and IL6/JAK/STAT3 signalling. These signalling pathways are widely known to be associated with the progression of various cancers, which might give us the direction of further research for GC. The hypoxia pathway has been reported to promote GC malignancy [57], migration and invasion of GC cells [58, 59] and high hypoxia-inducible factor-1 $\alpha$ expression is associated with poor prognosis of GC [60]. Epithelial-mesenchymal transition, angiogenesis, TGF- $\beta$ signalling and IL6/JAK/STAT3 are all involved in the carcinogenic process of GC [61, 62]. In conclusion, we identified and verified that three hub genes, FN1, COL1A1 and SERPINE1 are associated with the carcinogenic progression and poor prognosis of GC, and might be as good diagnostic and prognostic biomarkers of GC.

5. Conclusion

We identified FN1, COL1A1 and SERPINE1 as associated with the progression and poor prognosis of GC. These genes might be good diagnostic and prognostic biomarkers of GC.

Author contributions

ZQY and XJ contributed equally to preparation of the manuscript. ZQY performed the WGCNA analysis and drafted the manuscript. XJ contributed to WB, qPCR and IHC experimental verification. XJL helped collect GC surgical specimens. ZRL, SCH, WH and RJF contributed to interpretation or analysis of data. YLL, SYP and WFF contributed to revision for important intellectual content. SCH and XY contributed to the supervision of the study. All authors read and approved the final manuscript.

Data availability statement

The datasets generated and/or analyzed in the study are available in the National Center of Biotechnology Information’s GEO database (www.ncbi.nlm. nih.gov/gds/) with the accession no. GSE 54129 and GSE118916 and in the Genomic Data Commons Client (https://portal.gdc.cancer.gov/) for TCGA databases.

Footnotes

Acknowledgments

This research was supported by the National Key Research and Development Program of China (No. 2016YFC1302201), National Natural Science Foundation of China (No. 81970502, No. 81860107, No. 81260076, No.81460115) and Leading Talent Training Plan of the Gan-Po Outstanding Talents 555 Project of Jiangxi Province (2010-3-61).

Conflict of interest

The authors declare that they have no conflict of interests.

Supplementary data

Figure S1.

GO and KEGG pathway enrichment analysis of magenta module. (A), Visualization of GO cluster (including BP, CC and MF classification) with a bar chart. BP, biological process; CC, cellular component; MF, molecular function. (B), Visualization of KEGG pathway with a bubble chart.

References

Bray

et al., Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries, CA Cancer J Clin 68(6) (2018), 394–424.

Siegel

R.L.

Miller

K.D.

and Jemal

, Cancer statistics, 2020, CA Cancer J Clin 70(1) (2020), 7–30.

Maruyama

, The most important prognostic factors for gastric cancer patients: a study using univariate and multivariate analyses, Scandinavian Journal of Gastroenterology 22(sup133) (1987), 63–68.

Kim

J.-P.

et al., Significant prognostic factors by multivariate analysis of 3926 gastric cancer patients, World Journal of Surgery 18(6) (1994), 872–877.

Katai

et al., Five-year survival analysis of surgically resected gastric cancer cases in Japan: a retrospective analysis of more than 100,000 patients from the nationwide registry of the Japanese Gastric Cancer Association (2001–2007), Gastric Cancer 21(1) (2018), 144–154.

et al., Comprehensive epigenetic analyses reveal master regulators driving lung metastasis of breast cancer, J Cell Mol Med 23(8) (2019), 5415–5431.

et al., Bioinformatics approaches for anti-cancer drug discovery, Curr Drug Targets 21(1) (2020), 3–17.

Wei

D.Q.

Selvaraj

and Kaushik

A.C.

, Computational perspective on the current state of the methods and new challenges in cancer drug discovery, Curr Pharm Des 24(32) (2018), 3725–3726.

Huang

et al., Circulating miR-1246 targeting UBE2C, TNNI3, TRAIP, UCHL1 genes and key pathways as a potential biomarker for lung adenocarcinoma: integrated biological network analysis, J Pers Med 10(4) (2020).

10.

Cancer Genome Atlas Research

et al., The cancer genome atlas pan-cancer analysis project, Nat Genet 45(10) (2013), 1113–1120.

11.

Barrett

et al., NCBI GEO: archive for high-throughput functional genomic data, Nucleic Acids Research 37(suppl_1) (2008), D885–D890.

12.

Langfelder

and Horvath

, WGCNA: an R package for weighted correlation network analysis, BMC Bioinformatics 9(1) (2008), 559.

13.

Chen

et al., Identification of biomarkers associated with histological grade and prognosis of gastric cancer by co-expression network analysis, Oncol Lett 18(5) (2019), 5499–5507.

14.

Ren

Z.H.

et al., WGCNA co-expression network analysis reveals ILF3-AS1 functions as a CeRNA to regulate PTBP1 expression by sponging miR-29a in gastric cancer, Front Genet 11 (2020), 39.

15.

Zhang

et al., Prognostic value of sorting nexin 10 weak expression in stomach adenocarcinoma revealed by weighted gene co-expression network analysis, World J Gastroenterol 24(43) (2018), 4906–4919.

16.

Nangraj

A.S.

et al., Integrated PPI- and WGCNA-retrieval of hub gene signatures shared between barrett’s esophagus and esophageal adenocarcinoma, Front Pharmacol 11 (2020), 881.

17.

Liu

J.F.

et al., An integrated TCGA pan-cancer clinical data resource to drive high quality survival outcome analytics, Cancer Research 78(13) (2018), 2.

18.

Law

C.W.

et al., voom: precision weights unlock linear model analysis tools for RNA-seq read counts, Genome Biol 15(2) (2014), R29.

19.

Law

et al., RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR [version 3; peer review: 3 approved], 5(1408) (2018).

20.

Chung

et al., Median absolute deviation to improve hit selection for genome-scale RNAi screens, J Biomol Screen 13(2) (2008), 149–158.

21.

Rousseeuw

P.J.

and Croux

, Alternatives to the Median Absolute Deviation, Journal of the American Statistical Association 88 (1993), 1273–1283.

22.

Langfelder

and Horvath

, Fast R functions for robust correlations and hierarchical clustering, J Stat Softw 46(11) (2012).

23.

Kruger

, Editorial change at statistical applications in genetics and molecular biology, Stat Appl Genet Mol Biol 17(4) (2018).

24.

et al., clusterProfiler: an R package for comparing biological themes among gene clusters, OMICS 16(5) (2012), 284–287.

25.

Shannon

et al., Cytoscape: a software environment for integrated models of biomolecular interaction networks, Genome Res 13(11) (2003), 2498–2504.

26.

Chin

C.H.

et al., cytoHubba: identifying hub objects and sub-networks from complex interactome, BMC Syst Biol 8(Suppl 4) (2014), S11.

27.

Tang

et al., GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses, Nucleic Acids Res 45(W1) (2017), W98–W102.

28.

Gyorffy

et al., An online survival analysis tool to rapidly assess the effect of 22,277 genes on breast cancer prognosis using microarray data of 1,809 patients, Breast Cancer Res Treat 123(3) (2010), 725–731.

29.

Zhao

et al., Increased IGFBP7 expression correlates with poor prognosis and immune infiltration in gastric cancer, Journal of Cancer 12(5) (2021), 1343–1355.

30.

Robin

et al., pROC: an open-source package for R and S+ to analyze and compare ROC curves.

31.

Subramanian

et al., Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles, Proc Natl Acad Sci U S A 102(43) (2005), 15545–15550.

32.

Selvaraj

et al., Identification of target gene and prognostic evaluation for lung adenocarcinoma using gene expression meta-analysis, network analysis and neural network algorithms, J Biomed Inform 86 (2018), 120–134.

33.

Selvaraj

et al., Prognostic impact of tissue inhibitor of metalloproteinase-1 in non-small cell lung cancer: systematic review and meta-analysis, Curr Med Chem 26(42) (2019), 7694–7713.

34.

Chen

P.F.

et al., Co-expression network analysis identified CDH11 in association with progression and prognosis in gastric cancer, Onco Targets Ther 11 (2018), 6425–6436.

35.

Patarroyo

Tryggvason

and Virtanen

, Laminin isoforms in tumor invasion, angiogenesis and metastasis, Semin Cancer Biol 12(3) (2002), 197–207.

36.

Naba

et al., The matrisome: in silico definition and in vivo characterization by proteomics of normal and tumor extracellular matrices, Mol Cell Proteomics 11(4) (2012), M111 014647.

37.

Aguilera

K.Y.

et al., Collagen signaling enhances tumor progression after anti-VEGF therapy in a murine model of pancreatic ductal adenocarcinoma, Cancer Res 74(4) (2014), 1032–1044.

38.

Moreira

A.M.

et al., The extracellular matrix: an accomplice in gastric cancer development and progression, Cells 9(2) (2020).

39.

Khan

et al., Identification of novel drug targets for diamond-blackfan anemia based on RPS19 gene mutation using protein-protein interaction network, BMC Syst Biol 12(Suppl 4) (2018), 39.

40.

Liu

et al., Collagen 1A1 (COL1A1) promotes metastasis of breast cancer and is a potential therapeutic target, Discov Med 25(139) (2018), 211–223.

41.

Zhang

et al., COL1A1 promotes metastasis in colorectal cancer by regulating the WNT/PCP pathway, Mol Med Rep 17(4) (2018), 5037–5042.

42.

H.P.

et al., Collagen 1A1 (COL1A1) is a reliable biomarker and putative therapeutic target for hepatocellular carcinogenesis and metastasis, Cancers (Basel) 11(6) (2019).

43.

Ding

and Li

, Identification of COL1A1 and COL1A2 as candidate prognostic factors in gastric cancer, World Journal of Surgical Oncology, 2016, 14.

44.

Liu

et al., Identification of potential key genes associated with the pathogenesis and prognosis of gastric cancer based on integrated bioinformatics analysis, Front Genet 9 (2018), 265.

45.

Pankov

and Yamada

K.M.

, Fibronectin at a glance, J Cell Sci 115(Pt 20) (2002), 3861–3863.

46.

et al., Transcriptional activation of FN1 and IL11 by HMGA2 promotes the malignant behavior of colorectal cancer, Carcinogenesis 37(5) (2016), 511–521.

47.

Sponziello

et al., Fibronectin-1 expression is increased in aggressive thyroid cancer and favors the migration and invasion of cancer cells, Mol Cell Endocrinol 431 (2016), 123–132.

48.

Zhan

et al., Quantitative proteomics analysis of sporadic medullary thyroid cancer reveals FN1 as a potential novel candidate prognostic biomarker, Oncologist 23(12) (2018), 1415–1425.

49.

Suh

Y.S.

et al., Overexpression of plasminogen activator inhibitor-1 in advanced gastric cancer with aggressive lymph node metastasis, Cancer Res Treat 47(4) (2015), 718–726.

50.

Zhang

et al., MicroRNA-200c binding to FN1 suppresses the proliferation, migration and invasion of gastric cancer cells, Biomed Pharmacother 88 (2017), 285–292.

51.

T.P.

et al., Decreased expression of the long non-coding RNA FENDRR is associated with poor prognosis in gastric cancer and FENDRR regulates gastric cancer cell metastasis by affecting fibronectin1 expression, J Hematol Oncol 7 (2014), 63.

52.

Ghosh

A.K.

and Vaughan

D.E.

, PAI-1 in tissue fibrosis, J Cell Physiol 227(2) (2012), 493–507.

53.

Mazzoccoli

et al., ARNTL2 and SERPINE1: potential biomarkers for tumor aggressiveness in colorectal cancer, J Cancer Res Clin Oncol 138(3) (2012), 501–511.

54.

Dhanda

et al., SERPINE1 and SMA expression at the invasive front predict extracapsular spread and survival in oral squamous cell carcinoma, Br J Cancer 111(11) (2014), 2114–2121.

55.

Zhang

et al., Endothelial cells promote triple-negative breast cancer cell metastasis via PAI-1 and CCL5 signaling, FASEB J 32(1) (2018), 276–288.

56.

Nie

et al., Identification of hub genes correlated with the pathogenesis and prognosis of gastric cancer via bioinformatics methods, Minerva Med, 2019.

57.

Liu

et al., Hypoxia promotes gastric cancer malignancy partly through the HIF-1alpha dependent transcriptional activation of the long non-coding RNA GAPLINC, Front Physiol 7 (2016), 420.

58.

X.W.

et al., Hypoxia promotes migration and invasion of gastric cancer cells by activating HIF-1alpha and inhibiting NDRG2 associated signaling pathway, Eur Rev Med Pharmacol Sci 22(23) (2018), 8237–8247.

59.

et al., Hypoxia-inducible microRNA-224 promotes the cell growth, migration and invasion by directly targeting RASSF8 in gastric cancer, Mol Cancer 16(1) (2017), 35.

60.

Zhang

et al., Prognostic value of hypoxia-inducible factor-1 alpha and prolyl 4-hydroxylase beta polypeptide overexpression in gastric cancer, World J Gastroenterol 24(22) (2018), 2381–2391.

61.

Zhao

et al., IL-6 mediates the signal pathway of JAK-STAT3-VEGF-C promoting growth, invasion and lymphangiogenesis in gastric cancer, Oncol Rep 35(3) (2016), 1787–1795.

62.

Zhang

et al., Coexpression of FOXK1 and vimentin promotes EMT, migration, and invasion in gastric cancer cells, J Mol Med (Berl) 97(2) (2019), 163–176.

Weighted correlation network analysis identifies FN1,COL1A1 and SERPINE1 associated with the progression and prognosis of gastric cancer

Abstract

BACKGROUND:

OBJECTIVE:

METHODS:

RESULTS:

CONCLUSION:

Keywords

Abbreviation

1. Introduction

2.1 Data collection from the TCGA and GEO database

2.2 Data process for WGCNA

2.3 Gene co-expression network construction

2.4 Identification of module of interest

2.5 Identification of hub genes with a PPI network

2.6 Differential expression and survival analysis

2.7 Verification of differentially expressed and prognostic genes by GEO databases

2.8 Gastric tissue specimens and immunohistochemistry

Table 1 The primers used in the real-time quantitative PCR analysis

2.10 RNA extraction and real-time quantitative PCR analysis

2.11 Diagnostic efficiency evaluation for GC

2.12 GSEAs of three hub genes based on the TCGA RNA-seq data

3.1 Gene co-expression network and module analysis

3.2 Identification of modules associated with clinical characteristics

3.3 Enrichment analyses of the selected module

3.4 Identification of hub genes by PPI network construction

3.7 Hub gene validation in prognosis values and diagnostic evaluation

Table 2 Univariate and multivariate analysis of FN1, COL1A1 and SERPINE1 mRNA expression levels and clinical parameters in TCGA gastric cancer patients

4. Discussion

5. Conclusion

Author contributions

Data availability statement

Footnotes

Acknowledgments

Conflict of interest

Supplementary data

References

Table 1
The primers used in the real-time quantitative PCR analysis

Table 2
Univariate and multivariate analysis of FN1, COL1A1 and SERPINE1 mRNA expression levels and clinical parameters in TCGA gastric cancer patients