Basic Study Open Access
Copyright ©The Author(s) 2018. Published by Baishideng Publishing Group Inc. All rights reserved.
World J Gastroenterol. Nov 21, 2018; 24(43): 4906-4919
Published online Nov 21, 2018. doi: 10.3748/wjg.v24.i43.4906
Prognostic value of sorting nexin 10 weak expression in stomach adenocarcinoma revealed by weighted gene co-expression network analysis
Jun Zhang, Shuai Guo, Zhe Dong, Zhi-Chao Zheng, Yue Wang, Yan Zhao, Department of Gastric Cancer, Liaoning Cancer Hospital and Institute (Cancer Hospital of China Medical University), Shenyang 110042, Liaoning Province, China
Yue Wu, Department of Emergency, Sheng Jing Hospital of China Medical University, Shenyang 110042, Liaoning Province, China
Hao-Yi Jin, Pancreatic and Thyroid Surgery Department, Sheng Jing Hospital of China Medical University, Shenyang 110042, Liaoning Province, China
ORCID number: Jun Zhang (0000-0003-1739-6462); Yue Wu (0000-0002-2505-7842); Hao-Yi Jin (0000-0002-6685-1102); Shuai Guo (0000-0002-7418-6915); Zhe Dong (0000-0002-2181-8173); Zhi-Chao Zheng (0000-0003-1939-8138); Yue Wang (0000-0003-4536-9904); Yan Zhao (0000-0002-7760-916X).
Author contributions: Zhang J performed the majority of experiments and analyzed the data and drafted the manuscript; Zheng ZC, Zhao Y designed the research; Wu Y, Dong Z conducted the immunohistochemistry assays and assisted in writing the manuscript; Guo S, Wang Y collected and analyzed the data; Zhao Y provided critical revision of the manuscript for important intellectual content; Jin HY provided critical revision of the manuscript for important intellectual content.
Supported by Liaoning S&T Project, No. 2015020269.
Institutional review board statement: The study was reviewed and approved by the Faculty of Science Ethics Committee at Liaoning Cancer Hospital and Institute (Cancer Hospital of China Medical University) (20150308-2).
Conflict-of-interest statement: The authors declare that there are no conflicts of interest related to this study.
Data sharing statement: No additional data are available.
ARRIVE guidelines statement: The authors have read the ARRIVE guidelines, and the manuscript was prepared and revised according to the ARRIVE guidelines.
Open-Access: This article is an open-access article which was selected by an in-house editor and fully peer-reviewed by external reviewers. It is distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited and the use is non-commercial. See: http://creativecommons.org/licenses/by-nc/4.0/
Correspondence to: Yan Zhao, PhD, Professor, Vice-President, Department of Gastric Cancer, Liaoning Cancer Hospital and Institute (Cancer Hospital of China Medical University), No. 44 Xiaoheyan Road, Dadong District, Shenyang 110042, Liaoning Province, China. zhaoyan@cancerhosp-ln-cmu.com
Telephone: +86-24-31916823 Fax: +86-24-24315679
Received: August 2, 2018
Peer-review started: August 2, 2018
First decision: October 5, 2018
Revised: October 17, 2018
Accepted: October 21, 2018
Article in press: October 21, 2018
Published online: November 21, 2018

Abstract
AIM

To detect significant clusters of co-expressed genes associated with tumorigenesis that might help to predict stomach adenocarcinoma (SA) prognosis.

METHODS

The Cancer Genome Atlas database was used to obtain RNA sequences as well as complete clinical data of SA and adjacent normal tissues from patients. Weighted gene co-expression network analysis (WGCNA) was used to investigate the meaningful module along with hub genes. Expression of hub genes was analyzed in 362 paraffin-embedded SA biopsy tissues by immunohistochemical staining. Patients were classified into two groups (according to expression of hub genes): Weak expression and over-expression groups. Correlation of biomarkers with clinicopathological factors indicated patient survival.

RESULTS

Whole genome expression level screening identified 6,231 differentially expressed genes. Twenty-four co-expressed gene modules were identified using WGCNA. Pearson’s correlation analysis showed that the tan module was the most relevant to tumor stage (r = 0.24, P = 7 × 10-6). In addition, we detected sorting nexin (SNX)10 as the hub gene of the tan module. SNX10 expression was linked to T category (P = 0.042, χ2 = 8.708), N category (P = 0.000, χ2 = 18.778), TNM stage (P = 0.001, χ2 = 16.744) as well as tumor differentiation (P = 0.000, χ2 = 251.930). Patients with high SNX10 expression tended to have longer disease-free survival (DFS; 44.97 mo vs 33.85 mo, P = 0.000) as well as overall survival (OS; 49.95 vs 40.84 mo, P = 0.000) in univariate analysis. Multivariate analysis showed that dismal prognosis could be precisely predicted clinicopathologically using SNX10 [DFS: P = 0.014, hazard ratio (HR) = 0.698, 95% confidence interval (CI): 0.524-0.930, OS: P = 0.017, HR = 0.704, 95%CI: 0.528-0.940].

CONCLUSION

This study provides a new technique for screening prognostic biomarkers of SA. Weak expression of SNX10 is linked to poor prognosis, and is a suitable prognostic biomarker of SA.

Key Words: Stomach adenocarcinoma, The Cancer Genome Atlas, Weighted gene co-expression network analysis, Sorting nexin 10, Clinicopathological predictors, Disease-free survival, Overall survival

