Gene expression data and corresponding clinical information for GC, which included 375 tumor tissues and 32 normal tissues, were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Table 1, Table 2 shows the characteristics of patients with GC from the TCGA database. The gene expression of 174 normal tissues was downloaded from GTEx through UCSC XENA (http://xena.ucsc.edu). Fragments Per kilobase per Million (FPKM) RNAseq data were converted into transcripts Per Million reads (TPM), and log2 translated for subsequent analysis. All tissue samples with incomplete clinical data were excluded.

We downloaded TPM RNAseq data for tumor tissues (TCGA) and normal tissues (TCGA and GTEx) from the UCSC XENA. The differential expression between tumor and normal tissues was tested by Wilcoxon Rank Sum Test and visualized through boxplots and scatter plots. We also used Receiver Operating Characteristic (ROC) curve to determine the diagnostic value of BGN gene expression for GC.

Tumor and para-cancer biopsy tissues were collected from 12 consecutive patients that were diagnosed with GC for the first time from the Endoscopy Center of the First Affiliated Hospital of Guangxi Medical University. The body tissues were immediately immersed in RNA protection solution and rapidly stored in a refrigerator at −80°C. No patient was diagnosed with any other malignancy, nor had they received any treatment for the tumor.

Total RNA of tissues was extracted using Trizol reagent (R0016, Beyotime Biotechnology Co., Ltd., Shanghai, China, according to the manufacturer’s instructions. Complementary DNAs(cDNAs) were generated from 1 µg RNA PrimeScript™ RT Reagent Kit with gDNA Eraser (RR047A, Takara Bio, Inc.). RT-PCR was conducted via the FastStart Universal SYBR Green Master (ROX) (Roche) in the Applied Biosystems QuantStudio TM Real-PCR System (Q6). Human BGN primers were utilized, and the relative mRNA expression was determined using the comparative Ct method with Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as the reference gene. The primer sequences were as follows:

BGN-forward: 5′-TGA​CTG​GCA​TCC​CCA​AAG​AC-3′

BGN-reverse: 5′-GAG​TAG​CGA​AGC​AGG​TCC​TC-3′

GAPDH-forward: 5′-GTC​AGC​CGC​ATC​TTC​TTT-3′

GAPDH-reverse: 5′-CGC​CCA​ATA​CGA​CCA​AAT-3′

From January 2018 to September 2020, the tumors and adjacent tissues of 80 consecutive patients with GC after surgery in Suqian First People’s Hospital were collected. Patients who had received radiation or chemotherapy prior to surgery and had other malignancies were excluded from the study. After dewaxing, hydration, and thermal repair, the primary antibody against BGN (ab209234, Abcam, 1:2000) was incubated overnight at 4°C followed by incubation with detection polymer for 40 min at room temperature. 3,3′-Diaminobenzidine DAB (P0202, Beyotime Biotechnology co.) was used for signal detection. The images taken under the microscope were analyzed using the IHC profiler plugin of ImageJ software (Varghese et al., 2014). Finally, SPSS version 23.0 software was used to statistic the results.

According to the mean value of BGN expression, the data from the TCGA cohort were divided into high expression group and low expression group, and the DESeq2 package (Love et al., 2014) was used for differential analysis. DEGs were defined as having a p.adj <0.05 and |logFC| >1.5. The details of the DEGs were visualized using the volcano map.

After ID conversion of identified DEGs via or.Hs.eg.db package, further functional enrichment analysis was performed through clusterProfiler package (Yu et al., 2012). Enrichments that satisfied the following conditions were considered significant: p.adj<0.05, and q-value<0.2. DEGs results were employed for gene-set enrichment analyses (GSEA) and building gene-set enrichment plots against the Molecular Signatures Database (MSigDB) hallmark gene sets through the R package, clusterProfiler, and significance was set as an adjusted p < 0.05 and FDR<0.25.

After converting the level 3 HTSe1-FPKM format RNAseq data from the stomach adenocarcinoma (STAD) project of TCGA to TPM format, log2 conversion was performed. After normal tissue samples were removed, data from a total of 375 STAD samples were retained for subsequent analysis. The relative tumor infiltration levels of immune cell types were quantified using ssGSEA of clusterProfilerpackage (Yu et al., 2012) to quantify the relative tumor infiltration levels of immune cell types, and the marker genes of immune cell types for single-sample gene-set enrichment analysis (ssGSEA) were obtained from published signature gene lists (Bindea et al., 2013). Spearman’s Correlation Test was adopted to determine a correlation between BGN and the immune infiltration levels and the association of immune infiltration with the different expression groups of BGN.

For TCGA data, Wilcoxon signed Rank-Sum test and logistic regression analyses were used to evaluate the relationship between BGN expression and clinicopathological variables. Moreover, univariate and multivariate Cox regression analyses were used to compare the effects of BGN expression and other clinicopathological variables on the overall survival of GC patients. Multivariate Cox regression analysis was used to examine the independent factors affecting the prognosis of GC.

Furthermore, we collected clinicopathological data from 80 patients who underwent immunohistochemistry to evaluate the relationship between BGN expression and clinicopathological variables. Chi-square tests were used to evaluate the relationship between gender, pathological type, residual tumor status, and BGN expression. Fisher’s exact tests were used to evaluate the relationship between pathologic stage, T stage, N stage, primary treatment outcome, and BGN expression. Wilcoxon signed Rank-Sum test was used to evaluate the relationship between age and BGN expression.

The identified independent factors associated with GC prognosis were used to construct a nomogram that predicted the probability of 1-, 3-, and 5-years survival in patients with GC. The prognostic data were obtained from a study by Jianfang Liu(Liu et al., 2018). Nomogram was constructed by R package with the survival and rms package. The Harrell’s concordance index (C-index) was used to quantify the predictive accuracy, which ranges from 0.5 (no predictive power) to 1 (perfect prediction). Furthermore, calibration plots were generated to examine the performance characteristics of the predictive nomogram.

Significant differential expression of BGN was documented in most of the 33 cancers, including in STAD (Figure 1A). The expression of BGN in GC (375 cases from TCGA) was significantly higher than in normal tissues (32 para-cancer tissues from TCGA and 174 normal tissues from GTEx) (p < 0.001) (Figure 1B). Similarly, the comparison of 27 tumor tissues in TCGA with the corresponding para-cancer tissues also showed significant expression of BGN in tumor tissues (Figure 1C).

Furthermore, based on the expression profile of TCGA in tumor and normal tissues, a ROC curve of BGN for the diagnosis of GC was plotted. Figure 1D shows that in the prediction of tumor and normal outcomes, the variable BGN showed high accuracy (AUC = 0.945, CI = 0.915–0.975).

We further verified the BGN expression level using RT-PCR (Figure 2E, p = 0.0068) and IHC (Figures 2A–D). The results were consistent with those in the TCGA database, indicating significantly higher levels of BGN expression in GC than that in normal tissues.

The volcano map shows the expression of identified DEGs between groups with high and low BGN expression (Figure 3A). of all the 492 DEGs. Of them, 207 were up-regulated, and 285 were down-regulated genes.

In terms of Biological Process (BP), most of the DEGs were enriched in extracellular structure organization, extracellular matrix (ECM) organization, and skin development. In terms of cellular components (CC), DEGs were mostly enriched in the collagen-containing ECM, endoplasmic reticulum lumen, and ECM components. In terms of molecular functions (MF), the DEGs also showed significant association with ECM structural constituent, receptor-ligand activity, and glycosaminoglycan binding. Furthermore, they were found mainly enriched in three KEGG pathways, including protein digestion and absorption, ECM-receptor interaction, focal adhesion pathway (Figure 3B; Table 3). GSEA analysis revealed the following BGN-related enrichment pathways: collagen formulation, immunoregulatory interactions between a lymphoid and a non-lymphoid cell, focal adhesion, ECM glycoproteins, Wnt signaling, and signaling by vascular endothelial growth factor (VEGF), as shown in Figures 3C–H.

The BGN expression showed positive correlation with the enrichment of the NK cells (r = 0.620, p < 0.001) and macrophages (r = 0.550, p < 0.001) but negative correlation with the enrichment of Th17 cells (r = -0.250, p < 0.001) (Figures 4A–G; Table 4). The enrichment score of macrophages (High: 0.501 ± 0.061, Low: 0.44 ± 0.066, p < 0.001) and NK cells (High: 0.47 ± 0.031, Low: 0.433 ± 0.036, p < 0.001) in the group with high BGN expression was significantly higher than that in the group with low BGN expression, while the enrichment score of Th17 cells (High: 0.218 ± 0.111, Low: 0.266 ± 0.12, p < 0.001) in the group with high BGN expression was significantly lower than that in the group with low BGN expression (Table 5). The details of immune cell enrichment score in the BGN high expression group and low expression group are shown in Table 5.

BGN expression was remarkably correlated with histologic grade (Figure 5A, GI&G2 vs. G3, p < 0.001), histologic type (Figure 5B, Diffuse type vs. Tubular type, p < 0.001), histologic stage (Figure 5C, stage I vs. stage II and stage I vs. stage III, p < 0.001), T stage (Figure 5D, T1 vs. T2, T1 vs. T3, and T1 vs. T4, p < 0.001) and Helicobacter pylori (HP) infection (Figure 5G, yes vs. no, p < 0.05) in gastric cancer (GC). However, the following clinicopathological features showed no significant association with BGN expression: M stage, N stage, residual tumor, gender, age, primary therapy outcome, and Barrett’s esophagus (Figures 5F,H–L, p > 0.05).

The results in Table 6 showed that BGN expression was remarkably correlated with pathologic stage (p = 0.008), T stage (p = 0.001), histologic type (p < 0.001), and histological grade (p = 0.025) in 80 GC patients who underwent immunohistochemistry, but was not significantly associated with gender (p = 0.802), N stage (p = 0.232), residual tumor (p = 0.323), primary therapy outcome (p = 0.655), anatomic neoplasm subdivision (p = 0.905), and age (p = 0.600).

The results of survival analysis revealed significant association of greater BGN expression with poor Overall Survival (OS) in GC patients (Figure 6A, HR = 1.53 (1.09–2.14), p = 0.013), but no significantly association with Disease Specific Survival (DSS) (Figure 6B, HR = 1.43 (0.94–2.19), p = 0.095), and Progress Free Interval (PFI) (Figure 6C, HR = 1.27 (0.89–1.81), p = 0.189).

In order to eliminate the influence of other clinicopathological variables on OS of GC, multivariate Cox regression analysis was performed to identify independent factors affecting OS of GC. Table 7 and Figure 6D show that pathologic stage (stage I &II vs. stage III &IV, HR (95% CI) = 1.604 (1.022–2.517), p = 0.040), primary therapy outcome (CR vs. PD &SD &PR, HR (95% CI) = 4.594 (2.938–7.182), p < 0.001), age (≤65 vs. >65 years, HR (95% CI) = 1.654 (1.089–2.514), p = 0.018), histologic grade (G1 & G2 vs. G3, HR (95% CI) = 1.576 (1.014–2.451), p = 0.043), and BGN (low vs. high, HR (95% CI) = 1.798 (1.183–2.732), p = 0.006) had significant correlation with OS rates in patients with GC. However, BGN expression showed no association with poor DSS and DSS PFI (Tables 8; Tables 9).

A nomogram to predict 1-, 3-, and 5-years’ OS probability was constructed on the basis of multivariate Cox regression analysis. In it, five variables, namely pathologic stage, primary therapy outcome, age, histologic grade, and BGN expression level, were used. Figure 7A depicts 11 rows in the nomogram, with the rows ranging from 2 to 6 representing the above variables. The points of the five variables were added up to the total points, which were displayed in row 7 and corresponded to the linear predictor in the prediction of 1-, 3-, and 5-years survival probability in row 8. The C-index was used to quantify the predictive accuracy, ranging from 0.5 (no predictive power) to 1 (perfect prediction). The C-index of this nomogram was 0.728 (0.705–0.752), indicating that the prediction was in good agreement with the actual survival probability. The nomogram calibration plot (Figure 7B) also suggests that the nomogram was well-calibrated, with the mean predicted probabilities close to observed probabilities.