Identification and validation of anoikis-associated gene SNCG as a prognostic biomarker in gastric cancer

Identification of hub AG in GC

By setting the cutoff value at |log FC| > 1 and p < 0.05, 67 differentially expressed AGs (DEGs; including 26 downregulated and 41 upregulated AGs) in GC samples were selected from the differential expression analysis (Figure 1A). To identify the AGs associated with prognosis, the analysis of univariate Cox regression as carried out. and Through screening, we identified 20 AGs significantly associated with prognosis (Figure 1B). To further evaluate the functional clusters associated with patients with GC, WGCNA was also performed on the expression profile of 501 downloaded AGs. Based on the scale-free R² = 0.90, the best soft-threshold power was determined as β = 4 (Figure 1C), and three modules, including blue, grey, and turquoise, were identified in WGCNA (Figure 1D). The closest connection between the two clusters was presented as the blue module (r = −0.55, p < 0.05), therefore, we selected the AGs in the blue module for further analysis (Figure 1E). Based on the above analyses, overlapping gene SNCG were ultimately identified as hub AG and used to develop a risk signature in patients with GC (Figure 1F).

Figure 1. Characterization of candidate AGs in GC. (A) Volcano plot of differentially expressed AGs with FDR < 0.05 and |logFC| > 1. (B) Univariate Cox regression analysis showed prognostic AGs with p < 0.01. (C) Soft-thresholding powers scale-free fit index. (D) Clustering dendrogram of GC samples. (E) Heatmap showing the correlation between clinical traits and gene module. Each module was assigned with different colors. The correlation coefficient decreased in size from red to blue. (F) The Venn diagram of genes among DEGs, prognostic genes, and WGCNA lists.

Prognostic and clinical value analyses of SNCG in patients with GC

As a single diagnostic biomarker, the area under the ROC curve (AUC) of SNCG was 0.824, indicating that SNCG had a high predictive accuracy in patients with GC (Figure 2A). When evaluating gene expression level of SNCG in GC, we observed that compared to the tissues from health individuals, the tissues from GC patients exhibited significantly reduced SNCG expression (p < 0.05; Figure 2B). This result was also confirmed at the translational and transcriptional level by western blot (Figure 3A) and qRT-PCR (Figure 3B) analyses, respectively. For external validation, the results of immunohistochemistry obtained from the Human Protein Atlas (HPA) database were also employed (Figure 3C). To further explore the relationship of SNCG expression with the prognosis of GC, the analyses of KM survival and univariate Cox regression were carried out. The results indicated that SNCG was significantly connected to the overall survival (OS) of GC patients (p < 0.05; Figure 2C), meanwhile, the higher expression subgroups of SNCG had a significantly poorer OS (p < 0.001) in both TCGA (Figure 2D) and GSE84437 (Figure 2E) validation cohort. Interestingly, GC patients with age >65 had significantly lower SNCG expression level (Figure 2F). Furthermore, we evaluated the value for the prognosis of GC patients with diverse clinical features in clinics. Significant differences between the SNGC high- or low-expressed patients with age ≤65 years (Figure 2G), M1 stage (Figure 2H), and AJCC stage IV (Figure 2I) were observed. In these two subgroups, all high SNCG-expressed subgroup exhibited dramatically poor OS in comparison to those all low SNGC-expressed subgroup.

Figure 2. Functional analyses of SNCG. (A) ROC curves of SNCG. (B) Gene expression levels of SNCG in GC samples and normal control. (C) Univariate Cox regression of SNCG. Survival curve of SNCG in TCGA (D) and GSE84437 (E) cohorts. (F) Correlations between the expression of SNCG and age. The prognosis of SNCG under the stratifications of age {less than or equal to} 65 (G), M1 subtype (H), and AJCC stage IV (I).

Figure 3. Expression level of SNCG in GC. (A) Western blot analysis of SNCG. (B) QRT-PCR analysis of SNCG. (C) Immunohistochemistry data of SNCG from HPA database.

Functional enrichment analysis

To explore the potential role of SNGC, we carried out the analyses of GSEA and GSVA. GSEA analysis indicated that the SNCG with the highest enrichment score were significantly connected with pathways related to cell cycle and long-term depression (Figure 4A). GSVA analysis also confirmed that the high SNGC-expressed subgroups were significantly enriched in various cancer- and immune-associated pathways, such as mismatch repair, DNA replication, and intestinal immune network for IgA production pathways (Figure 4B), indicating that the activation of SNCG might participate in modulating cancer and immune progression.

Figure 4. Functional enrichment analysis of SNCG. (A) GSEA analysis of SNCG. (B) GSVA analysis of SNCG.

Immune infiltration analysis