Core tip: This study used The Cancer Genome Atlas (TCGA) and clinical data to identify sorting nexin (SNX)10 as associated with poor prognosis in stomach adenocarcinoma (SA). Firstly, we downloaded the TCGA-STAD RNA-seq data through R and identified the differentially expressed genes. Weighted gene co-expression network analysis was used to clarify the connection between modules and clinical information. We then chose SNX10 as the hub-gene of the tan-module after considering module membership, genes with high gene significance, degree and Molecular Complex Detection scores. Our center’s SA data were used to evaluate our hypothesis. Our findings demonstrated that weak expression of SNX10 is linked to poor SA prognosis.



INTRODUCTION

In spite of technological progress made in diagnosis and treatment, gastric cancer continues to be the fourth most prevalent cancer type with the second highest mortality[1]. The high occurrence of cancer-related death in stomach adenocarcinoma (SA) patients may be associated with specific biological characteristics of the disease, such as challenges associated with timely diagnosis, late clinical manifestation, and high frequency of invasion and metastasis[2]. Investigation of the underlying tumorigenesis and progression of SA, the most common pathological type of gastric cancer, is a prerequisite for developing new therapeutic agents[3,4]. Although some tumorigenesis and prognostic factors, including genes[5,6] and tumor microenvironment, were evaluated[7], the precise mechanisms and signaling pathways involved are still mostly unknown. Novel molecular biomarkers that can precisely indicate the stage of disease progression and clinical results could help provide personalized treatment and a better understanding of SA pathogenesis.

Rapid technological development of biochips and next-generation sequencing has opened up new dimensions in gene expression research for medical oncology[8]. The Cancer Genome Atlas (TCGA) is a large unified compilation of multi-platform molecular profiles along with clinicopathological annotation data. TCGA offers facilities with structured analysis of potential molecular processes of various clinicopathological parameters linked to cancer[9]. Weighted gene co-expression network analysis (WGCNA) was proposed for establishing co-expression modules or networks of genes to explore the correlation among different groups of genes, or within groups of genes, and clinical features[10]. These modules are built upon large gene expression profiles and the differences between the centrally located genes (hub genes) that drive the crucial cellular signaling pathways within important cell types[11,12]. The WGCNA method provides a functional interpretation tool for systems biology and has led to new insights into the pathophysiology of SA.

In the present study, we used gene expression data from TCGA to create a co-expression network using differentially expressed genes (DEGs) related to SA progression, and defined the gene clusters closely related to SA through WGCNA. Transcriptomics analysis enabled us to identify sorting nexin (SNX)10, a tumor suppressor, which showed over-expression in para-carcinoma tissues and gastric epithelial cells (GSE-1) compared to that in SA tissue and cells. We investigated the correlation between SNX10 expression and the related clinicopathological parameters in SA patients. The prognostic value of SNX10 in SA was also evaluated.

Our research provides a novel and extensive application platform for the identification of SA-related genes, which may be useful to identify novel molecular targets and develop competent therapeutic strategies.

MATERIALS AND METHODS
TCGA gene expression profiles

