3.1 Identification and functional analysis of DEGs

We analyzed and screened the DEGs based on different stromal scores and immune scores in low immune score samples and high immune score samples, low stromal score samples, high stromal score samples, respectively, and further plotted the heatmap of the top 20 DEGs with low and high expression to determine the gene expression in different components (Fig. 1a, b). We screened 2315 genes in total. A total of 1380 DEGs, including 268 down-regulated genes and 1112 up-regulated genes, were screened according to stromal scores. Similarly, there were 935 DEGs (169 down-regulated genes and 766 up-regulated genes) between high immune scores and low immune scores samples. A total of 656 genes were obtained, which are genes related to the TME (Fig. 1c, d). According to GO analysis (Fig. 2a, b), DEGs are mainly associated with a range of immune responses, such as leukocyte migration and T cell activation. Molecular function shows that DEGs play a potential immunological role in TME by extracellular matrix structural constituent, glycosaminoglycan binding, heparin binding and immune receptor activity. Collagen-containing extracellular matrix, the external side of plasma membrane and the MHC protein complex enhanced in the cellular component. As shown in KEGG analysis, most DEGs were involved in the intestinal immune network for Th17-cell differentiation, cytokine-cytokine receptor interaction, Th17-cell differentiation, Th1 and Th2 cell differentiation and IgA production (Fig. 2c, d). This is a further indication that there is some regulation and control of immune regulation of TME by DEGs.

Fig. 1

Differential genes in GC patients (a). Heatmap of differential genes between high and low scoring groups in immune cells (b). Heat map of differential genes between high and low scoring groups in stromal cells (c). Venn diagrams of DEGs in low-expressing stromal and immune components (d). Venn diagrams of DEG in high expressing stromal and immune components

Fig. 2

GO and KEGG enrichment analysis of DEGs (a). b 656 DEGs for molecular functions, cellular components, and biological processes of the top six most significant functions (c). d 656 DEGs for KEGG enrichment analysis

3.2 The relationship between scoring scores with clinicopathological features and survival

The corresponding clinical characteristics of 433 patients were analyzed in order to examine the association between the clinical characteristics of patients with GC and Estimate scores, Immune scores, and Stromal scores. The T-stage of the tumor was found to significantly differ from Stromal scores and Estimate scores (p = 0.0013, p = 0.0069, Fig. 3a, b). The Kaplan–Meier method was used to examine the relationship between each score and patient survival in more details. Stromal scores were significantly different and negatively correlated with patient survival [p = 1.8e−3, HR = 1.56, 95%CI (1.18, 2.08), Fig. 3c], while Immune scores were significantly positively correlated with patients’ OS [p = 6.4e−3, HR = 0.64, 95%CI (0.46, 0.88), Fig. 3d], suggesting that the lower the stromal cells content and the higher the immune cells content in tumor patients, the better prognosis. The Estimate scores and patient survival showed no difference (p = 0.22, Fig. S1). This indicated that the survival prognosis of patients was correlated with the stromal and immune components of TME.

Fig. 3

Stromal, immune and composite scores of DEG about GC. (a) Relationship between (b) clinicopathological features and (c-d) prognosis

3.3 PPI network and univariate COX regression analysis

To further analyze the interactions between DEGs in GC, we constructed a PPI network associated with 656 differential genes using the STRING database, the confidence of the minimum required interaction was 0.900. In the PPI network, a total of 627 nodes, 536 edges, 1.71 average node degree, 0.278 average local clustering coefficient and < 1.0e−16 PPI enrichment p value were generated. Then mapped the results of interactions between 263 associated proteins using Cytoscape software (Fig. 4a), and collated the top 30 differential genes in the PPI network in terms of the number of protein-adjacent nodes (Fig. 4b). Univariate COX regression analysis showed that 163 genes were associated with survival risk factors for GC (Fig. 4c). A cross-tabulation analysis of the PPI network results in combination with the results of the univariate COX regression analysis confirmed CTSK as the only overlapping gene (Fig. 4d).

Fig. 4

PPI network with univariate COX regression analysis (a). PPI network constructed by 656 DEGs, green indicates proteins with up-regulated expression and red indicates proteins with down-regulated expression (b). Top 30 genes by the number of adjacent nodes of the PPI network (c). Top 50 genes associated with GC prognosis by univariate COX regression analysis (d). Cross-tabulation analysis of the PPI network with univariate COX regression analysis yielded CTSK as the only key gene

3.4 Pan-cancer expression landscape of CTSK

Based on the analysis in the TCGA GTEx pan-cancer dataset (N = 19,131, G = 60,499), it could be observed that CTSK gene was significantly higher-expressed in 17 cancers such as LUAD, STES, STAD, KIRC, LUSC, LIHC, PAAD, CHOL than their adjacent normal tissues, and that another 10 tumors were observed to be significantly down-regulation, such as UCEC, BRCA, ACC, PRAD, THCA, KICH.CTSK expression was observed to be elevated in tumor tissues in 23 cancer species (Fig. 5). It can be speculated that CTSK may be a potential high-risk gene for cancer, and that its high expression may lead to cancer development.

Fig. 5