To evaluate the proportion of infiltrated immune cell in GC with a threshold of p < 0.05, we applied the CIBERSORT algorithm. As shown in Figure 5A, significantly differential infiltration of several immune cells was discovered between the normal control and GC samples. The results indicated that naïve B cells, M0 macrophages, and M1 macrophages were significantly upregulated in GC samples; however, the proportions of eosinophils, resting mast cells, monocytes, CD4 resting memory T cells, CD8 T cells, plasma cells, and memory B cell, were significantly increased in control samples. Finally, the association of the SNGC expression with infiltration of immune cells in GC was determined. Figure 5B shows that SNCG were strongly associated with the content of immune cells, including Tregs regulator T cells, gamma delta T cells, helper follicular T cells, naïve CD4 T cells, activated memory CD4 T cells, resting NK cells, activated NK cells, monocytes, resting mast cells, M0 macrophages, naïve B cells, and memory B cells, indicating that SNCG may be prognostic targets for GC immunotherapy.

Figure 5. Immune infiltration analysis. (A) The proportion of 22 types of immune cells between normal control and GC samples. (B) Correlation heatmap depicting correlations between infiltrated immune cells and SNCG in GC.

ceRNA interaction networks

To determine the potential underlying mechanisms of SNCG in GC, a ceRNA network centered on SNCG was generated. The results show that 4 pairs of miRNA-mRNA and 3 pairs of lncRNA-miRNA were included and visualized in Figure 6. SNCG was potentially regulated by miRNAs hsa-miR-497-5p, hsa-miR-107, hsa-miR-195-5p, and hsa-miR-103a-3p, while miRNA hsa-miR-103a-3p might be targeted by lncRNAs LA16c-306A4.2, LINC01070, and RP11-34P13.7.

Figure 6. ceRNA interaction network construction. Red represents mRNAs, yellow represents miRNAs, and blue represents lncRNAs.

Risk signature construction and clinical characteristic analysis

According to SNGC expression, a risk signature was generated in both TCGA and validation cohorts. Meanwhile, according to the median of calculated risk scores, we divided the GC patients into high- and low-risk score subgroups. The patients of the high-risk subgroup significantly associated with worse clinical outcomes in both the TCGA (Figure 7A) and validation (Figure 7B) cohorts (p < 0.05). Meanwhile, multivariate and univariate Cox regressions in the TCGA (Figure 7D and 7E) and validation (Figure 7F and 7G) cohorts showed a significant connection of the risk signature to the prognosis of patient, thus serving as an independent prognostic factor in GC. These indicate that the constructed risk signature is significantly related to the prognosis of GC. Additionally, significant association between patients aged ≤65 years with higher risk scores was also found (p < 0.05, Figure 7C).

Figure 7. Efficacy of the constructed model. Survival curve of the TCGA cohort (A) and GSE84437 cohorts (B). (C) Correlations between the risk score and age. Univariate (D) and multivariate (E) Cox regressions of clinicopathological features in TCGA cohort. Univariate (F) and multivariate (G) Cox regressions of clinicopathological features in GSE84437 cohort.

Finally, a nomogram plot was established by combining the risk signature and clinical features to exploit the prognostic value of the risk signature (Figure 8A). The results suggested that our established nomogram had a substantial agreement (Figure 8B), which confirmed that our constructed nomogram could potentially be applied for GC prognosis forecasting.

Figure 8. Development of a nomogram based on risk scores and clinical characteristics. (A) Nomograms to predict the overall survival of GC patients at 1-, 3-, and 5-year using data from TCGA. The red dashed line represented a sample of GC patient’s death probability by year 1, 3, and 5. (B) Calibration plot of the nomogram.

Data collection and processing

Normal gastric tissues and GC samples RNA sequencing datasets TCGA-GC (GC samples = 375, normal samples = 32) annd GSE84437 (GC samples = 433) with reliable sample sources were obtained from the databases of GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA (https://portal.gdc.cancer.gov/). All samples were obtained from Homo sapiens and GSE84437 was selected as the external validation. A total of 501 anoikis-associated genes (AGs) with a relevance score >0.4 were obtained from the GeneCards database. All public databases in this study were searched following the relevant guidelines, and no ethical approval was required from the Ethics Committee of the First People’s Hospital of Linping District.

Screening for overlapping hub AGs

To improve the highly connected genes discriminatory capacity, univariate Cox regression analysis, WGCNA, and differential expression analysis, were carried out to identify the differentially expressed genes associated with anoikis and were significantly related to GC prognosis. Differential expression analysis for AGs in control and GC samples were conducted using the R package “limma” with the criteria |log fold change (FC)| > 1 and p < 0.05. The packages “ggplot2” and “pheatmap” were applied to generate a volcano plot and heatmap, respectively. Meanwhile, to identify the AGs associated with prognosis, the analysis of univariate Cox regression with the cutoff criterion of p < 0.01 was conducted. The co-expression network between AGs and sample modules was constructed with the R package “WGCNA”. The best soft threshold was determined to be five as it met the minimum power value with a scale-free topological criterion of >0.90. Thereafter, a topological overlap matrix (TOM) transformed the weighted adjacency matrix, and dendrograms of TOM were constructed using the hierarchical clustering method. To avoid the generation of excessive modules, the major parameters were set as a minModuleSize of 50, deepSplit of two, and height cut-off of 0.25. Finally, AGs with high interconnections were classified into different patterns. The gene significance (GS) and module eigengenes (MEs) were calculated for each module. AGs overlapping between the analyses of WGCN, differential expression, and univariate Cox regression were visualized with package “VennDiagram” and considered as hub AGs for further analysis.

Clinical analyses of hub AG

The clinical value of hub AG for GC, as well as their specificity and sensitivity, were determined using Kaplan-Meier curves, univariate Cox regression, differential expression, and receiver operating characteristic (ROC) curve methods using various R packages. Based on the identified hub AG expression, GC patients were divided into two subgroups with high- and low-expression. Subsequently, the prognostic value of hub AG for GC patients with different clinical features were tested and visualized by the “survival” package in R. Statistical significance was set at p < 0.05. The HPA database was applied to determine the difference of hub AG at protein level.

Gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA)

