Mining the potential prognostic value of synaptosomal-associated protein 25 (SNAP25) in colon cancer based on stromal-immune score

Database and estimation of stromal and immune scores

The TCGA level 3 gene expression data and corresponding clinical information for colon cancer patients were obtained from Genomic Date Commons of the National Cancer Institute (http://portal.gdc.cancer.gov). Only patients with gene expression data, follow up information and clinicopathologic information were included in this study. For normalization, the RNA-seq data of all patients was transformed to transcripts per million (TPM) values (https://pubmed.ncbi.nlm.nih.gov/30379987/). The stromal and immune scores for each TCGA sample were conducted by R 3.6.2 using the R package “estimate” (https://pubmed.ncbi.nlm.nih.gov/24113773/).

We obtained gene expression profiles and clinical information of 430 colon cancer patients from TCGA database. Among them, 231 (53.7%) cases were male and 199 (46.3%) cases were female. The average age of patients at initial pathological diagnosis was 66.3 years (range: 31-90 years). Histologic diagnosis included 369 (85.5%) cases of colon adenocarcinoma and 57 (13.3%) cases of colon mucinous adenocarcinoma, 4 (0.9%) cases were not classified. The tumor stage was stage I in 17.4%, stage II in 38.6%, stage III in 29.1, stage IV in 14.4% of cases, and 2 cases (0.5%) were of unknown stage. The tumors were located in the left (40%) or right (55.3%) colon according to their anatomic neoplasm subdivision, with 4.7% unknown. Based on ESTIMATE algorithm, we obtained stromal scores (range: −2262.07∼1999.52) and immune scores (range: −954.97∼3035.59) for all these colon cancer patients.

Correlations between clinicopathologic data and stromal/immune score

The correlations between clinicopathologic data and stromal/immune score were analyzed by SPSS 22.0 software (SPSS, Inc., Chicago, IL, USA). Patients with colon cancer were divided into high stromal/immune score (the fourth quartile) and low stromal/immune score groups (quartile 1–3). The stromal/immune score of different clinicopathologic groups was compared by Mann–Whitney U test, and overall survival was estimated by the Kaplan–Meier method and compared by log-rank tests. A value of P < 0.05 was considered statistically significant.

Identification of differentially expressed genes (DEGs)

DEGs were identified based on immune and stromal scores (high stromal/immune score group vs. low stromal/immune score group) by Sanger_V1.0.8 software (https://shengxin.ren/softs/Sanger_V1.0.8.zip). Genes with log2 (fold change) >1.5 or < -1 combined with a P value < 0.01 were defined as DEGs. The volcano plot of the DEGs was drawn by Sanger_V1.0.8 software, and the venn diagram was drawn on a website Venny 2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

Gene ontology (GO) and pathway enrichment analyses

Cellular component (CC), molecular function (MF), biological process (BP) and pathway enrichment analyses were conducted using FunRich 3.1.3 (https://pubmed.ncbi.nlm.nih.gov/25921073/). A value of P < 0.05 was considered as the screening condition.

Survival analysis

A COX proportional hazards model was applied to illuminate prognostic DEGs of colon cancer obtained from TCGA and a P value < 0.01 was considered significant. An open source web tool PrognoScan (https://pubmed.ncbi.nlm.nih.gov/19393097/) was conducted to verify the survival outcomes between prognostic DEGs identified and colon cancer patients from GEO database.

Protein–protein interaction (PPI) network construction

The PPI network for DEG-encoded proteins was performed by STRING database (https://pubmed.ncbi.nlm.nih.gov/30476243/) and reconstructed by Cytoscape 3.7.2 (https://pubmed.ncbi.nlm.nih.gov/14597658/). The most significant modular analysis was identified by Molecular Complex Detection (MCODE) plugin of Cytoscape, and the plug-in Biological Networks Gene Oncology (BiNGO) of Cytoscape was applied to analysis GO term of hub genes.

Heatmap and clustering analysis

Heatmap and clustering analysis were completed by “heatmap” package.

Quantitative real-time PCR (qRT-PCR)

Total RNA was extracted with Trizol reagent (TaKaRa Bio Inc. Shiga, Japan) from colon cancer and adjacent normal tissues. cDNA was synthesized with PrimeScript™ RT reagent kit (TaKaRa, RR036A). Quantitative real-time PCR (qRT-PCR) was carried out using the TB Green™ Premix Ex Taq™ kit (TaKaRa, RR420A) on ABI step one Real-Time PCR system. The primers were as follows: AQP4, sense strand 5′- GAGCAGGAATCCTCTATC-3′, antisense strand 5′- AGTGACATCAGTCCGTTT-3′; SNAP25, sense strand 5′- GTAGTGGACGAACGGGAGC-3′, antisense strand 5′- CCTGTCGATCTGGCGATT-3′; GAPDH, sense strand 5′-GTCAACGGATTTGGTCTGTATT-3′, antisense strand 5′- AGTCTTCTGGGTGGCAGTGAT-3′.

The Institutional Medical Ethics Review Board of Taizhou first people’s hospital in Zhejiang Province approval to carry out the study within its facilities (Ethical Application Ref: 2019-KY009-03).

Gene set enrichment analysis (GSEA)

According to the expression level of AQP4 (or SNAP25), samples of the complete cohort from TCGA were divided into 2 groups, and implemented using GSEA by Sanger_V1.0.8 software. The KEGG gene set biological process database (c2.cp.kegg.v6) were chosen for enrichment analysis. Terms with both P value < 0.01 and false discovery rate (FDR) < 0.01 were identified.

Mining the immune-related mechanism of SNAP25

TISIDB (https://pubmed.ncbi.nlm.nih.gov/14597658/) is a friendly web portal integrates 988 immune-related anti-tumor genes from 7 databases (Ru et al., 2019). The correlations between immune features and any gene may be explored in 30 TCGA cancer types. In this study, the TISIDB database was used to investigate the associations between expression (or methylation) of SNAP25 and tumor-infiltrating lymphocytes (TILs). However, the relations between TILs and expression (or methylation) of AQP4 in colon cancer were not integrated in TISIDB database. A value of P < 0.05 was considered statistically significant.

The relationships between stromal/immune score and clinical features

The stromal and immune scores were variously distributed between adenocarcinoma and mucinous adenocarcinoma (Figs. 1A, 1C). Both stromal score and immune score of colon adenocarcinoma cases were significantly lower than those of mucinous adenocarcinoma cases (P = 0.001, 0.011, respectively). In addition, right colon tumors yielded higher immune scores than those left colon cases (Fig. 1D, P < 0.001), though no significant differences between left and right colon were found for the stromal scores (Fig. 1B, P = 0.818).

The potential association of overall survival and stromal/immune score was explored by classifying 430 colon cancer patients into high and low score groups based on their stromal or immune scores. As shown in Figs. 1E, 1F, the median overall survival of patients with a low stromal score was longer than those in high score group (2963 vs. 1930 days, log-rank test P = 0.038); consistently, the median overall survival of patients with a low immune score was longer than those in high score group (2894 vs. 2230 days, log-rank test P = 0.076), although there was no statistically significant difference. However, patients with both a high stromal score and a high immune score were found to have significantly worse survival than those with low scores (1891 vs. 2974 days, log-rank test P = 0.039) (Fig. 1G).

Identification of DEGs by stromal and immune scores in colon cancer

After standardization of the RNA-Seq data for all 430 colon cancer patients obtained from TCGA database, we identified 4881 and 1512 DEGs based on stromal and immune scores, respectively. As shown in the volcano plots of DEGs for stromal/immune score (Figs. 2A, 2B), 4,791 genes were up-regulated and 90 genes were down-regulated for the comparison based on stromal score. Similarly, 1113 genes were up-regulated and 399 genes were down-regulated based on high immune score group vs. low immune score group. Through Venn diagram (Figs. 2C, 2D) analysis, 563 shared up-regulated DEGs and 9 shared down-regulated DEGs from stromal score and immune score groups were identified and selected for subsequent analysis.

Functional and pathway enrichment analyses

Go and pathway enrichment analyses of the above 572 genes were performed (Figs. 2E, 2F, 2G, 2H). For CC, DEGs were mainly associated with plasma membrane, extracellular and extracellular space. With regard to MF, genes were mainly clustered in receptor activity, cell adhesion molecule activity and B cell receptor activity. DEGs in the BP category primarily enriched in immune response, cell communication and signal transduction. The pathway enrichment analysis showed genes were mainly enriched in epithelial-to-mesenchymal transition, peptide ligand-binding receptors and GPCR ligand binding.

Identification of prognostic DEGs in colon caner

The COX proportional hazard regression model was constructed to identify potential prognostic DEGs in colon cancer. Among the 563 shared up-regulated DEGs and 9 shared down-regulated DEGs, 70 up-regulated DEGs associated with poor outcomes were shown in Fig. 3.

PPI network construction and hub gene analysis among prognostic DEGs

To further explore the interplay among the 70 identified prognostic DEGs, we conducted a PPI network containing 30 nodes and 53 edges based on STRING tool and Cytoscape software (Fig. 4A). Module analysis using MCODE was constructed, and 15 hub genes were selected (Fig. 4A). CC, MF, BP analyses of the total 15 hub genes were performed using BiNGO (Fig. 4B). The 15 hub genes were mainly associated with plasma membrane, transporter activity, secretion and channel activity.

Heatmap and clustering analysis of 15 hub genes

The expression level of 15 hub genes in “dead” and “alive” groups was shown in Fig. 4C.

Verifying the survival outcomes of hub genes in the GEO database

The survival outcomes of 15 hub genes were identified on the PrognoScan web tool, which provided overall survival of GSE12945, GSE17536 and GSE17537 datasets for colorectal cancer. Then, AQP4 and SNAP25 were verified (Fig. 5) to be significantly associated with overall survival according to both the log-rank test and COX proportional hazards regression analysis (all P < 0.05).

The expression of two hub genes in 20 colon cancer and adjacent normal tissue sample set using qRT-PCR

The expression of AQP4 and SNAP25 was validated in 20 pairs of clinical tissues using qRT-PCR. Interestingly, the mRNA relative expression levels of both AQP4 and SNAP25 were significantly elevated in colon cancer tissues compared with adjacent normal tissues (P = 0.003, 0.001) (Fig. 6).

GSEA using TCGA database

To further investigate the underlying mechanism of AQP4 and SNAP25 in colon cancer, KEGG pathway enrichment analysis was performed by GSEA. For AQP4, “vascular smooth muscle contraction” (NES = 2.17, P < 0.001, FDR = 0.006) gene set was prominently enriched. For SNAP25, 57 gene sets were enriched, including 11 gene sets which were cancer-related processes (Fig. 7). Besides, high expression of SNAP25 might also be involved in “B cell receptor signaling pathway” (NES = 2.06, P = 0.002, FDR = 0.002), “cell adhesion molecules cams” (NES = 2.05, P = 0.008, FDR = 0.002), “chemokine signaling pathway” (NES = 2.07, P = 0.002, FDR = 0.002), “complement and coagulation cascades” (NES = 2.08, P < 0.001, FDR = 0.002), “T cell receptor signaling pathway” (NES = 2.06, P < 0.001, FDR = 0.002), “adipocytokine signaling pathway” (NES = 1.96, P < 0.001, FDR = 0.005), “aldosterone regulated sodium reabsorption” (NES = 2.02, P < 0.001, FDR = 0.003), “glycosaminoglycan biosynthesis heparan sulfate” (NES = 2.00, P < 0.001, FDR = 0.004), and “insulin signaling pathway” (NES = 1.94, P = 0.002, FDR = 0.007). This suggested that immunity and metabolism may be as well involved in the underlying mechanism of SNAP25 in colon cancer.

Regulation of immune molecules by SNAP25

The spearman’s correlations between lymphocytes and expression, methylation of SNAP25 were performed in TISIDB database (Fig. 8). The associations between the expression of SNAP25 and immune-related signatures of TILs types were shown in Fig. 8A, and the greatest correlations including natural killer cell (NK; r = 0.493, P <2.2e−16), macrophage (r = 0.45, P <2.2e−16), mast cell (r = 0.448, P <2.2e−16), and natural killer T cell (NKT; r = 0.447, P <2.2e−16) were shown in Fig. 8B. The correlations between methylation of SNAP25 and lymphocytes were shown in Fig. 8C, and Fig. 8D displayed the remarkable negative correlations including plasmacytoid dendritic cell (pDC; r =  − 0.419, P = 2.96e−14), type 1 T helper cell (Th1; r =  − 0.406, P = 4.06e−13), T follicular helper cell (Tfh; r =  − 0.391, P = 3.78e−12), and NKT (r =  − 0.391, P = 3.61e−12). Therefore, the potential underlying mechanism of SNAP25 in colon cancer may be involved in the regulation of the above TILs.