Gene expression of CTSK in 34 cancers Up-regulation of CTSK gene expression: GBM (p = 6.0e−45), GBMLGG (p = 1.8e−65), LGG (p = 3.6e−35), LUAD (p = 7.3e−8), STES (p = 1.1e−11) KIRP (p = 5.3e−3), STAD (p = 2.8e−20) HNSC (p = 2.9e−7), KIRC (p = 8.1e−3), LUSC (p = 4.2e−15), LIHC (p = 8.7e−6), PAAD (p = 1.2e−101), TGCT (p = 3.9e−4), ALL (p = 1.4e−6), LAML (p = 3.7e−42), PCPG (p = 0.04), CHOL (p = 2.4e−8), expression of downward adjustment: UCEC (p = 7.5e−11), BRCA (p = 1.8e−4), CESC (p = 1.5e−6), PRAD (p = 1.9e−15), BLCA (p = 5.6e−6), THCA (p = 7.7e−46), OV (p = 5.2e−60), UCS (p = 2.6e−15), ACC (p = 6.1e−12), KICH (p = 1.2e−7).*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

3.5 Relationship between CTSK expression and clinicopathological features and cancer progression

Based on the downloaded clinical data from GC patients, the correlation between patient gender, age, T, N clinicopathological characteristics and survival prognosis with CTSK expression was analyzed. The analysis showed no significant correlation between CTSK expression with gender, age and T-stage of patients (Fig. S2), but there was a significant association with N-stage of lymph node metastasis (p = 0.0018, Fig. 6a). Survival analysis shows that patients with high CTSK expression have 1.73-fold higher risk of death than those with low CTSK expression [p = 7.0e−5, HR = 1.73, 95%CI (1.32, 2.28), Fig. 6b]. In addition, high CTSK expression leads to a significantly higher infiltration of stromal cells and immune cells in TME (Fig. 6c). This suggests that CTSK is closely associated with the progression of GC, especially early lymph node metastasis and prognosis. In a further analysis, a survival analysis of the validation set GSE26253 [p = 8.4e−3, HR = 1.50, 95%CI (1.11, 2.03), Fig. 6d] and using the GEPIA2 (p = 0.011, HR = 1.8, Fig. 6e) and K-M plotter [p = 0.0075, HR = 1.74, 95%CI (1.15, 2.63), Fig. 6f] platform revealed that OS as well as recurrence free survival (RFS) was significantly lower for patients in the CTSK high expression group than in the CTSK low expression group.

Fig. 6

Differences in CTSK gene expression in clinical versus TME (a). N-stage (b). High CTSK expression leads to lower overall survival (c). Differences between high and low expression of CTSK genes in TME between stromal score, immune score and composite score. d Validating the RFS relationship of CTSK in dataset GSE26253. e Expression levels of the gene CTSK in GEPIA in relation to prognosis. f OS analysis of CTSK in GC patients mapped by K-M Plotter

3.6 The relationship between CTSK and GC TME and tumor immune response

CTSK genes were found to be involved in the regulation of numerous GC immunological responses and signaling pathways by GSEA enrichment analysis. Among them, the chemokine signaling pathway, WNT signaling pathway, TGF-βsignaling pathway, etc. were significantly correlated with high CTSK expression, while RNA degradation was negatively correlated with CTSK expression (Fig. 7a). After analyzing the infiltration of 22 TIICs in GC tumor tissue and the correlation between immune cells, we found that the interaction between T lymphocytes and Tumor-Associated Macrophages (TAM) was more prominent (Fig. 7b, c). When the distribution of TIIC in GC tumors samples with low or high expression of CTSK was further analyzed, it was discovered a considerably larger proportion of TAM in the high CTSK expression group than the low expression group (Fig. 7d), and we found that CTSK expression was closely associated with nine immune cell types (Fig. 7e). Of these, macrophage M2 cells, mast resting cells, monocytes, CD4 memory resting T cells, B memory cells had positive correlations, and four were negatively correlated, including CD8 T cells, follicular helper T cells, memory activated CD4 T cells, and macrophage M1 cells. The correlation between 47 common immune checkpoint genes and CTSK expression are shown in (Fig. 7f). From the above results, it can be concluded that CTSK played an important role in the TME and immune response of GC that may become a new immunotherapeutic target for tumor treatment.

Fig. 7

Potential pathways of CTSK and immune response relationships. a GSEA enrichment results of CTSK. b Distribution of 22 types of immune cell infiltration in each GC patient. c Interrelationship between 22 types of immune cells. d box plot of the proportional difference in the distribution between TIIC and the two groups with high and low CTSK expression. e Correlation between CTSK expression and 22 TIICs. f Correlation between CTSK and immune checkpoint genes. (*p < 0.05, **p < 0.01, ***p < 0.001)

3.7 Analysis of CTSK expression levels and drug sensitivity

To determine whether the CTSK gene could be used as a biomarker for responding to the efficacy of drug therapy in GC patients, the IC50 values of various antitumor drugs were estimated. From the results, we observed that patients with lower levels of CTSK expression were more sensitive to treatment with Bosutinib, Methotrexate, Epothilone, Paclitaxel, 5-Fluorouracil, Obatoclax Mesylate and cell cycle-specific drugs Vinorelbine that effectively block cell division, while patients with high levels of CTSK expression responded more positively to Cisplatin, Cytarabine, pazopanib and multi-target inhibitors such as Linifanib, Cabozantinib (XL-184), Dasatinib, Saractinib, Ponatinib (AP-24534). Overall, there was a correlation between CTSK expression and sensitivity to drug treatment in GC (Fig. 8).

Fig. 8

Analysis of the sensitivity of high versus low CTSK expression groups to different drugs