This article is mentioned in:

Introduction

Gastric cancer (GC) is one of the leading causes of cancer-related mortality worldwide, and is particularly prevalent in East Asian countries, including China, Japan and Korea (1). Each year, ~990,000 people are diagnosed with GC worldwide; ~738,000 of whom succumb to the disease (2). The high patient mortality rate is due to the fact that the clinical manifestations of GC usually only become apparent at an advanced disease stage, when the current available therapies will have a limited effect (3,4). Therefore, it is of utmost importance to understand the associated mechanisms and to identify biomarkers for the development of strategies for the screening, early detection and treatment of GC.

GC is a complicated and multifactorial disease, and environmental and genetic factors play important roles in its etiology (5). One of the characteristics of gastric malignant cells is metastasis, whereby cancer cells penetrate vascular channels and invade parenchymal tissue to form satellite tumors in distant organs (6). In this process, the extracellular matrix (ECM) and the basement membrane provide a protective barrier to prevent cancer cell invasion and metastasis (7). Similar to other malignancies, gene expression profiling using complementary DNA microarrays has been used to identify genes involved in gastric carcinogenesis, and to identify novel diagnostic and prognostic markers for GC (8–11). Recent studies have reported genetic alterations in GC, involving tumor suppressor genes, cell adhesion molecules, oncogenes and growth factors, such as p53, trefoil factor 1 and E-cadherin (10,12–15). However, these studies have yielded few useful biomarkers, most likely due to shortcomings concerning the experimental design, the validity of the supporting statistical analysis and the gene selection in the studies. Thus, the present study focused on the gene expression profiling of GC to identify novel biomarkers in this disease.

With the same gene expression profile, Wang et al performed gene set enrichment analysis and identified that increased INHBA expression was associated with poor survival in GC (16). A study by Liu et al demonstrated that the ECM-receptor and cell cycle pathways may play important roles in GC (17). In addition, a study using the same microarray data revealed high periostin expression in GC tissues, which was associated with gene groups that regulated the cell proliferation and cell cycle (18). The present study analyzed the differentially-expressed genes (DEGs) in GC using gene expression profiling. Comprehensive bioinformatics was used to analyze the significant pathways and functions, and to construct the gene co-expression network and sub-network to investigate the critical DEGs of GC. The study aimed to obtain a better understanding of the molecular circuitry in GC and to identify genes potentially useful as novel diagnostic or therapeutic markers for GC.

Materials and methods

Affymetrix microarray data

The gene expression profile of GSE19826 (16) was downloaded from the Gene Expression Omnibus database (19), which freely distributes high-throughput molecular abundance data, largely gene expression data generated by microarray technology. The platform information is as follows: GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix Inc., Santa Clara, CA, USA). In this dataset, 12 cancerous portions of gastric specimens (from Chinese patients) and 15 normal gastric tissues (controls) were included.

Data preprocessing and screening of DEGs

The preprocessed microarray data were obtained and then log2 transformation was performed on these data. The most popular method, the Linear Models for Microarray data (limma) package (20) in R language (21), was used to analyze the chip data. Upregulated and downregulated genes were identified between GC and normal controls. The false discovery rate (FDR) (22) was utilized for multiple testing correction using the Benjamini and Hochberg method (23). The threshold for the DEGs was set as |log2 fold change (FC)|>1.5 and FDR <0.05.

Hierarchical clustering