GSVA was performed for each gene set and scoring. According to the GSVA score matrix, the changes at the gene-level were converted into changes at the pathway-level by the R package “GSVA”, and the potential biological functions were ultimately evaluated. The C7 and Hallmark gene sets v7.4 were applied in the analysis of GSEA. Enriched gene sets were used to detect KEGG pathways. Gene sets with the adjusted p-value less than 0.05 were considered significantly enriched after 1000 substitutions.

Immune cell infiltration assessment

Using the CIBERSORT algorithm, the proportion of 22 types of immune cell infiltration was determined on TCGA-GC sequencing data. Using the package “vioplot”, we compared the difference of the infiltration of immune cells between GC and normal gastric tissues. Furthermore, the connection of the hub gene expression to the infiltration of immune cells was also quantified based on CIBERSORT analysis, while a correlation heatmap was visualized using the R package “corrplot”.

ceRNA interactions network construction

Interactions between hub AGs and miRNAs were predicted using miRWalk3.0 (http://mirwalk.umm.uni-heidelberg.de/). Meanwhile, the potential target lncRNAs of identified miRNAs were predicted based on the lncRNASNP3 (http://gong_lab.hzau.edu.cn/lncRNASNP3/#!/). The interactions that were fitted to the TargetScan, miRDB, and miRanda databases were selected to draw the ceRNA networks. Network of ceRNA was constructed and visualized using Cytoscape (version 3.7.1) software.

Identification of hub AG-associated risk signature

Based on the identified hub AG expression, we calculated the risk score of the hub AG-associated risk signature with the formula:

$$\text{risk score}=\sum \text{expgenei}\times \text{βi}$$

where i and β are regression coefficients, and expgenei is hub AG relative expression.

Based on the calculated risk scores, we separated the GC patients into two subgroups with high- and low-risk scores. To compare the ability of prognosis according to the risk level, the analyses of Cox regression and survival were carried out. For predicting the outcome of GC patients, a nomogram was constructed according to the calculated risk score and GC clinical characteristics, such as age, T stage, AJCC stage, gender, grade, N stage, and M stage, using the package “rms”. To validate our constructed nomogram, the bootstraps with 1000 resamples were employed. The discrimination and accuracy were estimated using calibration curves.

Cell culture

The GES-1 human gastric mucosal epithelial cells, MGC803 and SGC7901 GC cells were provided by Shanghai Zhong Qiao Xin Zhou Biotechnology Co., Ltd. The FBS (10%)- and penicillin/streptomycin (1%)-contained RPMI-1640 medium (Gibco) was employed to culture the cells at 37°C under 5% CO₂.

qRT-PCR

To isolate the RNA from cells, Trizol reagent (Invitrogen) was employed. The RNeasy mini kit (QIAGEN, Valencia, CA, USA) was used to isolate RNA and then RNA was reverse transcribed to cDNA using a Transcriptor First-strand cDNA synthesis kit (Roche). Next, a qRT-PCR was conducted using SYBR green master mix on an Applied Biosystems 7500 Real Time Cycler (Applied Biosystems). The PCR conditions are as follows: 95°C for 10 min, 40 cycles at 95°C for 15 sec and 60°C for 1 min, followed by a standard melting curve. All samples were running in triplicate and gene expression was measured using the 2^−ΔΔCT method. GAPDH was used as a normalization control. A list of the primer pairs and their sequences that were used in this study is presented in Supplementary Table 1.

Western blot analysis

To extract the total protein from cells, phenylmethanesulfonyl fluoride (Beyotime)-contained RIPA Lysis Bufer (Beyotime) was employed. For each sample, 60 μg proteins were separated using 10% SDS-PAGE gel. After transferring of the protein bands to a PVDF membrane, overnight incubation 4°C of the membrane with antibodies against SNGC (Abcam) and GAPDH (Abcam) was performed. Next, after 3-time washing using PBST and 1-hour incubation with DyLightTM800 4X PEG conjugated secondary antibody (Cell Signal Technology), the protein bands were identified by An Odyssey Infrared Imaging System.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.