The RNA sequencing dataset and corresponding clinical records of SA patients were obtained from TCGA (https://portal.gdc.cancer.gov/, Project ID: TCGA-STAD). The RTCGA Toolbox package in the R platform was used to download the raw dataset of RNA sequences[13]. A total of 406 samples were available in TCGA, which included 375 tumor and 31 normal tissue samples. Hub genes were identified from 375 tumor tissue samples (with complete clinical data) using WGCNA.

Screening for DEGs

Differentially expressed RNAs were identified by the package R[14]. The trimmed mean of M values was used to obtain gene read counts with one scaling normalized factor[15]. We calculated the significance of RNAs using the negative binomial model followed by correction of P values by the Benjamini-Hochberg method[16]. Significantly expressed RNAs were identified using adjusted P values of 0.05 and fold changes of 1.

Construction of WGCN

The WGCNA R package was used to establish co-expression networks of genes[10]. Towards this end, 375 tumor samples with expression of 6,232 DEGs and related clinical data were used. The weighted adjacency matrix was created from a power function that was dependent on a soft-thresholding parameter gram, deep split was taken as 4, and we considered the cutoff line as 0.15 for the module dendrogram and merged modules. The network was visualized by Cytoscape 3.5.1[17]. The network module was filtered and studied in Cytoscape using Molecular Compgram, deep split was taken as 4, and we considered the cutoff line as 0.15 for the module dendrogram and merged modules. The network was visualized by Cytoscape 3.5.1[17]. The network module was filtered and studied in Cytoscape using Molecular Complex Detection (MCODE).

Clinically significant module identification and module functional annotation

We identified biologically significant modules using Pearson’s correlation test to evaluate the association between modules and clinical features. Age, tumor stage and tumor grade were chosen as clinical traits. The modules showing the highest correlation with clinical features were chosen. The Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 (https://david.ncifcrf.gov/)[18] was used to categorize the gene data of the target module as follows: biological process FAT, cell component FAT, and molecular function FAT datasets of Gene Ontology (GO) functional[19,20] and pathway enrichment analysis using Kyoto Encyclopedia of Genes and Genomes (KEGG)[21]. P < 0.05 was considered statistically significant.

Identification of hub genes

We calculated the gene significance (GS), high module membership (MM) and MCODE score for different genes in the modules of interest. Based on GS, MM and MCODE score, hub genes were identified.

Ethics statement

The study was sanctioned by the Ethics Committee, Liaoning Cancer Hospital and Institute. Informed consent was obtained from each patient prior to surgery.

Patients and tissue samples

The study cohort comprised 362 patients with SA who underwent gastrectomy at the Liaoning Cancer Hospital and Institute (Shenyang, China) from January 2010 to December 2012. All patients underwent D2 lymph node dissection and gastrectomy at the first consultation. None of the enrolled patients were treated with preoperative chemotherapy or radiotherapy. All patients yielded sufficient paraffin-embedded tumor specimens. Patients with absolute adenocarcinoma, mixed carcinoma or signet ring cell carcinoma were excluded. Tumor staging was done as per the criteria detailed in the 8th edition of the TNM Staging Manual of American Joint Committee on Cancer/Union International Control Center (2017). The average age of patients was 59 years (range 33-80 years), and there were 248 men (68.5%) and 114 women (31.5%). Postoperative adjuvant chemotherapy was administered to patients with staging worse than IIA. The clinicopathological parameters analyzed included: age, gender, Bormann type, age, tumor size, tumor location, Lauren type, perineural invasion, tumor differentiation, T category, invasion of blood vessel, N category, and TNM stage (Table 1).

Table 1 Patient characteristics and univariate analysis, n = 362.
Characteristicn (%)DFS
OS
Time (mo)P valueF valueTime (mo)P valueF value
Age (median, yr)0.1382.2090.3730.795
≥ 60192 (53.0)37.2544.18
< 60170 (47.0)41.2646.29
Gender0.1941.6950.4630.540
Male248 (68.5)37.9444.58
Female114 (31.5)41.7246.46
Bormann type0.00027.3490.00023.230
I28 (7.7)64.0464.04
II161 (44.5)45.9651.32
III167 (46.1)22.8036.72
IV6 (1.7)16.0327.33
Tumor size0.00031.8650.00025.094
≥ 5 cm153 (42.3)30.5938.46
< 5 cm209 (57.7)45.3850.09
Location0.2121.5580.4030.912
Up76 (21.0)34.5542.11
Middle102 (28.2)40.7546.36
Low184 (50.8)40.1345.78
Lauren type0.2411.4280.2141.549
Intestinal166 (45.9)41.8747.77
Mixed carcinoma84 (23.2)39.0744.99
Diffuse112 (30.9)35.6242.07
Tumor differentiation0.00017.7110.00015.214
Moderate and high176 (48.6)44.8449.83
Poor186 (51.4)33.7440.76
Vessel invasion0.00021.8490.00020.983
Yes107 (29.6)29.6737.02
No255 (70.4)43.1048.59
Perineural invasion0.00017.7260.00015.346
Yes88 (24.3)28.1436.38
No274 (75.7)42.6648.00
T category0.00070.1210.00049.574
T168 (18.8)65.7265.75
T224 (6.6)60.2562.00
T312 (3.3)59.0461.33
T4258 (71.3)29.2937.46
N category0.000248.0920.000223.946
N0152 (42.0)61.9363.47
N162 (17.1)36.7348.52
N268 (18.8)24.0733.05
N380 (22.1)10.4917.52
TNM stage0.000386.6230.000231.577
I77 (21.3)64.4765.69
II87 (24.0)62.9964.49
III194 (53.6)19.0729.04
IV4 (1.1)5.509.50
SNX100.00021.8490.00018.278
Over- expression172 (47.5)44.9749.95
Weak-expression190 (52.5)33.8540.84

Follow-up evaluation was conducted every 3-6 mo after surgery for 5 years, which included physical examination, pulmonary, abdominal, and pelvic computed tomography scan, blood count, endoscopy, and liver function tests. The dates of first evidence of recurrence and death were recorded for each case. Survival was estimated as the duration between operation and final follow-up or demise. The final date of investigation was May 1, 2018. The tenure between the resection date and diagnosis of first relapse was considered to be disease-free survival (DFS). Overall survival (OS) was calculated as the time between surgery and death.

Immunohistochemistry

Tissue sections (5 μm thick) were fixed on slides coated with poly-L-lysine. SNX10 antibody (Abcam, 1:150, Cambridge, United Kingdom) was used for immunostaining of sections, and visualized by treatment with suitable biotin-conjugated secondary antibodies. Immmunoperoxidase was used for immunostaining of sections, and treatment with suitable biotin-conjugated secondary antibodies as well as immmunoperoxidase detection involving diaminobenzidine (DAB) substrate (Boster, China) and the Vectastain ABC Elite kit (Linaris, Wertheim, Germany) were used for visualization.

Evaluation of immunohistochemical staining

Immunohistochemical (IHC) staining of SNX10 was independently scored by two pathologists. Staining intensity and the ratio of SNX10-positive cells were measured to obtain the final score. This ratio was enumerated as average percentage of SNX10-positive cells in five different fields of 200× magnification. The positive cells were scored as follows: 0, negative; 1, ≤ 10% positive cells; 2, 10%-50% positive cells; 3, > 50% and ≤ 75% positive cells; and 4, > 75% positive cells. Staining intensity was scored between 0 and 3 (negative, weak, moderate, and strong, respectively). The final score was the product of the two individual scores. Median score of SNX10 = 4 served as the cut-off point. Patients were classified into either the over-expression or weak expression group.

Cell culture

GSE-1 is a human gastric epithelial cell line, and SGC-7901 and BGC-823 are human gastric cancer cell lines. Cells were cultured in RPMI 1640 medium enriched with 10% fetal bovine serum (Invitrogen, Carlsbad, CA, United States), 100 U/mL penicillin, and 100 μg/mL streptomycin (Invitrogen) at a temperature of 37 °C humidified with < 5% CO2. Cells were obtained from China Medical University (Shenyang, China). All experiments were performed in triplicate with independent cell cultures.

Real-time polymerase chain reaction

Total cellular RNA was extracted by TRIzol (Invitrogen) cell separation reagent and the PURELINK miRNA Isolation Kit (Invitrogen). The Promega miRNA First-Strand cDNA Core Kit (Promega, Madison, WI, United States) was used for reverse transcription of 500 ng RNA. Real-time polymerase chain reaction (PCR) (LightCycler 480; Roche AG, Basel, Switzerland) was performed using SYBR Master Mixture (Takara Bio, Kusatsu, Japan). Samples were run in triplicate. U6 DNA served as a loading control. Fold changes in mRNA expression between different cells were determined by 2−ΔΔCT normalization. The following forward and reverse expression primers were used for SNX10: 5’-GAAGGTCACTGTCTGGATTA-3’ and 5’-GTCTGCGTCTGGTTGTAA-3’, respectively.

Western blotting

Sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) (10%-15%) was used to resolve 20 μg cellular proteins and blotted onto polyvinylidende difluoride membranes. After transfer, the membranes were incubated in TBS-Tween buffer supplemented with 20 mmol/L Tris-HCl, 5% nonfat milk, 150 mmol/L NaCl, and 0.05% Tween-20 (pH 7.5) for 1 h at 21 °C, followed by treatment with primary antibodies against SNX10 (1:200; Abcam), and β-actin (1:4000; Santa Cruz Biotechnology, Santa Cruz, CA, United States) overnight at 4 °C, and finally incubated with secondary antibody (1:5000; Santa Cruz Biotechnology). β-Actin served as the endogenous loading control for protein quantity. The gray value of the expression of all proteins was measured using ImageJ (NIH, Bethesda, MD, United States). Data are presented as the mean of three independent observations.

Statistical analysis

Statistical analysis was performed using SPSS version 22.0. One-way analysis of variance was used to compare the mean values between groups. Data from the χ2 test or Fisher’s exact test reflected the correlation between clinicopathological variables and biomarkers (SNX10). Kaplan-Meier survival curves were generated from the data and univariate survival analysis was performed by log-rank test. We tested the effect of prognostic biomarkers on survival rate through the Cox proportional hazards model.

RESULTS
Differentially expressed mRNAs in SA

DEG screening was carried out using the TCGA-STAD dataset including 124375 SA samples and 31 para-carcinoma tissues. R language (Edge R) was used to estimate 6231 expression data sets. From the expression data of 24991 mRNAs, 6231 differentially expressed mRNAs were identified (Figure 1), which were further investigated by WGCNA.

Figure 1
Figure 1 Volcano plot of the differentially expressed mRNAs between stomach adenocarcinoma and para-carcinoma tissues. Red indicates high expression and green indicates low expression (|log2FC| > 1 and adjusted P value < 0.05). The X axis represents a false discovery rate (FDR) value and the Y axis represents a log2FC. Differential mRNAs were calculated by edgeR. There were 2456 highly expressed and 3775 lowly expressed mRNAs. This volcano plot was conducted by the ggplot2 package of R language. FDR: False discovery rate.
WGCN construction

To clarify the functional modules in SA patients, the expression values of 6231 DEGs were included for constructing a co-expression network with the WGCNA package. We selected β = 2 as the soft-thresholding power to construct a scale-free network (Figure 2A and B). Twenty-four distinct co-expression modules were identified that contained 41-1038 genes per module (Figure 2C).

Figure 2
Figure 2 Determination of soft-thresholding power in weighted gene co-expression network analysis and gene clustering dendrograms. A: Analysis of the scale-free fit index for various soft-thresholding powers (β); B: Analysis of the mean connectivity for various soft-thresholding powers; C: Gene clustering tree (dendrogram) obtained by hierarchical clustering of adjacency-based dissimilarity.
Identifying key clinically significant modules

Biologically, the most relevant module has the highest association with clinical features of SA patients. The tan module had the highest association with tumor grade (r = 0.24, P = 7 × 10-6; Figure 3A) and was selected for subsequent analyses. The tan module also had the second-highest association with tumor stage (r = 0.11, P = 0.03; Figure 3A). TOM plot function from the WGCNA package was used to generate a TOM plot (Figure 3B). Eigengene was treated as a representative profile to quantify module similarity. The dendrogram showed the correlation between modules (Figure 3C). The heat map revealed eigengene adjacency of modules (Figure 3D). MM and GS data (Supplementary Table 1) were exported for hub gene identification. Network analysis revealed that expression and regulation of tumor-related genes were altered in SA as opposed to non-tumor stomach tissues.

Figure 3
Figure 3 Identification of modules associated with clinical traits of stomach adenocarcinoma. A: Heat map of the correlation between module eigengenes and clinical traits of stomach adenocarcinoma (SA); B: Topological overlap matrix (TOM) plot; C: Visualization of hierarchical clustering dendrogram of the eigengenes. The eigen gene network represents the relationships among modules; D: Heat map of the eigengene adjacency. The color bars on the left and below indicate the module for each row or column.
Module functional annotation and hub gene identification

All 93 DEGs in the tan module were submitted to the online database DAVID to identify representative GO terms and KEGG pathways[22] for further elucidation of the functional properties of the DEGs. GO analysis showed that the tan module contained biological processes focused on immune response, inflammatory response and defense response to virus infection. These cellular components were focused on the external side of the plasma membrane, extracellular space and cells. Molecular function mainly focused on chemokine activity (Figure 4A, Supplementary Table 2). Figure 4B and Supplementary Table 3 show the most significantly enriched pathways of the tan module following KEGG pathway analysis. The DEGs were enhanced in the cytokine-cytokine receptor interaction, osteoclast differentiation, and Toll-like receptor signaling pathway. GS in the tan module of weighted co-expression was exported to Cytoscape software for WGCN construction (Figure 4C). We sorted the tan module DEGs according to MM, GS and MCODE. We selected the top 30 members of each group and took the intersection (Figure 4D, Supplementary Table 1). SNX10 was identified as the hub gene.

Figure 4
Figure 4 Module functional annotation and hub gene identification. A and B: Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for genes in the tan-module. The X-axis shows the number of genes and the Y-axis shows the GO and KEGG terms; C: Weighted co-expression network for tan-module genes. (nodes are colored according to the degree, and sized according to the Molecular Complex Detection (MCODE) score); D: Hub genes were screened with module membership (MM), gene significance (GS) and MCODE. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; GS: Gene significance; MM: Module membership; MCODE: Molecular Complex Detection.
Expression of SNX10 in SA cells and specimens

Western blotting and real-time PCR showed that SNX10 was weakly expressed in BGC-823 and SGC-7901 cells relative to control GES-1 cells (Figure 5A and B, P < 0.05). We examined levels of SNX10 by immunohistochemistry in 362 SA tissue samples to investigate the role of SNX10 in SA tissue. We found that 195 (53.9%) tissues showed weak SNX10 expression (Figure 5D), and 167 (46.1%) showed high expression (Figure 5E). The corresponding negative results are shown in Figure 5C.

Figure 5
Figure 5 Expression of sorting nexin 10 in stomach adenocarcinoma cells and tissue and its effect on survival. A: Sorting nexin 10 (SNX10) protein expression in stomach adenocarcinoma (SA) cell lines (SGC-7901 and BGC-823) and the human normal gastric epithelial cell line (GSE-1) by western blot. β-actin was used as a loading control; n = 3, aP < 0.05; B: SNX10 mRNA expression in SGC-7901, BGC-823 and GSE-1 by real time-PCR; n = 3, aP < 0.05; C: The representative negative result in SA tissue; magnifications are 200 ×; D: SNX10 was weakly expressed in SA tissue; magnifications are 200 ×; E: SNX10 had strong expression in SA tissue; magnifications are 200 ×; F: Disease free survival curves stratified by SNX10 expression (P = 0.00); G: Overall survival curves stratified by SNX10 expression (P = 0.00). DFS: Disease-free survival; OS: Overall survival; SNX10: Sorting nexin 10.
Correlations between clinicopathological factors and expression of SNX10

Table 2 summarizes the relationship between clinicopathological features and SNX10 expression. According to χ2 or Fisher’s exact test, expression of SNX10 was associated with T category (P = 0.042, χ2 = 8.708), N category (P = 0.000, χ2 = 18.778), TNM stage (P = 0.001, χ2 = 16.744) and tumor differentiation (P = 0.0000, χ2 = 251.930), which was consistent with our WGCNA forecast. The Pearson correlation between SNX10 and tumor differentiation was 0.834 (P = 0.000).

Table 2 Sorting nexin 10 expression and clinicopathologic characteristics, n (%).
CharacteristicSNX10
HighLowP valueχ2 value
172 (47.5)190 (52.5)
Age (median, yr)0.5990.342
≥ 6094 (49.0)98 (51.0)
< 6078 (45.9)92 (54.1)
Gender0.3850.755
Male114 (46.0)134 (54.0)
Female58 (50.9)56 (49.1)
Tumor size0.5660.330
≥ 5 cm70 (45.8)83 (54.2)
< 5 cm102 (48.8)107 (51.2)
Location0.4901.426
Up33 (43.4)43 (56.6)
Middle46 (45.1)56 (54.9)
Low93 (50.5)91 (49.5)
Bormann type0.1675.071
I13 (46.4)15 (53.6)
II85 (52.8)76 (47.2)
III73 (43.7)94 (56.3)
IV1 (16.7)5 (83.3)
Lauren type0.1723.521
Intestinal76 (45.8)90 (54.2)
Mixed carcinoma35 (41.7)49 (58.3)
Diffuse61 (54.5)51 (45.5)
Tumor differentiation0.000251.93
Moderate and high159 (90.3)17 (9.7)
Poor13 (7.0)173 (93.0)
Vessel invasion0.2641.246
Yes46 (43.0)61 (57.0)
No126 (49.4)129 (50.6)
Perineural invasion0.6570.198
Yes40 (45.5)48 (54.5)
No132 (48.2)142 (51.8)
T category0.0428.708
T140 (58.8)28 (41.2)
T211 (45.8)13 (54.2)
T38 (66.7)4 (33.3)
T4113 (43.8)145 (56.2)
N category0.00018.778
N092 (60.5)60 (39.5)
N125 (40.3)37 (59.7)
N228 (41.2)40 (58.8)
N327 (33.8)53 (66.3)
TNM stage0.00116.744
I44 (57.1)33 (42.9)
II53 (60.9)34 (39.1)
III73 (37.6)121 (62.4)
IV2 (50.0)2 (50.0)
Univariate and multivariate analysis of SA patients

Of the 362 patients, 345 had clinical follow-up data but 209 (57.7%) patients died before the end of follow-up. The median DFS was 33.0 mo (range 5-91 mo), and median OS was 45.50 mo (range 7-91 mo). The Kaplan-Meier method and log-rank test were used to study the clinical outcomes and the clinicopathological factors. Patients with weak expression of SNX10 had longer DFS (44.97 vs 33.85 mo, P < 0.05) and OS (49.95 mo vs 40.84 mo, P < 0.05) (Table 1, Figure 5F and G). Clinicopathological parameters with P < 0.05 in univariate analysis were used for multivariate analysis. Cox proportional hazards model was applied using a backward stepwise method. Higher TNM stage predicted shorter DFS and OS after adjusting for Bormann type, tumor size and other factors (DFS: P = 0.000, HR = 78.873, 95%CI: 36.367-171.057; OS: P = 0.000, HR = 76.703, 95%CI: 35.60-165.262; Table 3). Weak expression of SNX10 also served as an independent prognostic factor for prediction of DFS and OS (DFS: P = 0.014, HR = 0.698, 95%CI: 0.524-0.930; OS: P = 0.017, HR = 0.704, 95%CI: 0.528-0.940; Table 3).

Table 3 Multivariate analysis of significant prognostic factors for survival in stomach adenocarcinoma patients.
VariableDFS
OS
P valueHR95%CIP valueHR95%CI
Bormann0.4181.1090.864-1.4240.5311.0830.844-1.388
Tumor size0.3950.8840.665-1.1750.6250.9310.701-1.238
Tumor differentiation0.6860.8740.455-1.6780.6580.8590.440-1.679
Vessel invasion0.2921.4240.738-2.7460.4191.3200.674-2.586
Perineural invasion0.1710.8150.607-1.0930.2690.8640.628-1.139
TNM stage0.00078.87336.367-171.0570.00076.70335.600-165.262
SNX10 expression0.0140.6980.524-0.9300.0170.7040.528-0.940
DISCUSSION

The tumorigenesis and prognosis of SA are dependent on the interaction of several factors, such as overexpression of oncogenes, silencing of tumor suppressor genes, tumor biological characteristics and clinical factors. According to the genetic factors in SA occurrence, development and prognosis, RNA sequencing could further help to study the function of genes at the whole genome level. As a systematic biology method to describe how clinical features associate with genes, WGCNA was applied to investigate co-expression between SA tissue and normal samples. WGCNA, an algorithm for a scale-free network, is widely used in RNA sequencing and profiling microarray analysis[10]. The clustering criteria of WGCNA are focused on biological significance, in contrast to the other clustering algorithms that use geometric distance between data for clustering. A module in WGCNA depicted a group of genes in different tissue samples with similar trends in physiological processes or expression. After module identification, WGCNA calculates the stability of the module and indicates how that correlates with clinical features. The characteristic of a scale-free network is that it highlights several special nodes whose connectivity is significantly greater than that of other general nodes. The special nodes are called hub genes. Thus, it is the correlation between the hub genes of interest and clinical features that can be used to reveal the main causes of disease progression.

In the present study, DEGs were obtained by comparing SA and normal gastric epithelial tissue samples. WGCNA identified 24 distinct co-expression modules from DEGs. The highest positive correlation between the tan module and tumor grade was found through the gene module related to clinical features. Using comprehensive analyses of GS, MM, degree, and MCODE score, we inferred that SNX10 was the hub gene of the tan module. SNX10 is a member of the SNX family, which possesses an evolutionarily conserved PX or phosphoinositide-binding protein with a Phox homology domain. SNX10 is involved in different endosomal transport pathways by targeting endosomal membranes through the PX domain[23,24].

SNX10 only has one PX domain, and is the most basic structured member of the SNX family. Even so, its biological function is still largely unknown. Loss of SNX10 function may be an important therapeutic strategy for inflammatory bowel disease[25]. Chen et al[26] have suggested that SNX10 has interactions with V-ATPase. SNX10 has regulatory functions in membrane trafficking and endosomal stabilization, and SNX10 overdrive can lead to formation and accumulation of giant vacuoles that can be inhibited by brefeldin A[27]. Shen et al[28] have reported that SNX10 is significantly downregulated and acts as a suppressor in human colorectal cancer (CRC). It can increase activity of chaperone-mediated autophagy and decrease expression of p21Cip1/WAF1, thus regulating tumorigenesis and progression. Downregulation of SNX10 leads to loss of control of amino acid metabolism, which activates the mTOR signaling pathway and promotes tumorigenesis of CRC[29].

Gerber et al[30] have described two variants of SNX10 gene as candidates for CRC susceptibility. In addition, SNX10 exhibits the characteristics of a suppressor gene and acts as a potential marker of hepatic carcinoma, which may be regulated by miRNA-30d[31]. These findings show that inhibition of SNX10 expression may be closely related to tumorigenesis and tumor progression, but its clinical significance in SA remains unclear. Using western blotting and real-time PCR, we have revealed differential expression between SNX10 in SA cells and normal gastric epithelial cells shown through bioinformatics analysis. We also confirmed the correlation between SNX10 and tumor grade based on the biopsy database in our center. Hence, this study is believed to be the first to confirm SNX10 as a prognostic factor for DFS and OS in SA patients.

This study had some limitations. First, we did not analyze the effect of SNX10 on the phenotype of SA. The molecular regulatory mechanism of SNX10 also needs to be further explored in future studies. Second, this was a single-center retrospective study. Further studies should be conducted to validate our findings and highlight the mechanistic impact of SX10 on SA tumorigenesis and progression. In addition, sample size needs to be increased for further statistical analysis.

In summary, WGCNA as well as other methods were used to study RNA sequences and clinical data of SA patients from the TCGA database. We conclude that SNX10 is a hub gene associated with tumor grade and acts as an independent prognostic factor for DFS and OS in SA patients. SNX10 has the potential to become a novel prognostic indicator contributing to personalized therapy. However, more experiments are needed to validate the clinical and biological functions of SNX10.

ARTICLE HIGHLIGHTS
Research background

Stomach adenocarcinoma (SA) is by far the most prevalent pathologic version of gastric cancer, whose prognosis is influenced by the complex gene interactions involved in tumor progression. Guidelines have identified the correlation between clinical prognosis and tumor stage and grade. Detection of significant clusters of co-expressed genes or representative biomarkers associated with tumor stage or grade may prompt to highlight the mechanisms of tumorigenesis and tumor progression, and might be helpful to predict SA patient prognosis.

Research motivation

To detect significant clusters of co-expressed genes associated with tumorigenesis, which may help predict SA patient prognosis. The weighted gene co-expression network analysis (WGCNA) method provided a functional interpretation tool for systems biology and led to new insights into the pathophysiology of SA.

Research objectives

The aim of the present study is to reveal a novel biomarker of SA and evaluate the prognostic value of it in SA.

Research methods

The RNA-seq dataset and clinical dataset of SA in The Cancer Genome Atlas (TCGA) were used in this study. The WGCNA was used to identify meaningful modules and hub genes. A 326 patients database was used to evaluate the clinical significance of hub genes via survival analysis.

Research results

Differentially expressed genes (DEGs) (6231) were obtained through whole genome expression level screening. Gene modules (24) were identified using WGCNA, which were observed to be co-expressed. Pearson’s correlation analysis showed the tan-module to be the most relevant to tumor stages. In addition, we detected SNX10 as the hub gene of the tan-module. SNX10 expression was linked to TNM stage and tumor differentiation. Patients with high SNX10 expression trended to have longer disease-free survival (DFS) and overall survival in univariate analysis. Multivariate analysis also showed that dismal prognosis could be precisely predicted clinicopathologically using SNX10. However, more experiments are needed to validate the clinical and biological functions of SNX10.

Research conclusions

WGCNA as well as other methods were employed to study RNA-seq and the clinical data of SA patients from TCGA. SNX10 was considered as a hub gene associated with tumor grade and acted as an independent prognostic factor in SA patient DFS as well as overall survival. It has the potential to become a novel prognostic indicator, thus contributing to personalized therapy.

Research perspectives

Our research team will explore the molecule function of SNX10 by establishing animal models. We will also explore the mutual regulatory mechanism of SNX10 through mass spectrometry analysis and co-immunoprecipitation.

Footnotes

Manuscript source: Unsolicited manuscript

Specialty type: Gastroenterology and hepatology

Country of origin: China

Peer-review report classification

Grade A (Excellent): 0

Grade B (Very good): 0

Grade C (Good): C, C

Grade D (Fair): 0

Grade E (Poor): 0

P- Reviewer: D’Ugo DM, Lim SJ S- Editor: Ma RY L- Editor: Filipodia E- Editor: Bian YN

References
1.  Siegel RL, Miller KD, Jemal A. Cancer Statistics, 2017. CA Cancer J Clin. 2017;67:7-30.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 11065]  [Cited by in F6Publishing: 11897]  [Article Influence: 1699.6]  [Reference Citation Analysis (3)]
2.  Wang SM, Tie J, Wang WL, Hu SJ, Yin JP, Yi XF, Tian ZH, Zhang XY, Li MB, Li ZS. POU2F2-oriented network promotes human gastric cancer metastasis. Gut. 2016;65:1427-1438.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 43]  [Cited by in F6Publishing: 46]  [Article Influence: 5.8]  [Reference Citation Analysis (0)]
3.  Lordick F, Shitara K, Janjigian YY. New agents on the horizon in gastric cancer. Ann Oncol. 2017;28:1767-1775.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 61]  [Cited by in F6Publishing: 67]  [Article Influence: 11.2]  [Reference Citation Analysis (0)]
4.  Corso S, Giordano S. How Can Gastric Cancer Molecular Profiling Guide Future Therapies? Trends Mol Med. 2016;22:534-544.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 50]  [Cited by in F6Publishing: 38]  [Article Influence: 4.8]  [Reference Citation Analysis (0)]
5.  Zhang J, Wu Y, Lin YH, Guo S, Ning PF, Zheng ZC, Wang Y, Zhao Y. Prognostic value of hypoxia-inducible factor-1 alpha and prolyl 4-hydroxylase beta polypeptide overexpression in gastric cancer. World J Gastroenterol. 2018;24:2381-2391.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in CrossRef: 34]  [Cited by in F6Publishing: 37]  [Article Influence: 6.2]  [Reference Citation Analysis (0)]
6.  Fuse N, Kuboki Y, Kuwata T, Nishina T, Kadowaki S, Shinozaki E, Machida N, Yuki S, Ooki A, Kajiura S. Prognostic impact of HER2, EGFR, and c-MET status on overall survival of advanced gastric cancer patients. Gastric Cancer. 2016;19:183-191.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 70]  [Cited by in F6Publishing: 79]  [Article Influence: 9.9]  [Reference Citation Analysis (0)]
7.  Thompson ED, Zahurak M, Murphy A, Cornish T, Cuka N, Abdelfatah E, Yang S, Duncan M, Ahuja N, Taube JM. Patterns of PD-L1 expression and CD8 T cell infiltration in gastric adenocarcinomas and associated immune stroma. Gut. 2017;66:794-801.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 289]  [Cited by in F6Publishing: 332]  [Article Influence: 47.4]  [Reference Citation Analysis (0)]
8.  Vincent AT, Derome N, Boyle B, Culley AI, Charette SJ. Next-generation sequencing (NGS) in the microbiological world: How to make the most of your money. J Microbiol Methods. 2017;138:60-71.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 85]  [Cited by in F6Publishing: 71]  [Article Influence: 8.9]  [Reference Citation Analysis (0)]
9.  Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, Kovatich AJ, Benz CC, Levine DA, Lee AV. An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics. Cell. 2018;173:400-416.e11.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 1998]  [Cited by in F6Publishing: 1852]  [Article Influence: 308.7]  [Reference Citation Analysis (0)]
10.  Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 10254]  [Cited by in F6Publishing: 12973]  [Article Influence: 810.8]  [Reference Citation Analysis (0)]
11.  Botía JA, Vandrovcova J, Forabosco P, Guelfi S, D’Sa K; United Kingdom Brain Expression Consortium, Hardy J, Lewis CM, Ryten M, Weale ME. An additional k-means clustering step improves the biological features of WGCNA gene co-expression networks. BMC Syst Biol. 2017;11:47.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 183]  [Cited by in F6Publishing: 167]  [Article Influence: 23.9]  [Reference Citation Analysis (0)]
12.  Pei G, Chen L, Zhang W. WGCNA Application to Proteomic and Metabolomic Data Analysis. Methods Enzymol. 2017;585:135-158.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 149]  [Cited by in F6Publishing: 188]  [Article Influence: 23.5]  [Reference Citation Analysis (0)]
13.  Samur MK. RTCGAToolbox: a new tool for exporting TCGA Firehose data. PLoS One. 2014;9:e106397.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 107]  [Cited by in F6Publishing: 110]  [Article Influence: 11.0]  [Reference Citation Analysis (0)]
14.  Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139-140.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 22632]  [Cited by in F6Publishing: 24997]  [Article Influence: 1666.5]  [Reference Citation Analysis (0)]
15.  Fang SM, Hu BL, Zhou QZ, Yu QY, Zhang Z. Comparative analysis of the silk gland transcriptomes between the domestic and wild silkworms. BMC Genomics. 2015;16:60.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 57]  [Cited by in F6Publishing: 61]  [Article Influence: 6.8]  [Reference Citation Analysis (0)]
16.  Hochberg Y, Benjamini Y. More powerful procedures for multiple significance testing. Stat Med. 1990;9:811-818.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 1952]  [Cited by in F6Publishing: 1986]  [Article Influence: 58.4]  [Reference Citation Analysis (0)]
17.  Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498-2504.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 24663]  [Cited by in F6Publishing: 28534]  [Article Influence: 1426.7]  [Reference Citation Analysis (0)]
18.  Dennis G Jr, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, Lempicki RA. DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol. 2003;4:P3.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 5554]  [Cited by in F6Publishing: 5653]  [Article Influence: 269.2]  [Reference Citation Analysis (0)]
19.  Gene Ontology Consortium. The Gene Ontology (GO) project in 2006. Nucleic Acids Res. 2006;34:D322-D326.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 671]  [Cited by in F6Publishing: 780]  [Article Influence: 43.3]  [Reference Citation Analysis (0)]
20.  Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25-29.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 26268]  [Cited by in F6Publishing: 26170]  [Article Influence: 1090.4]  [Reference Citation Analysis (1)]
21.  Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27-30.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 18868]  [Cited by in F6Publishing: 20539]  [Article Influence: 855.8]  [Reference Citation Analysis (0)]
22.  Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44-57.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 24735]  [Cited by in F6Publishing: 26370]  [Article Influence: 1758.0]  [Reference Citation Analysis (0)]
23.  Worby CA, Dixon JE. Sorting out the cellular functions of sorting nexins. Nat Rev Mol Cell Biol. 2002;3:919-931.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 317]  [Cited by in F6Publishing: 317]  [Article Influence: 14.4]  [Reference Citation Analysis (0)]
24.  Cullen PJ, Korswagen HC. Sorting nexins provide diversity for retromer-dependent trafficking events. Nat Cell Biol. 2011;14:29-37.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 254]  [Cited by in F6Publishing: 260]  [Article Influence: 20.0]  [Reference Citation Analysis (0)]
25.  You Y, Zhou C, Li D, Cao ZL, Shen W, Li WZ, Zhang S, Hu B, Shen X. Sorting nexin 10 acting as a novel regulator of macrophage polarization mediates inflammatory response in experimental mouse colitis. Sci Rep. 2016;6:20630.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 30]  [Cited by in F6Publishing: 38]  [Article Influence: 4.8]  [Reference Citation Analysis (0)]
26.  Chen Y, Wu B, Xu L, Li H, Xia J, Yin W, Li Z, Shi D, Li S, Lin S. A SNX10/V-ATPase pathway regulates ciliogenesis in vitro and in vivo. Cell Res. 2012;22:333-345.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 54]  [Cited by in F6Publishing: 55]  [Article Influence: 4.2]  [Reference Citation Analysis (0)]
27.  Qin B, He M, Chen X, Pei D. Sorting nexin 10 induces giant vacuoles in mammalian cells. J Biol Chem. 2006;281:36891-36896.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 50]  [Cited by in F6Publishing: 55]  [Article Influence: 3.1]  [Reference Citation Analysis (0)]
28.  Zhang S, Hu B, You Y, Yang Z, Liu L, Tang H, Bao W, Guan Y, Shen X. Sorting nexin 10 acts as a tumor suppressor in tumorigenesis and progression of colorectal cancer through regulating chaperone mediated autophagy degradation of p21Cip1/WAF1. Cancer Lett. 2018;419:116-127.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 25]  [Cited by in F6Publishing: 26]  [Article Influence: 5.2]  [Reference Citation Analysis (0)]
29.  Le Y, Zhang S, Ni J, You Y, Luo K, Yu Y, Shen X. Sorting nexin 10 controls mTOR activation through regulating amino-acid metabolism in colorectal cancer. Cell Death Dis. 2018;9:666.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 14]  [Cited by in F6Publishing: 15]  [Article Influence: 2.5]  [Reference Citation Analysis (0)]
30.  Gerber MM, Hampel H, Zhou XP, Schulz NP, Suhy A, Deveci M, Çatalyürek ÜV, Ewart Toland A. Allele-specific imbalance mapping at human orthologs of mouse susceptibility to colon cancer (Scc) loci. Int J Cancer. 2015;137:2323-2331.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 4]  [Cited by in F6Publishing: 4]  [Article Influence: 0.4]  [Reference Citation Analysis (0)]
31.  Cervantes-Anaya N, Ponciano-Gómez A, López-Álvarez GS, Gonzalez-Reyes C, Hernández-Garcia S, Cabañas-Cortes MA, Garrido-Guerrero JE, Villa-Treviño S. Downregulation of sorting nexin 10 is associated with overexpression of miR-30d during liver cancer progression in rats. Tumour Biol. 2017;39:1010428317695932.  [PubMed]  [DOI]  [Cited in This Article: ]  [Cited by in Crossref: 12]  [Cited by in F6Publishing: 13]  [Article Influence: 1.9]  [Reference Citation Analysis (0)]