Hierarchical clustering methodology is a powerful data mining approach that has been extensively applied to identify groups of similarly expressed genes or conditions from gene expression data. In order to reveal sets of samples in which the closest groups were adjacent, two-way hierarchical clustering analysis (24) was performed on genes and conditions using Euclidean distance (25) by the ‘pheatmap’ package (http://cran.r-project.org/web/packages/pheatmap/index.html) in R language. The result was represented by a heatmap.

Co-expression network construction of DEGs

From the perspective of systems biology, functionally-related genes are frequently co-expressed across a set of samples (26). COXPRESdb (http://coxpresdb.jp) provides co-expression associations for multiple species of mammals, as comparisons of co-expressed gene lists can increase the reliability of gene co-expression determinations (27). The gene co-expression network was constructed to assess the functional associations between co-expressed genes of DEGs using COXPRESdb, in which genes were indexed by their Entrez Gene IDs. To obtain the co-expression associations, a Pearson Correlation Coefficient >0.6 was chosen as the threshold.

Selection of modules in co-expression network

Gene products in the same module often have the same or similar functions, and they work together to perform one bio-function (28). Therefore, the network was visualized using Cytoscape (29) and module division was made by using the plugin ClusterOne (30) in Cytoscape (parameters: Minimum size, 3; overlap threshold, 0.8), then module function was annotated using another plugin-Bingo (31) and the significant function of each module was achieved.

Function and pathway enrichment analysis of DEGs in modules

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed for the DEGs in the co-expression network using the online tool, DAVID (32). P<0.05 was used to indicate statistical significance.

Results

DEG screening

Following data preprocessing, 42,450 genes were mapped to the probes; the gene expression profile after normalization is shown in Fig. 1. The black lines in each of the boxes, representing the medians of each dataset, are almost in a straight line, indicating a good degree of standardization. Compared with the normal tissues, a total of 596 genes were differentially expressed, consisting of 182 upregulated and 414 downregulated genes.

Figure 1. Normalized expressed value data. The black line in each box represents the median of each set of data, which determines the degree of standardization of data through its distribution.

Hierarchical clustering

The hierarchical clustering algorithm was used to group the genes and samples on the basis of similarities of gene expression. In the results shown in Fig. 2, one normal sample was grouped into the region of GC samples, suggesting that 93.33% of samples were classified correctly. Thus, the DEGs screened had significant expression patterns that could distinguish the disease samples from the normal controls.

Figure 2. Heat-map overview of the two-way hierarchical clustering analysis of the differentially-expressed genes. Blue coloration indicates decreased expression in GC and red coloration indicates increased expression in GC. Samples (12 GC and 15 controls) are represented by columns. The red box represents the normal sample that was enriched into the cancer samples group.

Co-expression network construction and module selection

A total of 57 co-expressed gene pairs were determined between DEGs. In Fig. 3, the upregulated and downregulated genes tended to connect up, respectively. The study identified 4 modules from the network (Fig. 4). The function of the DEGs in each module is presented in Table I. Module 1 had one upregulated gene, PDZ and LIM domain 7 (PDLIM7). Genes in modules 2 and 4, which mostly belonged to the collagen (COL) family, were significantly associated with phosphate transport.

Figure 3. Co-expression network of differentially-expressed genes. Green nodes represent downregulated genes. Red nodes indicate upregulated genes.

Figure 4. Selected modules from the gene co-expression network. Green nodes represent downregulated genes. Red nodes indicate upregulated genes.

Table I. Functions of the genes in the modules.

Table I.

Functions of the genes in the modules.

GO-ID	corr P-value	N	Description	Genes in test set
Module 1
31214	4.30×10−2	2	Biomineral formation	CASR, PDLIM7
1503	4.30×10−2	2	Ossification	CASR, PDLIM7
Module 2
6817	8.57×10−8	5	Phosphate transport	COL6A3, COL1A2, COL1A1, COL5A2, COL5A1
15698	7.47×10−7	5	Inorganic anion transport	COL6A3, COL1A2, COL1A1, COL5A2, COL5A1
6820	1.25×10−6	5	Anion transport	COL6A3, COL1A2, COL1A1, COL5A2, COL5A1
6811	8.30×10−4	5	Ion transport	COL6A3, COL1A2, COL1A1, COL5A2, COL5A1
6810	3.10×10−2	5	Transport	COL6A3, COL1A2, COL1A1, COL5A2, COL5A1
48513	1.76×10−4	6	Organ development	FAP, COL6A3, FBN1, COL1A2, COL1A1, COL5A2
48731	9.03×10−4	6	System development	FAP, COL6A3, FBN1, COL1A2, COL1A1, COL5A2
48856	1.64×10−3	6	Anatomical structure development	FAP, COL6A3, FBN1, COL1A2, COL1A1, COL5A2
7275	4.45×10−3	6	Multicellular organismal development	FAP, COL6A3, FBN1, COL1A2, COL1A1, COL5A2
51234	1.03×10−2	6	Establishment of localization	FAP, COL6A3, COL1A2, COL1A1, COL5A2, COL5A1
32502	1.29×10−2	6	Developmental process	FAP, COL6A3, FBN1, COL1A2, COL1A1, COL5A2
51179	1.62×10−2	6	Localization	FAP, COL6A3, COL1A2, COL1A1, COL5A2, COL5A1
32501	2.24×10−2	6	Multicellular organismal process	FAP, COL6A3, FBN1, COL1A2, COL1A1, COL5A2
Module 3
6936	2.46×10−3	2	Muscle contraction	MYL1, KBTBD10
3012	2.46×10−3	2	Muscle system process	MYL1, KBTBD10
Module 4
6817	7.88×10−4	3	Phosphate transport	COL6A3, COL12A1, COL1A1
15698	2.15×10−3	3	Inorganic anion transport	COL6A3, COL12A1, COL1A1
6820	2.47×10−3	3	Anion transport	COL6A3, COL12A1, COL1A1
22610	2.64×10−2	3	Biological adhesion	COL6A3, COL12A1, THBS2
7155	2.64×10−2	3	Cell adhesion	COL6A3, COL12A1, THBS2
6811	3.28×10−2	3	Ion transport	COL6A3, COL12A1, COL1A1
48513	1.16×10−2	4	Organ development	FAP, COL6A3, COL12A1, COL1A1
48731	2.64×10−2	4	System development	FAP, COL6A3, COL12A1, COL1A1
48856	3.20×10−2	4	Anatomical structure development	FAP, COL6A3, COL12A1, COL1A1
7275	3.46×10−2	4	Multicellular organismal development	FAP, COL6A3, COL12A1, COL1A1
51234	4.04×10−2	4	Establishment of localization	FAP, COL6A3, COL12A1, COL1A1
32502	4.11×10−2	4	Developmental process	FAP, COL6A3, COL12A1, COL1A1
51179	4.52×10−2	4	Localization	FAP, COL6A3, COL12A1, COL1A1

[i] Corr P-value, corrected P-value; N, nodes; GO-ID, gene ontology identification.

Pathway annotation of the DEGs in modules

Two pathways were found to be enriched (Table II), the ECM-receptor interaction and focal adhesion pathways. Of these, the ECM-receptor interaction pathway was most significantly enriched (P=9.44×10−5), and 7 genes [collagen, type VI, α3 (COL6A3), COL1A2, COL1A1, COL5A2, thrombospondin 2 (THBS2), COL11A1 and COL5A1] were predicted to participate in the pathway.

Table II. Significant pathways of DEGs in the selected modules.

Table II.

Significant pathways of DEGs in the selected modules.

Term	Count	FDR	Genes
hsa04512: ECM-receptor interaction	7	9.44×10−5	COL6A3, COL1A2, COL1A1, COL5A2, THBS2, COL11A1, COL5A1
hsa04510: Focal adhesion	8	9.10×10−4	PGF, COL6A3, COL1A2, COL1A1, COL5A2, THBS2, COL11A1, COL5A1

[i] Term represents the pathway name. Count represents the number of DEGs enriched in each pathway. DEGs, differentially-expressed genes; FDR, false discovery rate; ECM, extracellular matrix.

Discussion

GC is the fourth most frequently occurring malignant tumor worldwide, with high incidence and mortality rates. Therefore, it is of great importance to conduct research on the treatment of GC (33). Major efforts are being made to understand GC at a molecular level (34). Since microarrays can simultaneously investigate the expression levels of thousands of genes in the human genome, use of the technique has been widely applied in the identification of disease biomarkers (26,35). In the present study, a total of 596 DEGs were identified in the GC samples compared with the normal controls. Furthermore, the co-expression interaction network of DEGs was construction and 4 modules were identified. The upregulated PDLIM7 gene was enriched in module 1, while 7 other upregulated genes (COL6A3, COL1A2, COL1A1, COL5A2, THBS2, COL11A1 and COL5A1) were involved in the ECM-receptor interaction pathway. The COL family of genes were mainly enriched in module 2, for which the function was phosphate transport.

COL1A1 and COL1A2 encode the α1 and α2 chains of type I collagen, respectively (36). Collagen is the main constituent of the ECM component in tumors, and a number of collagen types have been found in GC tissues (37). The major constituents of the ECM are collagens, adhesive glycoproteins and proteoglycans (38). Specific interactions between cells and ECM-mediated cell-surface-associated components and transmembrane molecules result in the control of cellular activities, such as adhesion and migration (39). Matsui et al showed that collagen degradation, which was an essential step in the tumor cell invasion of the surrounding tissues, was increased in GC tissues (40). Su et al reported that COL1A1 and COL1A2 were commonly upregulated in GC, and were associated with invasion and metastasis (8). In line with this previous study, the present results showed that COL1A1 and COL1A2 were upregulated in GC, suggesting that they play an important role in cancer cell invasion and metastasis in this disease. On the other hand, the COL family genes were mainly enriched in modules 2 and 4, for which the function is phosphate transport. COL6A3 was clustered into module 2. COL6A3 encodes one of the three α chains of type VI collagen. Another significant DEG that was enriched in GC was COL11A1, another member of the COL family, which encodes one of the two α chains of type XI collagen. Using microarray technology, COL6A3 and COL11A1 levels have been proven to be elevated in GC endothelium when compared with normal endothelium (41,42). The present study demonstrated that COL6A3 and COL11A1 were upregulated and participated in the ECM-receptor interaction pathway, which was in line with these previous studies. Taken together, the results indicated that the COL family in the present study may be molecular biomarkers for GC.

THBS2, which has demonstrated functions as a potent inhibitor of tumor growth and angiogenesis, is a disulfide-linked homotrimetric glycoprotein that mediates cell-to-matrix and cell-cell interactions (43). Stamper et al reported that genes associated with ECM-receptor interactions, including TBHS2, underwent significant changes in expression when comparing craniosynostosis patients and controls (44). In addition, Yasui et al suggested that changes in the ECM could be induced by the degradation of collagen I, which was of great importance to the infiltration and metastasis of cancer cells in GC (45). In line with this previous study, the present results also indicated that TBHS2 was upregulated in GC compared with normal controls, suggesting that TBHS2 may play a role in ECM changes and promote GC progression.

PDLIM7 is a family of proteins composed of PDZ and LIM domains that have been proposed to direct protein-protein interactions. Wu et al demonstrated that the LIM domains of Enigma recognized tyrosine-containing motifs with specificity residing in the target structures and the LIM domains (46). Another study showed that receptor tyrosine kinases play essential roles in the control of cancer cell growth and differentiation (47). In the present study, PDLIM7 was found to be upregulated, showing enrichment in module 1, and interacted with other DEGs identified in the study. Another hub gene in module 1 was adenosine deaminase, RNA-specific, B1 (ADARB1) (Fig. 4). ADARB1, also known as ADAR2, encodes the enzyme responsible for pre-mRNA editing of the glutamate receptor subunit B by site-specific deamination of adenosines (48). A previous study demonstrated that the dysregulation of adenosine to inosine in human cancers possibly contributed to the altered transcriptional program required to sustain carcinogenesis (49). Moreover, Camarata et al reported that PDLIM7 could regulate T-box protein 5 transcriptional activity, which is involved in the transcriptional regulation of genes required for mesoderm differentiation (50). In this context, we speculate that PDLIM7 may play a crucial role in GC development via the interaction with ADARB1.

In conclusion, the present study investigated the critical genes in GC based on microarray data. The target genes COL1A1, COL1A2, COL6A3, THBS2, COL11A1, PDLIM7 and ADARB1 were involved in the progression of GC. COL6A3, COL1A2, COL1A1, THBS2 and COL11A1were identified to be involved in the ECM-receptor interaction pathway. Furthermore, the genes of the COL family were associated with phosphate transport. COL1A1 and COL1A2 may play an important role in tumor invasion and metastasis in GC. TBHS2 may impact ECM changes and promote GC progression. Moreover, PDLIM7 may play a crucial role in GC development via the interaction with ADARB1. The genes identified in GC tissues in the present study may prove to be molecular biomarkers for this disease, although further studies must be performed to confirm these results.

Acknowledgements

The authors wish to express their warm thanks to Fenghe (Shanghai) Information Technology Co., Ltd. (Shanghai, China), whose ideas and assistance provided a valuable added dimension to this study.

References