Abstract Background Gastric cancer (GC) is characterized by significant intertumoral heterogeneity, which often leads to the development of resistance to platinum-based chemotherapy. Combining platinum drugs with other therapeutic strategies may improve treatment efficacy; however, the mechanisms underlying platinum resistance in GC remain unclear. Methods Key genes related to platinum resistance in GC were selected from the platinum resistance gene database and GC resistance datasets. The Similarity Network Fusion (SNF) algorithm was employed, along with prognosis-related methylation data and somatic mutation data, to classify the molecular subtypes of GC based on GC platinum resistance genes. Gene expression profiles, prognosis, immune cell infiltration, chemotherapy sensitivity, and immunotherapy responsiveness were comprehensively evaluated for each subtype. Localization and functional evaluation were conducted at the single-cell and spatial transcriptomics levels, and predictive models were developed using machine learning techniques. These functional differences in platinum resistance gene models were further explored in GC. Moreover, experimental validation was conducted to elucidate the mechanisms of key genes involved in platinum resistance in GC. Results Stomach adenocarcinoma (STAD) patients were classified into three subtypes using the SNF algorithm and multiomics data. Patients with subtype CS2 exhibited a significantly poorer prognosis than those with subtypes CS1 and CS3 (p < 0.05). Subtype CS1 was characterized as immune-deprived, CS2 as stroma-enriched, and CS3 as immune-enriched. Patients with subtype CS2 also exhibited the most adverse therapeutic responses to docetaxel, cisplatin, and gemcitabine. Single-cell analysis revealed high enrichment of M1 module cells with elevated expression of resistance genes, including the transcription factor KLF9. Spatial transcriptomic analysis further confirmed the independent spatial distribution of malignant cells with high expression of drug resistance genes (DRGs). Predictive models based on machine learning demonstrated excellent prognostic performance. Patients in the high DRG group also exhibited poorer responses to immunotherapy. Cellular experiments revealed that KLF9 overexpression significantly inhibited the proliferation of AGS cells (p < 0.05), reduced their resistance to platinum-based drugs, and markedly decreased the levels of inflammatory cytokines in them. Conclusion KLF9 was identified as a promising therapeutic target for overcoming platinum resistance in GC, warranting further investigation into its role and potential clinical applications. Electronic supplementary material The online version contains supplementary material available at 10.1186/s12967-025-06725-7. Keywords: Gastric cancer, Platinum resistance, Similarity network fusion, Spatial transcriptomics, KLF9 Background Gastric cancer (GC) is the fifth most prevalent malignancy worldwide. Due to lifestyle differences, 71.4% of GC cases and 70.1% of GC-related deaths occur in Asia [[40]7, [41]34]. Most patients with GC are diagnosed at advanced stages, with metastatic GC representing approximately 23% of cases in China and having a 5-year survival rate of only 9.4% [[42]37]. Despite advances in molecular-targeted therapies, chemotherapy remains the standard treatment for unresectable advanced GC [[43]49]. Platinum-based chemotherapy regimens, such as XLOEX and SOX, have significantly improved disease-free survival and overall survival (OS) in patients with GC. However, resistance to cisplatin (CP) often results in rapid disease progression [[44]18]. The mechanisms underlying CP resistance in GC are multifactorial, involving genetic and epigenetic modifications, altered signaling pathways, and disrupted cellular metabolism [[45]16, [46]45]. For example, overexpression of ABC transporters like P—glycoprotein enhances drug efflux, reducing intracellular platinum concentration. Enhanced DNA repair mechanisms, such as upregulated BRCA1, help tumor cells survive platinum-induced DNA damage [[47]16]. Tumor microenvironment factors, including ECM enrichment and hypoxia, also play roles. ECM—rich areas may impede drug delivery and promote resistance. Hypoxia alters gene expression, boosting tumor growth and metastasis. Studies have found potential biomarkers for platinum resistance, such as DUSP5P1 and specific lncRNAs like POU3F2 [[48]41]. These discoveries offer new ways to predict treatment response and develop targeted therapies to overcome platinum resistance. Molecular subtyping of cancers based on their molecular characteristics is essential for a deeper understanding of tumor heterogeneity and its impact on diagnosis, prognosis, and treatment [[49]10]. Traditional classification methods, which rely on the anatomical origin of tumors, often fail to capture molecular alterations that characterize disease progression. More accurate cancer subtypes can be identified by considering specific genetic and proteomic alterations within tumor cells. Molecular subtyping has revealed distinct subgroups with similar molecular features and biological behaviors, providing valuable insights into the mechanisms of cancer development and facilitating the formulation of more targeted and personalized treatment strategies. Molecular classification enables more precise prognostic predictions and facilitates tailored treatment approaches, such as targeted therapies that may be more effective for specific subtypes. Moreover, it offers critical insights into cancer research, assisting in elucidating the mechanisms underlying cancer initiation and progression. In this study, we integrated single-cell, transcriptomic, epigenomic, and somatic mutation data obtained from patients with GC and platinum resistance data sourced from the GC platinum resistance database. Using the Similarity Network Fusion (SNF) algorithm, we identified the molecular subtypes of GC and analyzed biological differences among subgroups to characterize key evolutionary patterns in GC progression. We also used spatial transcriptomics for spatial localization. In addition, we explored potential clinical treatment strategies based on specific subtype molecular characteristics, including chemotherapy and immunotherapy, paving the way for personalized medicine for patients with GC. Finally, we conducted a comprehensive functional validation of the key transcription factor KLF9 and explored its underlying mechanisms in GC (Fig. [50]1). Fig. 1. [51]Fig. 1 [52]Open in a new tab Flowchart of the present study Methods Data collection and processing Transcriptomic data for GC were sourced from the TCGA database ([53]https://portal.gdc.cancer.gov) along with corresponding clinical data. Methylation data, somatic mutation data, and copy number variation (CNV) data were downloaded from the UCSC Xena database ([54]https://xenabrowser.net/). The [55]GSE186205, [56]GSE62254, [57]GSE84437, and [58]GSE15459 datasets, containing patient survival data, were retrieved from the GEO database. High-throughput transcriptome sequencing data were converted into transcripts per million, followed by normalization and deduplication of expression profiles. The [59]GSE186205 dataset, derived from GC cell lines, specifically contained information related to platinum resistance. Differential gene expression analysis was performed using the “limma” package to identify differentially expressed genes associated with platinum resistance in GC. The gene set was intersected with the platinum resistance gene set obtained from the HGSOC-Platinum database ([60]http://ptrc-ddr.cptac-data-view.org/) to identify key platinum resistance genes in GC. The “ComBat” function within the sva package was used to integrate the [61]GSE62254 and [62]GSE84437 datasets into a metavalidation set for subgroup validation. The [63]GSE15459 and [64]GSE62254 datasets served as model validation sets to assess model performance. Identification of molecular subtypes related to platinum resistance The Multi-Omics integration and VIsualization in Cancer Subtyping (MOVICS) package [[65]23] was used to analyze gene matrices related to platinum resistance in conjunction with DNA methylation and somatic mutation data. Differential methylation analysis was performed on CpG sites, and univariate Cox regression analysis was conducted using TCGA survival data to identify CpG sites associated with OS. For somatic mutation data, the top 50 genes with the highest mutation rates were selected based on ranking mutation frequency. The results of the analyses of gene expression, DNA methylation, and somatic mutations were further investigated. The optimal number of clusters was determined after the initial round of feature selection. This was achieved using the “getClustNum” function within the MOVICS package, which integrates multiple metrics, including the Clustering Prediction Index (CPI), gap statistics, and silhouette scores, to robustly estimate the appropriate number of subgroups [[66]21]. In addition, the clustering structure was validated against previously established classifications of GC cells. Subsequently, clustering analysis was performed using the “getMOIC” function, incorporating 10 clustering algorithms: CIMLR, ConsensusClustering, SNF, iClusterBayes, PINSPlus, moCluster, NEMO, IntNMF, COCA, and LRA. These algorithms were specified as inputs to the “methodslist” parameter, maintaining the default parameters within the MOVICS package. The robustness of clustering was enhanced using the “getConsensusMOIC” function, which synthesizes the results of multiple algorithms through consensus clustering. During this process, the “distance” metric and “linkage” criterion were set to “Euclidean” and “average,” respectively. This approach ultimately yielded reliable clustering results. Biological characteristics and stability validation of platinum-resistant molecular subtypes The GC molecular subtypes were subjected to pathway enrichment analysis using the GSVA package [[67]13]. Enrichment data for transcription factors (TFs) were derived from a study by Lu et al. [[68]22], which included 23 TFs related to the induction/suppression of target genes and 71 candidate regulatory factors implicated in tumor chromatin remodeling. The distribution of immune checkpoints among the identified subtypes was systematically analyzed, and the immune and stromal scores of tumor tissues were estimated using the ESTIMATE package [[69]46]. Furthermore, the enrichment of 24 immune cell types within the tumor microenvironment [[70]43] was assessed using the GSVA package. To validate the stability of the identified subtypes, subtype-specific biomarkers were applied to a validation cohort to confirm clustering outcomes. Subsequently, whether the consensus clustering results were consistent with those obtained from the Nearest Template Prediction (NTP) and PAM classifiers was assessed. Subtype mutation analysis Somatic CNV data were visualized using the maftools package [[71]24], with annotations based on the hg19.mat reference genome. The GISTIC2.0 algorithm [[72]25] was applied to assess genomic alterations at the subtype level and identify significant differences in chromosomal deletions and amplifications among the subtypes. Tumor mutation burden (TMB) was calculated by measuring the number of nonsynonymous mutations per million base pairs. In addition, the MOVIC package’s built-in function was used to calculate the fraction of genome altered (FGA) [[73]29] due to copy number amplifications or deletions. Processing of single-cell and dropout data The single-cell dataset [74]GSE183904 included 40 samples, i.e., 11 samples of normal tissues and 29 samples of tumor tissues and peritoneal tumor tissues. Genes with expression levels greater than 300 were selected, and mitochondrial gene expression was filtered to below 15%. Doublet cells were removed using the Scrublet package [[75]38], and RNA contamination was eliminated using the DecontX package [[76]44], resulting in a total of 77,850 high-quality cells. After dimensionality reduction and clustering, batch correction was performed using the Harmony algorithm. Annotation and visualization were conducted using the CellMarker 2.0 database, based on the literature [[77]15]. In addition, malignant cells were identified using the InferCNV package. Spatial transcriptomics data were obtained from [78]GSE251950, specifically [79]GSM7990475 and [80]GSM7990482. The data were processed according to the recommended guidelines ([81]https://satijalab.org/seurat/articles/spatial_vignette.html) using the Seurat 4.0 package in R. Data normalization was performed using the SCTransform function, followed by dimensionality reduction and clustering via PCA and UMAP. Gene expression features were visualized using the SpatialFeaturePlot function. In addition, spatial transcriptomics data were analyzed using the Scanpy package in Python 3.10 [[82]42], and cell interactions and communications were explored via the Stlearn package. Acquisition of drug-resistant genes for each subtype and enrichment scoring Key differentially expressed genes were identified through differential analysis across transcriptomic subtypes. Five scoring methods were applied to evaluate single-cell data: AUCell, UCell [[83]1], singscore [[84]11], ssGSEA [[85]4], and AddModuleScore [[86]35]. After normalization, the scores were merged, resulting in a final score. Samples were categorized into high and low-expression groups based on the median value of the score before proceeding to subsequent analyses. Consensus non-negative matrix factorization (cNMF) cNMF is a tool used to infer gene expression patterns from single-cell RNA sequencing (scRNA-Seq) data. It takes a count matrix (N cells × G genes) as input and produces a gene expression program matrix of size (K × G) and a cell-program assignment matrix of size (N × K), indicating the program usage for each cell. Here, cNMF was applied to cluster malignant cells in GC and identify the cell modules most strongly associated with platinum resistance. The STARTRAC package [[87]48] was used to calculate the ratio of observed to expected cell numbers (known as the Ro/e index), which allowed comparison of differences in distribution among cell subpopulations. Cell–cell communication Cell–cell communication was analyzed using the R package CellChat (version 1.6.0) [[88]17] based on single-cell data and our previously established cell classification. The built-in reference database CellChatDB.human was used to investigate intercellular interactions based on 32 signaling pathways. In addition, the Nichenet package [[89]8] was used to gather protein–protein interaction data as well as data for the interactions of regulatory factors from multiple databases. This allowed the identification of downstream signaling targets of endothelial cells as ligands and the cross-validation of these findings via CellChat. SCENIC analysis The SCENIC package (version 3.1.4)(Van [[90]36]) was used to analyze the TF regulatory networks in high- and low-stemness groups. Enrichment analysis was conducted using the gene-motif rankings from the hg19-tss-centered-10 kb database in RcisTarget to investigate the regulatory networks between transcription start sites and genes in different endothelial cell populations. The results were visualized using the pheatmap package. Development and validation of prognostic models Based on the selected prognostic genes, we aimed to develop a consensus inverse regularized least squares model with high accuracy and stability by integrating 10 machine learning algorithms and 101 algorithm combinations. The combined algorithms included Random Survival Forest (RSF), Elastic Net (Enet), Lasso, Ridge, Stepwise Cox, CoxBoost, Cox Partial Least Squares Regression (plsRcox), Supervised Principal Components, Generalized Boosting Model, and Survival Support Vector Machine (survival-SVM). The workflow was as follows: (a) differential analysis of transcriptomic subtypes was conducted to identify key differentially expressed genes, followed by univariate Cox analysis using a bootstrap kernel to select key prognostic genes; (b) the key differentially expressed genes obtained were then subjected to 101 algorithm combinations to fit predictive models based on the leave-one-out cross-validation framework from the TCGA-LIHC cohort; (c) the models were validated based on an external validation dataset (HCCDB18); and (d) the Harrell’s Concordance Index (C-index) was calculated across all validation datasets for each model, and that with the highest average C-index was considered the optimal one. Using the optimal cut-off method to stratify high- and low-risk groups. Immune infiltration and immune therapy prediction in prognostic models The ESTIMATE algorithm was applied to expression data to assess stromal and immune cell scores in malignant tumor tissues [[91]46]. Differences in the distribution of immune cells between high- and low-expression prognostic groups were calculated using four methods within the IOBR [[92]47] package: CIBERSORT [[93]27], EPIC [[94]28], MCPCounter [[95]5], and TIMER [[96]20]. The IMvigor210 [[97]26] dataset, which included data from a Phase II clinical trial of atezolizumab for treating bladder cancer, was used to evaluate the immune therapy efficacy of the model. The SubMap algorithm [[98]14] was subsequently applied to predict the response to immune therapy based on the prognostic model, while the TIP algorithm was used to assess functional differences in T cells between groups. Screening of potential therapeutic agents for patients with high-expression DRGs Expression data for human cancer cell lines (CCLs) were obtained from the Broad Institute Cancer Cell Line Encyclopedia (CCLE) [[99]31]. Drug sensitivity data for CCLs were sourced from the CTRP v.2.0 ([100]https://portals.broadinstitute.org/ctrp and PRISM Repurposing datasets (19Q4; [101]https://depmap.org/portal/prism/). Then, the area under the curve (AUC) was calculated for both high- and low-expression groups to further assess drug sensitivity. Cell culture AGS cells, purchased from Beijing, China, were cultured in RPMI-1640 medium (10% FBS) at 37 °C with 5% CO[2]. qRT-PCR Total RNA was isolated from AGS cell lines utilizing an RNA extraction kit (Aibotek Biotechnology, Wuhan), followed by quantification of mRNA expression levels through primer-specific amplification. GAPDH was utilized as the internal reference. Primers design see Table [102]1. Table 1. The primer sequences employed in qRT-qPCR Primer name Sequence (5′– > 3′) IL1-β Forward AACCTCTTCGAGGCACAAGG Reverse TGGCTGCTTCAGACACTTGA IL-6 Forward TTCGGTCCAGTTGCCTTCTC Reverse TCACCAGGCAAGTCTCCTCA TNF-α Forward GCCCATGTTGTAGCAAACCC Reverse GGAGGTTGACCTTGGTCTGG KLF9 Forward ACAGTGGCTGTGGGAAAGTC Reverse GAAAGGGCCGTTCACCTGTA GAPDH Forward ATGTTGCAACCGGGAAGGAA Reverse GCAGGAGGCATTGCTGATGA [103]Open in a new tab ad-RNA transfection Gene silencing was performed using the HiPerFect Transfection Reagent (Qiagen, Germany). AGS cell lines were transfected with 5 nM KLF9-specific ad-RNA according to the manufacturer’s protocol, with non-targeting ad-RNA transfection and sham transfection serving as negative controls. Following a 24-h incubation period post-transfection, knockdown efficiency was quantitatively assessed through reverse transcription qRT-PCR analysis. Migration assays Migratory ability of AGS cells were detected by using the 24-well migration inserts (Corning, USA). For the migration assay, the upper wells were filled with 8 × 104 cells in 200 μL of RPMI-1640 (0.5% FBS) and 600 μL RPMI-1640 (10% FBS) was added to the lower chamber. Then we fixed the chambers with methanol after incubating it for 48 h. We then cleared up the upper chamber with the cotton swabs gently and stained it with 0.5% crystal violet. Cell scratch Seed AGS cells onto the bottom of a six-well plate and mark scratch lines to create wounds. Remove the old culture medium and divide the cells into two groups: si-NC, ad-NC (negative control) and ad-KLF9. Place the plate in a 37 °C, 5% CO2 cell culture incubator. At specified time points such as 0 h and 48 h, remove the cells from the plate and observe the width of the scratch at the same position under a microscope, taking pictures. Finally, use Image J software to analyze the distance of cell migration and the area of the scratch. Clone formation Transfection to overexpress KLF9, use high-concentration serum (10–20%) with ensured batch consistency. Gradient dilute cells to 100–500 cells/well, gently pipette to disperse them, and seed into 6-well plates with three replicates per group. Culture in a 37 °C, 5% CO[2] incubator for 10–14 days, performing half-medium replacement every 3–4 days to maintain stable culture conditions. Following incubation, wash twice with PBS, fix with 4% paraformaldehyde for 30 min, stain with 0.1% crystal violet for 20 min, rinse under running water, and air-dry. Quantify clones containing ≥ 50 cells by imaging both entire plates and individual wells. Statistical analysis All data processing and statistical analysis were conducted using R software (version 3.6.1) and GraphPad Prism. Student’s t-test and One-way analysis of variance (ANOVA) were used to determine differences between groups, and a p-value < 0.05 indicated statistical significance. Results Identification of three key subgroups associated with platinum resistance in GC via a multiomics approach Analysis of the [104]GSE18605 dataset resulted in the identification of a total of 3283 differentially expressed genes (p < 0.05). Following intersection with the known platinum resistance-associated genes (PRAGs), a total of 191 key PRAGs in GC were identified (Fig. [105]2A). By integrating the PRAG matrix with prognosis-related methylation sites and the top 50 somatic mutation sites, we constructed a matrix for 350 TCGA-STAD (GC) patients (Fig. [106]2B). The CPI and gap statistics were used to estimate the optimal number of clusters, and a k value of 3 was selected for subsequent analysis (Supplementary Fig. [107]1A-B). The integration of 10 different algorithms revealed that k = 3 significantly differentiated patients with STAD (Fig. [108]2C, D). Furthermore, survival analysis showed significant prognostic differences among the three subgroups (p < 0.001), with subgroups CS2 and CS1 exhibiting the most adverse prognosis and most favorable prognosis, respectively (Fig. [109]2E). Fig. 2. [110]Fig. 2 [111]Open in a new tab Multiomics-based identification of subgroups related to platinum resistance in GC. A Venn diagram of key genes related to platinum resistance in GC; B Heatmap of the matrix of platinum resistance–related genes showing differential expression, prognosis-related methylation sites, and top 50 somatic mutation sites; C Clustering of patients based on the integration of 10 multiomics clustering algorithms; D Consensus clustering matrix of the three prognostic subtypes based on the 10 algorithms; E Survival curves illustrating the differences among the three GC subtypes Biological functional differences among three GC subtypes Currently, GC subtypes are distinguished based on molecular expression levels and caseomics. Therefore, in this study, a new classification of GC s was attempted based on previous research. This classification considered the high heterogeneity of GC in the immune microenvironment, stromal microenvironment, genomic integrity, and oncogenic features. By analyzing 15 pathways, including immune, stromal, DNA damage repair, and oncogenic pathways, three distinct subtypes of GC, i.e., CS1, CS2, and CS3, were classified via GSVA [[112]21]. Subtype CS1, characterized as immunity-deprived (ImD), primarily features the upregulation of DNA damage repair and cell cycle pathways and the downregulation of immune, stromal, PI3K-Akt, Wnt, and TGF-β pathways. Subtype CS2, characterized as stroma-enriched (StE), primarily shows the upregulation of extracellular matrix (ECM), PI3K-Akt pathways, and other oncogenic pathways and the downregulation of Wnt, TGF-β, DNA damage repair, and cell cycle pathways. Subtype CS3, characterized as immunity-enriched (ImE), primarily shows the upregulation of immune, DNA damage repair, and cell cycle pathways. Prognostic analysis revealed that CS2 and CS3 had the poorest and most favorable prognoses, respectively, consistent with previous studies (Fig. [113]3A). Chromatin remodeling regulators related to cancer and 23 published TFs were analyzed to further study transcriptomic differences (Fig. [114]3B). Notably, the results showed that EGFR, AR, PGR, ESR1, STAT3, and FGFR1 were significantly activated in CS2, whereas GATA6, KLF4, FOXA1, and PPARG were specifically enriched in CS1. These differences in the activities of chromatin remodeling regulators further highlight the differential regulation among CSs, indicating that epigenetically driven transcription networks may be key distinguishing factors for these molecular subtypes. Considering the critical role of immune microenvironments in tumorigenesis and progression, we further evaluated differences in cellular infiltration among subtypes, and the results showed significantly increased stromal levels in CS2. Moreover, discernible differences were detected in the expression of immune checkpoints and immune cell subtypes (Fig. [115]3C), CS1 showed the lowest Immune Score (aligning with ImD status), CS2 showed the highest Stromal Score (aligning with StE status). Both CS2 and CS3 demonstrated elevated Immune Scores, with CS3 attaining maximal values. This stratification reveals distinct immune microenvironment profiles: while CS3 displays robust immune activation, CS2 shows concurrent stromal dominance with moderate immune infiltration. To substantiate the robustness of our molecular classification, the top 100 differentially expressed genes among the three CSs were analyzed within the META-STAD cohort (Fig. [116]3D), which confirmed the stability of our classification system. NTP was used to classify samples from external cohorts into the identified CSs, with CS2 in the META-STAD cohort exhibiting the worst prognosis among all subtypes (p < 0.001) (Fig. [117]3E). Simultaneously, the CS data were shown to be highly consistent with the results of the NTP and clustering algorithm (PAM) (p < 0.001) (Fig. [118]3F–H). Finally, GSEA analysis based on the differentially expressed genes in the three CSs indicated that CS2 exhibited high levels of oxidative stress, EMT, apoptosis, and KRAS pathway activation (F[119]ig. [120]3I). These findings not only deepen our understanding of the biological functional differences among GC subtypes but also hold promise for the development of targeted therapeutic strategies. Fig. 3. [121]Fig. 3 [122]Open in a new tab Biological functional differences among the three GC subtypes. A Subgroup enrichment analysis of GC-related biological pathways; B Regulatory activity profiles of 23 transcription factors (TFs) (top) and potential chromatin remodeling regulators associated with the three GC subtypes (bottom); C Immune characteristics in the TCGA-STAD cohort. The annotations at the top of the heatmap show immune enrichment scores for tumor-infiltrating lymphocytes and stromal enrichment scores, whereas those at the bottom indicate the enrichment levels for immune cells related to the tumor microenvironment; D Validation of the GC subtype classification in the META-STAD cohort; E Survival analysis of CSs in the META-STAD cohort; F, G Analysis of the consistency of CS data with NTP (F) and PAM (G) results in the TCGA-STAD cohort; H Analysis of the consistency of NTP results with PAM results in the META-STAD cohort; I GSEA enrichment analysis of the CS2 subgroup Drug prediction and mutation analysis of GC subtypes A comprehensive analysis of the drug sensitivity profiles of CS1, CS2, and CS3 to cisplatin, docetaxel, and gemcitabine was conducted to determine the therapeutic vulnerabilities of these three subtypes. The results showed that CS2 had the poorest sensitivity to these three drugs, consistent with previous findings, indicating that CS2 may be a drug-resistant subtype (Fig. [123]4A–C). Simultaneously, the genomic landscape indicated different TMBs among CSs, and the relative contributions of the four mutation signatures also varied. Intriguingly, the selected differentially expressed genes showed variation in mutation frequencies (Fig. [124]4D). Our analysis revealed variations in TMB among the CSs as well as divergence in the relative contributions of the four predominant mutation signatures. GSEA enrichment analysis suggested that CS2 is mainly enriched in pathways related to ECM construction, which may be one of the reasons for platinum resistance (Fig. [125]4E). Finally, GISTIC2.0 was used to calculate CNV frequencies among the three CS groups, and the results indicated different CNV mutations on various chromosomes (Fig. [126]4F–H). In summary, our drug prediction and mutation analysis of GC subtypes sheds light on the complex genetic and molecular landscape that governs drug sensitivity. These insights are pivotal for the rational design of personalized treatment strategies and may pave the way for more effective, subtype-specific therapies to treat GC. Fig. 4. [127]Fig. 4 [128]Open in a new tab Drug prediction and mutation analysis of GC subtypes. (A-C) CS2 shows the poorest sensitivity to cisplatin (A), docetaxel (B), and gemcitabine (C). D Differences in the genomic landscape of the three CS subtypes (from top to bottom), shown by tumor mutation burdens (TMBs), relative contributions of four mutation signatures, and selected differentially mutated genes; E GSEA differential enrichment analysis of the three CS subtypes; F–H CNV GISTIC score mutation profile for CS1 (F), CS2 (G), and CS3 (H) Single-cell landscape of GC The immune landscape of GC at the single-cell level was further investigated. After dimensionality reduction and clustering of the [129]GSE183904 database, 28 clusters were delineated (Fig. [130]5A). Through meticulous annotation, 10 cell subtypes including T cells, B cells, plasma cells, mast cells, monocytes, macrophages, fibroblasts, endothelial cells, epithelial cells, and pericytes were identified. Using the infercnv tool, malignant epithelial cells were then identified (Fig. S2A), expanding our cellular repertoire to 11 distinct subtypes (Fig. [131]5B). Subsequently, each cell subtype was subjected to functional enrichment analysis to validate the accuracy of our annotations (Fig. [132]5C). Then, data were extracted from tumor samples, and key differentially expressed genes among transcriptional subtypes were intersected with platinum resistance genes from the [133]GSE186205 data to identify key platinum resistance–related genes. Enrichment analysis based on five single-cell scoring methods revealed that these genes were most significantly enriched in malignant cells (Fig. [134]5D). Finally, based on the median drug resistance score, malignant cells were categorized into two drug-resistant groups (DRGs): one showing high levels of resistance and other showing low levels of resistance. Fig. 5. [135]Fig. 5 [136]Open in a new tab Single-cell landscape of GC. A Dot plot for single-cell annotation; B UMAP of single-cell subtypes; C Heatmap of marker genes for each cell type; D Dot plot for the platinum resistance gene set scoring in single cells cNMF evaluation of key modules for platinum resistance in malignant cells Malignant cell subpopulations in GC were analyzed via cNMF, and six distinct feature modules were identified, along with another module referred to as “other” (Fig. [137]6A). Cell markers for the subpopulations facilitated their unequivocal identification (Supplementary Fig. [138]2B). Subsequently, SCENIC was used to compute TFs for each module, revealing that the M1 subpopulation was significantly enriched in the high DRG group but was absent in the low DRG group (Fig. [139]6B). The M1 subpopulation was also shown to have the highest scores (Fig. [140]6C). In addition, the R/O index confirmed the significant enrichment of M1 in the group with high resistance (Fig. [141]6D). Overall, these findings suggest that M1 is a key subpopulation driving platinum resistance in GC. Finally, SCENIC analysis revealed significant transcriptional differences between modules (Fig. [142]6E), with M1 overexpressing TFs such as KLF9, FOS, and ATF3, which are known to mediate inflammatory responses. Fig. 6. [143]Fig. 6 [144]Open in a new tab Identification of key modules for platinum resistance. A Identification of key malignant cell subpopulations using cNMF. B Bipartite heatmap displaying differences in the activities of identified transcription factors across groups. The top 15 regulators with the highest regulator-specific score are plotted for each module. The proportion of active regulators differed by at least 20% between groups, and the AUC score differences of the regulator-specific scores between groups are significant. C Violin plot of DRG scoring. D Differences in distribution between the high- and low-resistance DRGs as revealed by the R/O index. E Top 10 differential transcription factors in M1–6 modules Differences in cell communication between subpopulations with high and low drug resistance Cell communication dynamics between subpopulations with high or low drug resistance and other cells were investigated to elucidate the mechanisms underlying drug resistance and identify potential therapeutic targets. The results revealed that malignant cells with high resistance DRG exhibited significant communication with endothelial cells, whereas malignant cells with low-resistance DRG showed weaker interactions (Fig. [145]7A). In addition, pathways related to GRN, KIT, SEMA3, LIFR, and complement were highly enriched in malignant cells with high DRG, and proliferation and inflammatory pathways such as VEGF, SPP1, and TNF were also highly expressed in this group (Fig. [146]7B, C). When malignant cells acted as ligands in cell communication, those in the high DRG group showed higher expression of VEGFA-VEGFR1/2 and SEMA3B-(NRP1/2 + PLXNA2)-related ligands with endothelial cells and MIF-ACKR3-related ligands with stromal cells than the malignant cells in the low DRG group (Fig. [147]7D). Conversely, when malignant cells acted as receptors, the MDK-related ligands were only expressed between the high DRG malignant cells and endothelial cells, highlighting the important role of endothelial cells in tumor resistance (Fig. [148]7E). Fig. 7. [149]Fig. 7 [150]Open in a new tab Differences in cell communication between subpopulations with high and low drug resistance. A Differences in cell communication based on ligand pair counts. This panel highlights the differences in the number of ligand–receptor pairs involved in cell communication between subpopulations with high and low drug resistance. B Differences in the strength of cell communication pathways between subpopulations with high and low drug resistance. This panel shows how the intensities of communication pathways vary between the two groups. C Heatmap of cell communication showing the differential expression of communication pathways between subpopulations with high and low drug resistance, highlighting key differences in cellular interactions. D Ligand-based cell communication in malignant cells. This panel shows the ligand expression patterns in malignant cells acting as ligands in cell communication, highlighting the significant ligands involved in interactions between malignant cells with high drug resistance and other cell types. E Receptor-based cell communication in malignant cells. This plot shows the expression of receptor-related ligands in malignant cells acting as receptors in cell communication, emphasizing interactions between malignant cells with high drug resistance and endothelial cells Immune factor regulation and spatial localization of subpopulations with high and low drug resistance The differences in immune factor regulation between subpopulations with high and low drug resistance were further explored, The CS2 subtype demonstrated marked overexpression of immunomodulatory ligands (CCL5, VEGFB, TGFB1), coupled with global hypomethylation relative to other subtypes. Comparative genomic analysis revealed elevated somatic copy-number alteration frequencies (amplifications/deletions) in CS3 versus other subgroups, whereas CS2 exhibited reduced amplification frequency. These findings collectively suggest stromal enrichment in CS2 may modulate tumor cell plasticity through microenvironmental crosstalk (Fig. [151]8A). Spatial transcriptomics analysis revealed that the high resistance DRG tended to cluster locally and was surrounded by wall cells, suggesting the potential formation of niches in these malignant cells within the tumor. These cells were also associated with immune cell infiltration, particularly T cells, at their periphery (Fig. [152]8B–E). Subsequent spatial cell communication analysis also indicated significant differences in the levels of communication between malignant cells with high or low drug resistance and other cell types, with differential expression of ligand–receptor pairs (Supplementary Fig. [153]3A-D). Fig. 8. [154]Fig. 8 [155]Open in a new tab Immune factor regulation and spatial localization of subpopulations with high and low drug resistance. A Regulation of immune modulatory factors. From left to right: mRNA expression; relationship between expression and methylation; amplification frequency; and deletion frequency of 75 immune subtype genes. B Spatial distribution heatmap of malignant cells with high drug resistance in the [156]GSM7990475 sample. C Spatial distribution of cell subpopulations in the [157]GSM7990475 sample. D Spatial distribution heatmap of malignant cells with high drug resistance in the [158]GSM7990482 sample. E Spatial distribution of cell subpopulations in the [159]GSM7990482 sample Construction of a prognostic model of platinum resistance and evaluation of its effectiveness A prognostic model was constructed based on differentially expressed genes associated with platinum resistance. Using bootstrap kernel-based univariate Cox regression analysis, 28 key prognostic genes were identified (Fig. [160]9A and Supplementary Table [161]1). Subsequently, the prognostic model was built using 101 different algorithms. RSF was selected as the best model based on its C-index (0.711), which was significantly higher than that obtained from the second-best model (0.662) (Fig. [162]9B). The model components included SOX9, TPM2, KLF9, MYL9, PRR15L, CHST15, CHST3, AHNAK2, SOD3, MYADM, and THBS2. Furthermore, the selected model showed superior prognostic performance than other clinical indicators in both the training and validation cohorts, as evidenced by the higher C-index value (Fig. [163]9C–E). Patients that were assigned a higher risk score exhibited significantly worse prognosis (p < 0.001) (Fig. [164]9F–H). The prognostic AUC curve showed 1-, 3-, and 5-year AUC scores of 0.99, 0.99, and 0.95, respectively (Fig. [165]9I). In addition, as the risk score increased, the hazard ratio (HR) for the patients also gradually increased (Fig. [166]9J). Calibration analysis demonstrated good prognostic performance of the model at 1, 3, and 5 years (Fig. [167]9K). Decision curve analysis further confirmed the superior performance of the model over other clinical indicators (Fig. [168]9L–N). Finally, the nomogram further validated the model’s effectiveness (p < 0.01) (Fig. [169]9O). These results collectively validated the robustness and clinical utility of our prognostic model of platinum resistance, highlighting its potential to refine risk stratification and assist in the development of personalized therapeutic approaches for GC treatment. Fig. 9. [170]Fig. 9 [171]Open in a new tab Construction of the prognostic model and evaluation of its performance. A Forest plot of prognosis-related genes, with risk-associated genes and protective genes indicated in red and blue, respectively. B Screening of the optimal predictive model using 101 different algorithms. C–E Comparison of the C-index values of the models based on the TCGA dataset (C), [172]GSE62254 dataset (D), and [173]GSE15459 dataset (E). F–H Kaplan–Meier survival curves for high- and low-risk groups based on the TCGA dataset (F), [174]GSE62254 dataset (G), and [175]GSE15459 dataset (H). I AUC scores of the selected model at 1, 3, and 5 years. J Restricted cubic spline plot demonstrating HR in the Cox proportional hazards model as a function of continuous variables. K Decision curve analysis of the model. L–N Calibration curve for the model and other clinical factors at 1 year (L), 3 years (M), and 5 years (N). O Prognostic nomogram for survival prediction. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001) Evaluation of the prognostic model for immunotherapy and drug response prediction The IMvigor210 cohort was employed to assess the utility of the predictive model in predicting responses to immunotherapy. The results demonstrated that patients in the high DRG group had significantly worse OS and poorer prognoses in both disease stages I–II and III–IV than those in the low DRG group (p < 0.05) (Fig. [176]10A–C). Moreover, the proportion of progressive disease or stable disease (SD) cases was higher for patients exhibiting malignant cells with high drug resistance (p < 0.001, Fig. [177]10D, E), indicating poorer immunotherapy efficacy in them. Submap analysis further revealed a higher proportion of SD in the low DRG group, suggesting better immunotherapy efficacy in this group (Fig. [178]10F). Drug response prediction based on the CTRP and PRISM databases suggested that patients exhibiting high resistance DRG exhibited poorer responses to drugs such as dasatinib and erlotinib (Fig. [179]10G, H), indicating potential resistance to these and possibly other drugs. In addition, analysis of TMB revealed higher mutation frequencies in patients with high resistance DRG (Supplementary Fig. [180]4 A), while no significant differences in tumor neoantigen burden (TNB) were observed between the groups (Supplementary Fig. [181]4B). Fig. 10. [182]Fig. 10 [183]Open in a new tab Immunotherapy and drug response predictions. A Overall survival (OS) is significantly worse in the high DRG group within the IMvigor210 cohort. B Prognosis is poorer for the high DRG group at disease stages I–II (C) and stages III–IV (D) in the IMvigor210 cohort. Differences in immunotherapy response rates between the high and low DRG groups in the IMvigor210 cohort. E Bar chart illustrating immunotherapy response rates in the high- and low DRG groups within the IMvigor210 cohort. F Submap-based prediction of immunotherapy responses in the high- and low DRG groups. G, H Drug response predictions based on the CTRP database (G) and PRISM database (H) Functional enrichment and immune infiltration analysis Pancancer enrichment analysis revealed a significant enrichment of inflammation-associated pathways in high DRG tumors, including EMT, JAK-STAT, TNF-α, and KRAS signaling pathways (Fig. [184]11A). These findings suggested that such tumors may harbor niche-like features with prominently expressed inflammatory pathways. Single-gene pathway enrichment analysis demonstrated that genes such as KLF9 were also enriched in pathways related to EMT, TNF-α signaling, immune responses, and JAK-STAT signaling (Fig. [185]11B). Immune infiltration analysis at the single-gene level revealed a positive correlation between most genes (excluding SOX9 and PRR15L) and cytotoxic immune cells (Fig. [186]11C). By intersecting model genes with M1 module TFs, KLF9 was identified as a candidate gene. Spatial localization analysis further revealed that KLF9 was highly enriched in malignant cell subpopulations within the high DRG group, suggesting that KLF9 may serve as a critical mediator of drug resistance. Fig. 11. [187]Fig. 11 [188]Open in a new tab Functional enrichment and immune infiltration analysis. A Pancancer functional enrichment analysis based on the model. B Functional enrichment analysis for a single gene. C Immune infiltration analysis for a single gene. D Spatial distribution heatmap of KLF9 + malignant cells in the [189]GSM7990482 sample. E Spatial distribution of KLF9 + cell subpopulations in the [190]GSM7990482 sample Overexpression of KLF9 inhibits the progression of GC and reduces chronic inflammation The inhibition of KLF9 reduced chronic inflammation during HCC progression. Subsequently, the analysis of the TCGA database revealed that KLF9 exhibits high sensitivity and specificity, whereas TCGA and GEO data mining revealed low expression of this gene (Fig. [191]12A). Given the significant resistance of the AGS cell line to cisplatin, it was hypothesized that the reduced expression of KLF9 may be a key determinant of this chemoresistant phenotype. Surprisingly, KLF9 was shown to be downregulated in GC (Fig. [192]11A). Overexpression of KLF9 resulted in the inhibition of GC invasion and proliferation (Fig. [193]12B–D). Moreover, the mRNA expression levels of inflammatory cytokines were significantly reduced upon KLF9 overexpression, with notable inhibition of IL-6, IL-1β, and TNF-a (Fig. [194]12E). These findings revealed KLF9 as a critical regulator of GC progression, with implications for both inflammation management and tumor suppression, thereby reinforcing its potential as a therapeutic target in GC treatment. Fig. 12. Fig. 12 [195]Open in a new tab Overexpression of KLF9 inhibits GC progression. A Expression and diagnostic value of KLF9. B–D Plate cloning experiment (B), Transwell assay (C), and cell scratch assay (D) revealing the effect of KLF9 overexpression on GC. E KLF9 overexpression inhibits the expression of chronic inflammatory factors Discussion In the present study, GC was successfully classified into three distinct subtypes, i.e., CS1, CS2, and CS3, using the SNF algorithm and integrating GC-related genes, methylation sites, and somatic mutation data associated with platinum resistance. These subtypes exhibit distinct biological characteristics: CS1 is characterized by cell proliferation and immune deficiency, CS2 demonstrates the upregulation of carcinogenic pathways, such as ECM remodeling and the PI3K-Akt pathway, and CS3 shows enhanced activity in immune-related pathways. Among these subtypes, CS2 shows the worst prognosis, a finding validated in independent datasets. This highlights the critical role of the ECM in tumor resistance. The ECM, a key noncellular component of the tumor microenvironment [[196]33] (TME), plays a pivotal role in maintaining tissue architecture and function. As a physical barrier, the ECM can dissolve drugs or delay their delivery, promoting tumor resistance by modulating drug transport, facilitating immune evasion, and influencing transmembrane signaling processes [[197]36]. On a physical level, an enriched stromal matrix obstructs drug delivery and disrupts tumor vasculature, contributing to therapeutic resistance [[198]37]. In addition, the ECM in the interstitial space increases fluid pressure due to the physical constraints posed by the tumor mass, significantly impairing drug efficacy [[199]38], and excessive cancer cell proliferation enhances fluid flow from the tumor to the surrounding tissue, further hindering drug transport. Our subsequent analysis revealed that model-related genes were highly enriched in pathways related to epithelial–mesenchymal transition (EMT) and oxidative stress. Spatial transcriptomic analysis indicated that malignant cells exhibiting high DRG were spatially localized, whereas immune cell subpopulations were distributed peripherally without infiltrating these regions, underscoring the ECM-enriched characteristics of the CS2 subtype. Drug sensitivity predictions further confirmed the resistance of CS2 cells to cisplatin, docetaxel, and gemcitabine. In addition, mutation frequency analysis revealed that CS2 exhibited a lower mutation frequency than the other two subtypes, suggesting that CS2-mediated resistance is primarily driven by ECM-induced drug sequestration rather than genetic mutations. Identifying these patients and targeting key genes to inhibit stromal matrix production may significantly improve OS in this patient group. Several therapeutic strategies have been proposed to overcome the treatment challenges posed by the ECM—rich tumor microenvironment of the CS2 subtype in gastric cancer. Targeting the ECM directly is a promising approach. For instance, hyaluronidase (HAase), an enzyme that degrades hyuralonic acid, a major component of the ECM, has shown potential in clinical trials to improve drug delivery and tumor uptake in ECM—rich cancers [[200]50, [201]51]. Additionally, inhibitors of collagen crosslinking, such as lysyl oxidase (LOX) inhibitors, may reduce the physical barriers imposed by the ECM, thereby enhancing drug penetration. Another strategy involves modulating the EMT process, as EMT is closely associated with ECM remodeling and drug resistance. Drugs that inhibit EMT may synergize with conventional chemotherapies to enhance treatment efficacy. Furthermore, addressing the interstitial fluid pressure (IFP) elevated by the ECM—rich stroma is crucial [[202]32]. Angiopoietin—2 inhibitors have been shown to normalize tumor vasculature and reduce IFP, potentially improving drug distribution within tumors [[203]2, [204]39]. Given the CS2 subtype's characteristics, such as its ECM—enriched environment and spatial localization of malignant cells, therapies targeting the ECM and related processes are likely to be particularly effective. These approaches have the potential to significantly enhance treatment outcomes for patients with the CS2 subtype of gastric cancer, offering a promising avenue for future research and clinical translation. Further single-cell analyses showed that patients affected by the CS2 subtype exhibited significantly high concentrations of M1 macrophages [[205]12], particularly those with high DRG, consistent with our hypothesis. The M1 subpopulation showed high expression of TFs such as KLF9, FOS, and ATF3, which mediate inflammatory responses and participate in ECM remodeling. Cell communication analysis identified the MIF–ACKR3 receptor–ligand pair as a key mediator of interactions between stromal and tumor cells. Previous studies have demonstrated that MIF–ACKR3 plays a role in the crosstalk between esophageal cancer and fibroblasts [[206]39, [207]40] as well as in hematologic malignancies, where it facilitates tumor–stroma communication and promotes tumor progression. Our prognostic model confirmed the association of high DRG expression with worse survival outcomes and reduced response to immunotherapy [[208]41, [209]42], a result that was further validated in the IMvigor210 cohort. Interestingly, the gene encoding KLF9, a key TF in the M1 module, was identified as one of the model genes. KLF9, a key member of the KLF family, plays a significant role in the proliferation and metastasis of GC and other malignancies [[210]9]. Studies have shown that KLF9 is closely associated with tumor progression. It can markedly enhance GC cell proliferation and increase the expression of pro-inflammatory cytokines. Moreover, silencing KLF9 can suppress the resistance of GC cells to platinum-based therapy, indicating that KLF9 may act as a pro-carcinogenic factor in GC [[211]40]. However, some studies have also revealed that KLF9 can function as a tumor suppressor in certain types of cancer. For instance, in hepatocellular carcinoma, KLF9 can inhibit tumor metastasis by regulating the epithelial-mesenchymal transition (EMT) process [[212]40]. This suggests that the role of KLF9 in GC may be complex and its function may be influenced by the tumor microenvironment and cell type. Regarding inflammatory pathways, KLF9 may influence GC progression by modulating the NF-κB signaling pathway. NF-κB plays a central role in inflammatory responses and immune reactions, and its abnormal activation is linked to the development of various cancers, including GC [[213]6]. KLF9 may interact with NF-κB p50/p65 and bind to the GT box on the MMP9 promoter, thereby inhibiting MMP9 expression and affecting tumor cell invasion and metastasis [[214]3, [215]52]. In addition, KLF9 may also regulate the STAT3 signaling pathway. STAT3 is crucial for inflammatory and immune responses, and can promote tumor cell proliferation, survival, and angiogenesis. Although there is limited research on KLF9 and STAT3 specifically in GC, studies in other cancers have shown that KLF9 can negatively regulate STAT3 activity and suppress tumor progression [[216]30]. KLF9 may also be involved in regulating the inflammatory microenvironment. Its expression may be regulated by cytokines in the inflammatory microenvironment, and it can, in turn, influence the behavior and function of cells within this microenvironment [[217]19]. For example, in other tumor models, KLF9 can modulate the polarization of tumor-associated macrophages, shifting them from the pro-tumoral M2 type to the anti-tumoral M1 type. This implies that KLF9 may participate in reshaping the inflammatory microenvironment in GC through a similar mechanism, thereby affecting tumor progression and prognosis. In-depth research on the role of KLF9 in GC and its interaction with inflammatory pathways is crucial. It can help uncover the mechanisms underlying GC development, identify new therapeutic targets, and improve the survival rate and quality of life for GC patients. Further experimental studies and clinical trials are needed to validate and expand these findings, paving the way for the application of KLF9 in GC therapy. This study has several limitations. First, while we elucidated mechanisms underlying the poor prognosis of the CS2 subtype, the biological distinctions between CS1 and CS2 require deeper characterization. Second, although cross-cancer analyses suggested the clinical potential of our immunotherapy prediction model, inherent biological disparities in tumor microenvironment composition, key immune evasion mechanisms, and therapeutic response patterns may limit direct translation to GC. Finally, the precise mechanisms by which KLF9 mediates platinum resistance in GC warrant further investigation. To address these gaps, we will prospectively validate findings in a dedicated GC treatment cohort and employ internal clinical specimens to systematically explore KLF9's dual roles in platinum resistance and immunotherapy responsiveness, supported by functional validation using patient-derived models. Conclusion In summary, a prognostic model based on genes associated with platinum resistance was developed by integrating bulk transcriptomics, spatial transcriptomics, and single-cell data. Multiomics analysis was used to systematically explore how key GC subpopulations mediate platinum resistance. Experimental validation of critical genes further supports our findings, offering novel insights for the development of therapeutic targets to overcome platinum resistance in GC. Supplementary Information Below is the link to the electronic supplementary material. [218]Supplementary Figure 1 (PNG 577 KB)^ (576.9KB, png) [219]Supplementary Figure 2 (PNG 1711 KB)^ (1.7MB, png) [220]Supplementary Figure 3 (PNG 607 KB)^ (606.7KB, png) [221]Supplementary Figure 4 (PNG 113 KB)^ (112.1KB, png) [222]Supplementary Table 1 (XLSX 13 KB)^ (12.2KB, xlsx) Abbreviations GC Gastric cancer SNF Similarity network fusion STAD Stomach adenocarcinoma DRGs Drug resistance genes CNV Copy number variation MOVICS Multi-Omics integration and VIsualization in Cancer Subtyping CPI Clustering prediction Index TFs Transcription factors NTP Nearest template prediction TMB Tumor mutation burden TNB Tumor neoantigen burden scRNA-Seq Single-cell RNA sequencing RSF Random survival forest Enet Elastic net plsRcox Cox Partial least squares regression C-index Harrell’s Concordance Index CCLs Human cancer cell lines CCLE The Broad Institute Cancer Cell Line Encyclopedia TME Tumor microenvironment EMT Epithelial–mesenchymal transition HAase Hyaluronidase LOX Lysyl oxidase Author contributions PZ, HI, and LW conceived the study. PZ, LW, KX, YX, HL, KC, YH, JZ, RY, HC, HS, and YH drafted the manuscript. YH and LW performed literature search and collected the data. LW completed in vitro experiments. LW analyzed and visualized the data. PZ, HI, and LW helped with the final revision of this manuscript. All authors reviewed and approved the final manuscript. Funding This study was supported by the open competition mechanism to select the best candidates for key research projects of Ningxia Medical University (XJKF240313). Availability of data and materials The datasets analyzed in this study can be found in GEO ([223]https://www.ncbi.nlm.nih.gov/geo/) and Xena ([224]https://xena.ucsc.edu/). The datasets used and analyzed in this study are available from the corresponding authors of this study upon reasonable request. Declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Competing interests The authors declare no conflict of interest. Footnotes Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Pengcheng Zhang and Lexin Wang have contributed equally to this work. Contributor Information Hang Song, Email: hangsong@ahtcm.edu.cn. Peng Wang, Email: wisdom-wp@tmmu.edu.cn. Yajuan Fu, Email: 20210202@nxmu.edu.cn. References