Identification and validation of platinum resistance signature in gastric cancer

Abstract

Background

Platinum was the first drug with proven activity against gastric cancer (GC), the combination with fluoropyrimidine is the standard first-line systemic therapy for patients of GC. However, a major cause of treatment failure still is the existence of drug resistance. The purpose of this study is to identify and validate the platinum-related genes in GC and to construct a multi-gene joint signature for predicting the prognosis of GC patients.

Methods

Based on 326 platinum-related genes from GeneCards, GO and KEGG analysis were applied for differentially expressed genes in GC, UniCox regression analysis was used to select effective genes and Lasso-Cox regression was utilized to construct a prognosis model. Stratified analysis, CNV landscape, TMB and MSI status, HLA gene expression, GSEA and GSVA analysis, immune activities, immunotherapy sensitivities were evaluated in the resistant high and low groups. Drug resistant cell lines, PDO and PDX models were used to validate this signature.

Results

GO analysis of 140 differentially expressed genes were involved in many processes and KEGG pathways were enriched in platinum resistance and cancer. UniCox regression analysis was screened out 21 genes and conducted a platinum resistance scoring model. Stratified analysis indicated that the drug resistance score had a good predictive value in subgroups divided by T-stage, age and race. CNV changes were more occurred in the score-high group, and most model genes were negatively correlated with TMB, MSI and HLA gene expression. The immune score in resistant group was significantly higher, within more mast cell, regulatory T cell and dendritic cell infiltrated in. In vitro and in vivo models showed that 21 platinum resistance genes had varying degrees of upregulation under CDDP chemotherapy pressure.

Conclusions

The 21 gene-signature for platinum was developed to predict response to platinum chemotherapy for GC patients. It is worthwhile to further evaluate the molecular biology and the clinical applications of this signature.

Introduction

Gastric cancer (GC), a prevalent and aggressive disease worldwide, is often diagnosed in locally advanced or metastatic stage and, therefore, has a poor prognosis [1]. Chemotherapy is essential for individual ineligible for surgical resection, experiencing postoperative recurrence or metastasis, or having residual tumors [2]. Currently, platinum-based chemotherapy remains the standard management for GC patients, but the benefit of the treatment is often hampered by rapid development of drug resistance [3]. Despite initial rates of response to platinum-based chemotherapy can be higher than 50%, nearly all patients develop chemotherapy resistance and the median survival period is just extended to 10–12 months [4].

Platinum has become a very relevant drug in the treatment of patients of GC and other tumors [5]. For better outcome, it is requisite to clear up the molecular mechanisms underlying the phenomena of drug resistance, as they are the main cause of treatment failure and tumor progression. The discovery of these mechanisms and related biomarkers will aid in identifying individuals unlikely to benefit from platinum-based treatments and facilitate the development of novel regimens to overcome this resistance.

Platinum resistance in GC manifests as a sequential progression within tumor cells and microenvironment. As cancers advance to metastatic stages, their capacity to develop drug resistance intensifies, primarily attributed to the heightened expression of oncogenes or reduced levels of onco-suppressor factors [6]. To explore platinum signature genes in GC, we analyzed the differentially expressed platinum related genes by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, and established a scoring model integrating platinum related genes of GC. Gene Set Enrichment Analysis (GSEA), Gene Set Variation Analysis (GSVA), immune activity analysis and therapeutic drugs prediction of platinum resistance score were analyzed. Drug resistant GC cells, Patient Derived Organoid (PDO) and Patient Derived Tumor Xenograft (PDX) models were employed for the validation of this scoring model. Schematic overview of the study was shown in Fig. 1. In summary, we conducted a comprehensive platinum-related gene signature to further elucidate the mechanisms underlying the drug resistance in GC.

Fig. 1

Flowchart of the study. 375 stomach adenocarcinoma (STAD) and 32 normal tissue samples were included in this study. Drug resistance model was developed and validated. Further, GSEA, GSVA, clinical prognosis and immune activity analysis of drug resistance score were analyzed. Resistant cells, PDO and PDX were then used to validate the efficacy of drug resistance model

Materials and methods

Data processing

From TCGA (The Cancer Genome Atlas, https://portal.gdc.cancer.gov/), transcriptome profiling raw data of 375 stomach adenocarcinoma (STAD) and 32 normal tissue samples were downloaded with ‘count’ datatype (Supplementary Table 1). The somatic mutation of TCGA-STAD with ‘masked somatic mutation’ data were download and visualized by R package maftools [7]. The ‘mastered copy number segment’ data of TCGA-STAD were download by R package TCGAbiolinks [8] and analyzed the gene copy number with GISTIC 2.0 through GenePattern (https://cloud.genepattern.org) [9]. In addition, the clinical data of TCGA-STAD patients, including age, overall survival (OS), progression free survival (PFS), disease specificity survival (DSS), disease free survival (DFS), follow-up time, stage, etc. were downloaded from UCSC Xena. Finally, 350 patients were used as the training set of prognosis model by matching gene expression and clinical data of STAD patients. GSE15459 dataset including 200 GC patients from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) was acquired as validation cohort.

Differential expression analysis of platinum-resistance related genes

A total of 326 platinum-resistance related genes were obtained by searching the keyword ‘platinum-resistance’ in GeneCards database [10]. In order to analyze the impacts of these genes on GC, R package limma [11] was used to screen significant differential genes between STAD and normal samples in the dataset. The threshold values were set as the absolute value of fold change (FC) > 1.5 and P < 0.05. The differentially expressed genes related to platinum-resistance were obtained by crossing the differentially expressed genes and platinum-resistance genes and displayed by volcano map.

Functional enrichment analysis

Gene Ontology (GO) functional annotation [12] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis [13] were implied for the significant differential expression genes related to drug resistance by R package clusterProfiler [14].

Cox regression analysis

Univariate Cox (UniCox) regression was firstly employed to screen candidate genes from platinum-resistance related genes, then, the least absolute shrinkage and selection operator (Lasso-Cox) regression was applied to narrow down the genes screened from UniCox regression by R package glmnet. The expression value of each normalized gene weighted by the penalty coefficient of the characteristic gene was used to establish a drug resistance scoring formula to predict the prognostic of STAD:

$\begin{aligned} &\:\text{r}\text{e}\text{s}\text{i}\text{s}\text{t}\text{a}\text{n}\text{c}\text{e}\:\text{s}\text{c}\text{o}\text{r}\text{e}\:=\\&\:\sum\:_{\text{i}}\text{C}\text{o}\text{e}\text{f}\text{f}\text{i}\text{c}\text{i}\text{e}\text{n}\text{t}\:\left({\text{g}\text{e}\text{n}\text{e}}_{\text{i}}\right)\text{*}\text{m}\text{R}\text{N}\text{A}\\&\:\text{E}\text{x}\text{p}\text{r}\text{e}\text{s}\text{s}\text{i}\text{o}\text{n}\:\left({\text{g}\text{e}\text{n}\text{e}}_{\text{i}}\right) \end{aligned}$

Patients were classified into high and low score groups by median.

Prognostic analysis of drug resistance score

For evaluating the prognostic value of drug resistance score, 350 samples of STAD with OS were divided into high and low score groups according to the best cut-off value, and Kaplan-Meier analysis was performed using R package survival and survminer.

For testifying the ability of drug resistance score combined with clinicopathological characteristics to predict the OS of GC patients, univariate and multivariate Cox analysis was applied. Subsequently, the clinical predictive nomogram was constructed by R package rms using the drug resistance score and the selected clinicopathological characteristics [15]. By comparing the predicted value of the nomogram with the observed actual survival rate, a calibration curve was generated to evaluate the performance of the nomogram. R package ggDCA [16] was operated to conduct the decision curve analysis (DCA) and draw clinical impact curve to assess whether model-based decision is beneficial to patient prognosis.

Clinical stratification analysis

To explore the predictive value of drug resistance score in the subgroups, the patients were stratified according to age, sex, race and stage. The box chart of patient scores in each subgroup was drawn using R package ggplot2 [17].

Copy number analysis

To analyze and study the changes of somatic cell copy number between high- and low-score groups, ‘masked copy number segment’ data were combined with GISTIC 2.0 to identify copy number variation (CNV), determine gene deletion or amplification and get GISTIC score by setting the confidence level of 0.9.

Prediction of treatment response

Wilcoxon test was employed to explore the correlation between the resistance score and potential immunotherapeutic markers, including (immune checkpoint blockade) ICB related genes [18, 19], IFN- γ pathway markers [20], m6A and m5C regulatory factors [21, 22]. The GC cell line drug dataset was obtained from GDSC database [23], and compare the sensitivity of patients with high and low scores to therapeutic drugs by R package oncoPredict [24].

Identification of drug resistance score related biological characteristics

To investigate the differences in biological processes between patients with high and low drug resistance scores, Gene Set Enrichment Analysis (GSEA) by R package clusterProfiler [14] was used to evaluate the changes as the reference gene sets were “c5. go. v7.4. entrez. gmt” and “c2. cp. kegg. v7.4. entrez. gmt” from the MSigDB [25]. Also, 50 marker pathways of the reference gene set “h. all. v7.4. symbols. gmt” were downloaded, the enrichment fraction of each sample in each hallmark was calculated. Kaplan Meier analysis was used to validate the carcinogenic marker pathways’ prognostic value. P < 0.05 is considered statistically significant.

Immune infiltration analysis

To explore the infiltration of TME (tumor microenvironment) cells in two score groups, 28 reported immune cell types were tested by single sample gene set enrichment analysis (ssGSEA) algorithm [26], and the immune cell composition was visualized through box graph. Wilcoxon test was used to calculated the difference proportion of immune cells. The heat map of immune cell infiltration level was drawn by the R package corrplot, the relationship between drug resistance score and immune cells was employed by Spearman correlation analysis, the prognostic value of these significant immune cells was performed by Kaplan Meier survival analysis, and the correlation between hub gene levels and immune scores was achieved by R package estimate [27].

Genome variation analysis

To explore drug resistance score related somatic mutations and tumor mutation burden (TMB), microsatellite instability (MSI), human leukocyte antigen (HLA) genes, R package maftools [7] and cBioPortalData [28] were used to describe the landscape in two groups. R package ggplot2 was used to draw mutation lollipop diagrams or bar graphs to describe the correlation between the model genes and mutation sites or TMB, MSI, HLA genes.

Clinical relevance analysis

R package ggplot2 and pROC [29] were implied to describe the relevance of model genes and various clinical variables to evaluate the effect of drug resistance score in the model on clinical classification of GC patients.

Cell culture

GC cell lines MKN45 and NCI-N87 were purchased from American Type Culture Collection (ATCC, VA, USA) and their cisplatin (CDDP, Calbiochem, CA, USA) resistance sublines NCI-N87/DR and MKN45/DR were established in our institute. All experimental cells were cultured in DMEM supplemented with 10% FBS, 100 U/mL penicillin at 37 °C in 5% CO2 incubator. Short tandem repeat (STR) analysis with a test date of October 2023 was applied to the identification of the cell lines used. Each cell line was tested for mycoplasma contamination before use.

IC50 determination

Cells with 4000/well were seeded into 96-well plates the day before, the medium was replaced by fresh medium with concentration gradient of 0, 5, 10, 25, 50 and 100 µM cisplatin on the day. After incubation at 37 °C for 48 h, cell viability was performed by CCK-8 assay. The half maximal inhibitory concentration (IC50) value was defined as the concentration that resulted in a 50% reduction in cell growth compared with control. Three independent experiments were performed.

RNA extraction and qRT-PCR

Total RNA was extracted with TRIzol (Invitrogen, CA, USA) and cDNA was synthesized by reverse transcription with 1 µg RNA using Superscript III reverse transcriptase (Invitrogen). Quantitative RT-PCR (qRT-PCR) was performed by SYBR Green PCR master mix (ThermoFisher, MA, USA) on the ABI 7900. GAPDH was served as an internal reference and 2^−ΔΔCt method was used to calculate relative mRNA expression. Supplementary Table 2 provides the primer sequences.

Patient derived organoid (PDO) culture

The tissues from GC patients were cut into 1 mm³ pieces and digested with 2 mM EDTA at 4 °C for 3 h with shaking. Glands were seeded in Matrigel (BD Biosciences, CA, USA) and overlaid with medium containing advanced DMEM/F12 (Life Technologies, CA, USA, the same below) supplemented with B27 (1: 50), N2 (1: 100), N-acetylcysteine (50 mM), recombinant epithelial growth factor (EGF, 50 ng/mL), Noggin (100 ng/ml), R-Spondin (1 mg/mL), Glutamax-I Supplement (1: 100), HEPES (10 mM) and Penicillin/Streptomycin (400 mg/mL). PDO were passaged once a week with a split ratio of 1:2/1:3 and used between passage 2 and 4. Treatment with CDDP was performed 24 h after seeding.

Patient derived xenograft (PDX) construction

Fresh GC tissues were cut into 1–2 mm³ pieces and inoculated subcutaneously in NOD scid gamma (NSG) mice under anesthesia with isoflurane. At about 1000 mm³, tumors were extracted for serial transplantation into NSG mice and used for the experiment when tumor volume reached about 100 mm³. PDX-bearing mice were randomly divided into groups to receive CDDP (4 mg/kg, twice a week) or PBS (control) via intraperitoneal injection (i.p.) for two weeks. Tumor volume (volume = length × width× width/2) was measured every four days. After experiment completion, mice were euthanized, followed by tumor dissection and processing for qRT-PCR analyses.

Statistical analysis

Data processing and analysis were performed using R software (v4.1.1). For comparing two groups of continuous variables, the independent Student’s t-test was used for normally distributed data, and the Mann-Whitney U test for non-normally distributed data. For categorical variables, the Chi-square test or Fisher’s exact test was applied. P < 0.05 was statistically significant.

Results

Development of platinum resistance related genes in GC

Total of 326 platinum resistance related genes from GeneCards were used to intersect with TCGA-STAD expression profile, then, 140 significant differential genes were obtained through differential gene analysis between tumor and normal samples, including 93 up-regulated and 47 down-regulated genes (Fig. 2A and B).

Fig. 2

Pattern diagram and biological process of cisplatin resistance related genes. (A) The volcano map of differentially expressed genes between TCGA-STAD tumor and normal samples. Red: up-regulated differentially expressed genes; Blue: down-regulated differentially expressed genes; Gray: the genes that are not significantly differentially expressed. (B) The heat map of differential gene expression TCGA-STAD tumor and normal samples. C-E. BP (biological processes, C), CC (cell components, D), MF (molecular functions, E) analysis in GO terms of differentially expressed genes. F. KEGG enrichment analysis of differentially expressed genes. Node color indicates the expression level of genes, quadrilateral color indicates Z-score corresponding to biological functions

For survey of the discrepancies in biological processes between tumor and normal samples, GO functions of 140 differentially expressed genes were annotated and the genes involved in many processes were traced. In terms of biological processes, these genes were mainly enriched in gland development, signal transmission by p53 class mediator, signal transmission in response to DNA damage, mitotic cell cycle phase transition and aging (P < 0.05, Fig. 2C); among the cell components, the genes were concentrated in chromosomal region, translation regulator complex, nuclear chromosome, condensed chromosome and lateral element (P < 0.05, Fig. 2D); for molecular functions, the genes were accumulated in platelet derived growth factor receiver binding, RNA polymerase II specific DNA binding, translation factor binding, DNA binding platelet derived growth factor binding and DNA binding translation activator activity (P < 0.05, Fig. 2E). In the meantime, these genes were enriched in KEGG pathways such as platinum drug resistance, pancreatic cancer, prostate cancer, melanoma, endocrine resistance (P < 0.05, Fig. 2F).

Construction and validation of cisplatin resistance scoring model

Subsequently, the differential expression genes related to cisplatin resistance were imported into UniCox regression analysis and 21 effective prognostic genes (EZH2, TMEM88, KIT, PDGFD, TNFAIP2, PDGFRA, COL3A1, SNAI1, SFRP4, RBMS3, E2F2, NTRK3, VCAN, DNMT1, GLI1, PDGFRB, FANCC, KLK4, NUAK1, NCALD, TMEM108) were selected with P < 0.05. Then, Lasso-Cox regression on 21 genes was utilized to obtain the penalty coefficient of each gene (Fig. 3A). The prognostic correlations of these 21 genes were also analyzed, and the results showed that most of them had significant positive correlations (r > 0.3, P < 0.05, Fig. 3B).

Fig. 3

Lasso-Cox regression analysis. (A) Screening of characteristic genes of drug resistance score. (B) Correlation analysis of expression levels of characteristic genes. Red indicates positive correlation and blue indicates negative. (C) Risk factor correlation diagram of drug resistance score. Upper: the drug resistance score arrangement diagram of patients, yellow represents high drug resistance score and blue represents low; Lower: correlation analysis between survival time, survival status of patients and drug resistance score, yellow represents dead and blue represents alive

Further, the risk factor correlation chart was drawn according to the drug resistance score of GC patients. The results showed that the number of deaths in the high score population is significantly more than that in the low group and the survival time of the high score population is shorter (Fig. 3C and Supplementary Fig. 1A). These results indicate that the drug resistance score is a potentially reliable predictor.

Clinical prognosis analysis of drug resistance score

Based on drug resistant score model, survival analysis was further performed. Compared with patients with low drug resistance score, the patients in high score group had a significantly worse prognosis (P < 0.05, Fig. 4A). At the same time, time ROC curves were drawn according to the survival information combined with drug resistance score (Fig. 4B). The area under the curve (AUC) values of 1-, 2- and 3-year survival were predicted to be 0.540, 0.575 and 0.536, respectively. Further, the stratification analysis indicated that the drug resistance score displayed a good predicted value in the subgroups divided by age, sex, race and stage (P < 0.05, Supplementary Fig. 1B).

Fig. 4

Clinical prognosis analysis of drug resistance score. (A) Survival curve of drug resistance score, the blue represents the patients with low score while the yellow represents high ones. (B) 1-, 2- and 3-year survival time ROC curve of drug resistance score in GC patients. (C) Univariate and multivariate analysis of drug resistance score and prognostic indicators. (D) Nomograms of prognosis in GC patients. (E) Nomogram model correction curve. (F) The survival prediction DCA curve, the line color represents the survival prediction of different characteristics

The dataset GSE15459 was used to verify the impact of drug resistance scores on the prognosis of patients with GC. The OS between patients with low and high scores displayed a significant difference in GSE15459 (P < 0.001, Supplementary Fig. 2A). Moreover, the AUC values predicted by the drug resistance score for 1-, 3- and 5-year survival were 0.619, 0.626 and 0.679 (Supplementary Fig. 2B). These results demonstrate that the resistance score has predictive value.

To figure out whether the drug resistance score prediction model is an independent prognostic indicator in GC, univariate and multivariate analysis were conducted and the result showed that drug resistance score and age were prognostic factors of GC (low risk score: HR < 1, P < 0.05; age > = 65: HR > 1, P < 0.05, Fig. 4C). The nomogram was constructed based on drug resistance score and age (Fig. 4D). After adjusting the nomogram model, we found that the predicted survival of the patients was very close to the actual survival of the patients (Fig. 4E). In order to assess the model accuracy and predict the benefits of patients’ intervention in line with the prediction results, we compared the drug resistance score, age, stage and decision curve analysis (DCA) of drug resistance score and clinical characteristics. The prediction line remained above the standard line, indicating that these factors could bring good benefits to patients. The decision based on the nomogram model was above all clinical characteristics which suggested the model was most beneficial to the prognosis of GC patients (Fig. 4F).

Clinical stratification analysis of drug resistance score

In addition, the patients were stratified according to age, sex, race and stage, then the distribution of patient scores in each clinical variables was observed by the box charts (Fig. 5A-G). Stratified analysis showed that the drug resistance score predicted well in subgroups divided by T-stage (P = 2.9e-05), age (P = 0.02) and race (P = 0.019), and also illustrated that the drug resistance score was a reliable prognostic feature.

Fig. 5

Clinical stratification and genome variation analysis. A-G. The box charts of drug resistance score with different classification, age, race and sex. H-I. Correlation of TMB (H) and MSI (I) with drug resistance score. J. The correlation analysis of HLA gene and model gene expression level. Node size represents significance and node color represents correlation

Then, we analyzed the landscape of CNV in drug resistance score high and low groups (Supplementary Fig. 2C-D). Interestingly, compared with the low score group, more CNV changes were occurred in the high, especially the copy number loss feature.

Further, the mutation lollipop diagrams of model genes were drawn to show the details of mutation sites (Supplementary Fig. 3A-S), and the survival analysis was conducted after the patients were grouped according to whether each gene was mutated (Supplementary Fig. 4A). In addition, the correlation analysis between model genes with TMB or MSI was conducted, and the results showed that most model genes had a significant negative relationship with TMB and MSI (P < 0.05, Fig. 5H and I). Further, the correlation between model gene and HLA gene expression was also analyzed and the results displayed that HLA gene was significantly negatively correlated with FNACC and SANI1 (P < 0.05), weekly positively correlated with GLI1, KIT, KLK4, RBMS3 and TMEM88 (P < 0.1, Fig. 5J). The weak correlation may be due to the presence of other important factors affecting the relationship.

GSEA and GSVA analysis of drug resistance score

Next, we conducted GSEA analysis on the patients in the low and high score groups, and the results showed that the biological processes such as adenylate cyclase binding, keratan sulfate catabolic process, anchored component of synaptic membrane, regulation of dopaminergic neuron differentiation, voltage gated sodium channel complex and other biological processes were promoted in the patients in the high rating group (Fig. 6A), as well as the cation channel activity and cation channel complex were displayed separately based on their significant P value and enrichment Score (Supplementary Fig. 4B-C); however, these biological processes such as kinetochore organization, double strand break repair via break induced replication, kinetochore assembly, DNA strand elongation involved in DNA replication, positive regulation of establishment of protein localization to telomer were inhibited in the high score group, kinetochore organization and chromium remolding at centromere were displayed individually (Fig. 6B and Supplementary Fig. 4D-E). Meanwhile, we also analyzed the activity of pathways in patients with low and high scores and the results indicated that the biological processes including renin angiotensin system, arrhythmogenic right ventricular cardiomyopathy, neuroactive ligand receptor interaction, dilated cardiomyopathy, hypertrophic cardiomyopathy were promoted in patients with high score, neuroactive ligand receptor interaction and renin angiotensin system were shown individually (Fig. 6C and Supplementary Fig. 4F-G); while homologous recombination, aminoacyl tRNA biosynthesis, mismatch repair, citrate cycle TCA cycle, RNA polymerase and other pathways were inhibited, as well as homologous recombination and aminoacyl tRNA biosynthesis were shown (Fig. 6D and Supplementary Fig. 4H-I).

Fig. 6

GSEA and GSVA analysis. A-B. GSEA-GO analysis: the promoted (A) and inhibited (B) biological processes in the samples with high drug resistance scores. C-D. GSEA-KEGG analysis: the promoted (C) and inhibited (D) biological processes in the high drug resistance score samples. E. The GSVA score of hallmarks between samples with high and low drug resistance scores, *P < 0.05; **P < 0.01; ***P < 0.001. F. The correlation between hallmarks and drug resistance score

To further survey the difference of cancer marker pathways between low and high score samples, GSVA analysis were carried on and great majority pathways were significantly enriched in these high score samples (P < 0.05, Fig. 6E), including hallmark UV response, hallmark TGF-beta signaling, hallmark myogenesis, etc., which are related to well-known carcinogenic pathways. The correlation between the drug resistance score and the GSVA score of cancer marker pathway of each patient revealed that most cancer marker pathways were positively correlated, especially with hallmark myogenesis, hallmark hedgehog signaling, hallmark epithelial mesenchymal transition, hallmark UV response and hallmark angiogenesis (P < 0.05, Fig. 6F). To evaluate the prognostic value of cancer marker pathways, Kaplan-Meier survival analysis was employed and different OS probabilities were observed between groups distinguished by these carcinogenic marker pathways, like, hallmark hypoxia, hallmark p53 pathway, hallmark PI3K-AKT-mTOR signaling and hallmark DNA repair (Supplementary Fig. 4J-M). To sum up, the drug resistance score is closely related to the carcinogenic pathways.

Immune activity analysis of drug resistance score

Tumor microenvironment (TME) is mainly composed of inflammatory cells, immune cells, interstitial cells, fibroblasts, various cytokines and chemokines, which is a comprehensive system. Immune cell infiltration plays an essential part in tumor research and prognosis prediction. We conducted ssGSEA to analyze the infiltration level of immune cells between low and high score samples, and the results showed that 19 types of immune cells, including activated CD4 + T cell, activated CD8 + T cell, and activated B cell were enriched in the high score group (P < 0.05, Fig. 7A). Moreover, drug resistance score showed a positive correlation with immune score (P < 0.05, Fig. 7B). Thirteen TME cell types were intersected by Venn diagram overlapped by differentia analysis, correlation analysis and survival analysis, as mast cell, effector memory CD4 + T cell, macrophage, plasmacytoid dendritic cell, immature dendritic cell, type 1 T helper cell, central memory CD4 + T cell, eosinophil, central memory CD8 + T cell, monocyte, gamma delta T cell, natural killer T cell, effector memory CD8 + T cell (both 19 cell types were screened from differentia analysis and survival analysis and fully coincided, Fig. 7C). These results suggest that the drug resistance score has a significant correlation with tumor immune infiltration and plays an important role in the stratification of GC patients. We also analyzed the immune score between low and high groups, and found that the immune score in drug resistance score high group was significantly higher than that in low group (P < 0.05, Fig. 7D).

Fig. 7

Immune activity analysis of drug resistance score. (A) The box diagram of TME cell infiltration scores in different drug resistance score samples. (B) The correlation analysis between drug resistance score and immune score. (C) Venn diagram of differential TME cells (yellow), TME cells related to drug resistance score (green) and TME cells affecting GC patient prognosis (blue). (D) The immune score of the samples with low and high drug resistance score. (E) The correlation analysis between the expression level of immune phenotype related genes and drug resistance score. F-G. TME cell content analysis in the high (F) and low (G) drug resistance score samples

Meanwhile, the correlation between drug resistance score and immune cell infiltration were also analyzed, and indicated that drug resistance score was positively correlated with most immune cells, which indicated that there were more immune cells in patients with high score (Fig. 7E). By correlation analysis, 23 TME cell types were screened and a positive correlation with drug resistance score was displayed (r > 0, P < 0.05, Fig. 7E). By calculating the correlation of immune cells in two groups, memory B cell, activated CD4 + T cell, gamma delta T cell, activated dendritic cell, natural killer T cell, macrophage, T follicular helper cell, regulatory T cell, MDSC, natural killer cell, central memory CD4 + T cell, activated CD8 + T cell, effector memory CD8 + T cell, type 1 T helper cell, monocyte, activated B cell and immune B cell showed a positive correlation in the high score samples (r > 0, Fig. 7F); gamma delta T cell, mast cell, macrophage, central memory CD8 + T cell, MDSC, regulatory T cell, T follicular helper cell, central memory CD4 + T cell, natural killer T cell, natural killer cell, monocyte, activated CD8 + T cell, effector memory CD8 + T cell, type 1 T helper cell, activated CD4 + T cell, activated dendritic cell, eosinophil, effector memory CD4 + T cell, activated B cell, immune B cell, memory B cell and type 2 T helper cell showed positive correlation in the low score samples (r > 0, Fig. 7G). Although the majority of immune cell types were similar, their degrees of correlation and contents of immune cells were dissimilarity. We also calculated the correlation between model gene and immune cells in two groups. The results showed that in the high score group, mast cell, plasmacytoid dendritic cell, T follicular helper cell, regulatory T cell and type 1 T helper cell were significantly correlated with the expression levels of characteristic genes (P < 0.05, Supplementary Fig. 5A), while the correlation was not prominent in the low score group (Supplementary Fig. 5B).

Clinical correlation analysis of model genes

Through the Human Protein Atlas (HPA), immunohistochemical information of 21 characteristic genes were checked and EZH2, TNFAIP2, PDGFRA, SNAI1, SFRP4, TMEM108 proteins were found to be highly expressed in GC tissues. Compared with normal tissues, the IHC staining of tumor tissues were much deeper, suggesting that these proteins were significantly overexpressed in GC (Fig. 8A). Although there is no significant elevation in the expression of COL3A1 [30], E2F2 [31], and GLI1 [32] in GC, some studies have shown their association with GC. We will conduct further in-depth analysis of these molecules in the future to provide stronger evidence to support our conclusions.

Fig. 8

Correlation between model genes and clinical variables. (A) IHC of characteristic genes including EZH2, SFRP4, E2F2, COL3A1, TMEM108, TNFAIP2, PDGFRA, SNAI1 and GLI1 in GC tumor and normal tissues from HPA database. (B) Gene expression of m6A regulators in two group samples. (C) The expression levels of m5C regulators in the two groups. (D) IFN-γ pathway markers in high and low drug resistance score samples. (E) ICB related gene expression levels in the two groups. *P < 0.05; **P < 0.01; ***P < 0.001

The correlation between model genes and various clinical variables were analyzed and the corresponding box graphs and ROC curves were draw (Supplementary Fig. 5C-I). These findings indicated that the model gene had good performance in distinguishing clinical variables such as stage and race of patients (P < 0.05). The results indicate that this model possesses a good predictive value in subgroups divided by stage and race, and that the drug resistance score is a reliable prognostic feature.

Drug resistance score May predict the immunotherapy sensitivity in GC patients

Considering the significant correlation between CD8 + T cells and drug resistance score and the pivotal role of m6A and m5C methylation in impairing the anti-tumor ability of CD8 + T cells, CD8 + T cell related IFN-γ pathway markers and m6A/m5C regulatory factors were analyzed subsequently in the different score groups. The results depicted that most of them were significantly related to the drug resistance score. Most m6A regulatory factors, such as CBLL1, CPSF6, EIF3A, EIF3B and EIF3H, were low expressed in the high score group (P < 0.05, Fig. 8B). Similarly, most m5C regulatory factors, such as ALYREF, DNMT1 and DNMT3B, were low expressed in the high score group (P < 0.05, Fig. 8C). For IFN-γ pathway, markers such as JAK2, IFNGR1 and SOCS3 were upregulated in the high score group (P < 0.05, Fig. 8D). It is reported that the levels of immune checkpoint blockade (ICB) associated genes are important to therapeutic response of immune checkpoint inhibitors (ICI). Targeting promising checkpoints has become a hopeful strategy for cancer treatment. To evaluate the potential of drug resistance score to predict the response of GC patients to immunotherapy, 50 immunomodulators in the two groups were determined. As shown in Fig. 8E, more than half were significantly related to drug resistance score, including CD200, CD276, CD28, CD40, etc., which were highly expressed in the high score group (P < 0.05). These findings indicate that the drug resistance score holds a great potential in immunotherapy response evaluation for GC patients.

Subsequently, in order to explore the difference in drug sensitivity between the patients in the high and low score groups, the IC50 values of the two groups of patients for different therapeutic drugs were predicted and the results showed that Temsirolimus_1016, QL-VIII-58_1166, Sepantronium bromide_268 were more sensitive in patients with high score (P < 0.05, Supplementary Fig. 6).

Identification of platinum resistance genes in vitro and in vivo

GC cell lines N87 and MKN45 were used to induce CDDP resistance in vitro, and the IC50 of parental and drug-resistant strains were calculated. The IC50 for N87/DR cells was 36.92 ± 1.01 µM, while for N87 cells was 19.10 ± 1.15 µM. The IC50 of MKN45/DR and MKN45 cells were 19.10 ± 0.91 µM and 9.49 ± 0.95 µM, respectively (Fig. 9A). QRT-PCR was used to detect the mRNA expression levels of 21 platinum resistance genes. Compared to the parent strain, the resistant strain showed significant upregulation of 11 genes in the N87 cell line and 7 genes in MKN45 cell line (Fig. 9B). Not all 21 genes were validated in our own induced drug-resistant GC cell lines, and we analyzed possible reasons for this result. First, the emergence of drug resistance in tumor cells is a complex and multifactorial phenomenon involving the interaction of many genes and the dynamic influence of regulatory networks. Second, because these 21 genes were obtained by bioinformatics analysis, and the background information of different cell lines was quite different, this was also the experimental phenomenon that caused the low expression of some predicted genes. In future studies, we will further analyze the genes that are not significantly expressed to enrich our research.

Fig. 9

Identification of platinum resistance genes in vitro and in vivo. (A) IC50 detection of N87 and MKN45 cells, and related drug-resistant strains. (B) The mRNA expression levels of 21 platinum resistance genes by real time PCR in parental and resistant cells. (C) Representative images of two GC organoids treated with CDDP in Day 0 and Day 2. (D) 21 platinum resistance gene levels in PDO-1T and 2T treated with CDDP for 2 days, compared with non-treated PDOs. (E) Photographs of tumors derived from nude mice injected with PDX-4918 treated with or without CDDP, and tumor growth kinetics in the two groups. (F) The expression levels of platinum resistance genes in PDX-4918 model. (G) Heatmap of 21 platinum resistance gene levels in drug resistant cell lines, PDOs and PDX models. The data are the mean ± SD of three independent experiments. *P < 0.05; **P < 0.01; ***P < 0.001

Then, we prepared organoids from resected primary human GC tissues and added 10 µM CDDP during the sphere growth period. After two days, the total spheres of PDO-1T and − 2T with CDDP were significantly smaller than the control groups (Fig. 9C). QRT-PCR of 21 platinum resistance genes showed that the expression levels were increased with varying degrees after CDDP was added to two PDOs (Fig. 9D). Due to the large heterogeneity of each patient, some genes did not show significant differential expression in the experimental validation, but we will continue to study these genes in the future when we have the energy and enough time to better explore the potential mechanism of platinum resistance in GC.

Given the effect of 21 platinum related genes on CDDP in resistant cells and in PDO models, we further validated the outcome in PDX model. When the subcutaneous tumor volume reaches about 100 mm³ after transplantation, mice were randomly divided into 2 groups (3 mice for each group) and treated with CDDP or PBS via intraperitoneal injection at 4 mg/kg every 4 days. All mice were subsequently euthanized and xenografts were harvested. Tumor volume was recorded every 4 days after inoculation, and the tumor size and growth rate of the medication group were inhibited, showing significant differences compared to the control group (Fig. 9E). QRT-PCR results indicated that 19 of the 21 platinum resistance genes in CDDP group were significantly upregulated (Fig. 9F).

Drug resistant cell lines, PDO and PDX models all showed varying degrees of gene upregulation under CDDP chemotherapy pressure, playing an important role in GC chemotherapy resistance (Fig. 9G).

Discussion

Since cancer is a leading cause of mortality worldwide, chemotherapy has been widely used as one of the most important strategies in its treatment [33]. Platinum is a common chemotherapeutic agent in cancer suppression that its broad-spectrum function against various kinds of cancers including breast, bladder, ovarian, testicular, esophageal, lung, gastric and colorectal cancers have been investigated [34,35,36,37]. The platinum-containing drugs including cisplatin, oxaliplatin and carboplatin have been widely utilized.

The anticancer function of platinum drugs is based on its interaction with DNA and forming intra strand and inter strand DNA crosslinks [38, 39]. As well as causing DNA damage, platinum owns the capability of increasing oxidative damage to mediate mitochondrial dysfunction in accelerating apoptosis in tumor cells [40]. Nevertheless, intrinsic or acquired resistance to platinum-based combinations is still a major cause of the failure of clinical effectiveness of chemotherapy. Therefore, it is of great significance to straight out the causes underlying this phenomenon in order to overreach it and to uncover better ways of fighting cancer.

In this study, we analyzed the mRNA expression of 326 platinum resistance related genes in the TCGA-STAD dataset. UniCox regression analysis showed that 21 genes were associated with the survival of GC. Dataset GSE15459 was verified the stratification of drug resistance scores on GC patient prognosis. The correlation analysis between the expression of model genes with CNV, TMB, MSI or HLA gene expression were conducted. GSEA and GSVA analysis indicated that the drug resistance score is closely related to the carcinogenic pathways. Moreover, the immune score in resistant-high group was significantly higher than that of resistant-low group, and such as mast cell, plasmacytoid dendritic cell, T follicular helper cell and regulatory T cell were more infiltrated in. Based on the prediction of Temsirolimus_1016, QL-VIII-58_1166 and Sepantronium bromide_268 were more sensitive in patients with high score, the drug resistance score also has great potential in evaluating the response of immunotherapy to GC patients.

Multiple factors can influence treatment outcomes of platinum-based chemotherapy [3, 34, 41]. Based on our results, immune system showed a crucial role in the effectiveness. We analyzed the composition and activation state of different immune cell types within the tumor microenvironment between score-high and -low groups, and found more presence of specific immune cell subsets in the high group, such as cytotoxic T cells and natural killer cells. High levels of infiltrated immune cells indicated a more active immune response against cancer cells, which can suggest increased sensitivity to platinum-based chemotherapy.

Finally, the identification of platinum resistance signature was conducted at three levels, CDDP resistant cell lines, PDO model and PDX model. Varying degrees of upregulation of signature genes were observed under CDDP chemotherapy pressure, indicating that this signature may play an essential part in GC chemotherapy resistance. Not all of the 21 genes were validated in these models. This may be attributed to the complex, multifactorial nature of drug resistance emergence in tumor cells, which involves the interaction of numerous genes and the dynamic influence of regulatory networks. Consequently, the expression of specific genes in models varies due to differing microenvironmental factors. Furthermore, since these 21 genes were identified through bioinformatics analysis, significant differences in background information among various cell lines were observed, leading to experimental phenomena where some predicted genes did not exhibit significant changes.

In this study, the proposed model demonstrated predictive capabilities; however, certain limitations remain. The research was conducted with a relatively small sample size, which may affect the generalizability of the findings. Therefore, further validation in larger, prospective clinical cohorts is essential to confirm the model’s efficacy in real-world settings. Additionally, the heterogeneity of sample sources and patient characteristics may limit the model’s universality. Furthermore, the AUC values for survival prediction are relatively modest, indicating potential for improvement in predictive accuracy. Future research will focus on several enhancement strategies, including increasing the sample size and incorporating diverse patient populations to improve the model’s reliability and applicability, integrating additional potential biomarkers to enhance predictive capability, optimizing the model construction algorithm, and systematically comparing the model with existing models of platinum resistance.

Collectively, this platinum resistance signature could help to determine whether GC is resistant or sensitive to platinum drug, and may be a useful tool for clinical stratification. Nonetheless, integrating this biomarker into clinical practice may face several challenges. A primary obstacle is clinical acceptance, as both physicians and patients may initially exhibit skepticism towards the new model, necessitating time for adaptation and acceptance of this innovative approach to treatment guidance. Furthermore, data standardization emerges as a critical issue, given that disparate testing standards across hospitals and laboratories can affect the consistency of results. Therefore, establishing a unified standardized testing protocol is imperative to ensure data comparability across different institutions. By thoroughly addressing these challenges and implementing appropriate measures, we aim to facilitate the clinical translation of this drug resistance gene signature, thereby enabling more precise treatment plans for patients.

Conclusions

Summary, we developed the 21 gene-signature for platinum based on transcriptome, which could be applied to predict response to platinum chemotherapy for GC. The resistant samples classified by this signature showed multidimensional resistance related characteristics compared with the sensitive samples. It is worthwhile to further evaluate the clinical applications of this signature, which may assist clinicians to make a suitable strategy for GC patients.