Integrated single-cell sequencing, spatial transcriptome sequencing and bulk RNA sequencing highlights the molecular characteristics of parthanatos in gastric cancer

Sample source and parthanatos-related gene source

Based on the TCGA platform, we downloaded and curated multi-omics data of pan-cancer in humans, including CNV, SNV, and methylation data, as well as mRNA expression profiles and corresponding clinical information at the transcriptome level. The processing methods for these data were similar to previous studies [17, 18]. Visualization of the results was achieved using the R language and TBtools software. In addition to the pan-cancer data, we focused on in-depth analysis and exploration of the transcriptome data of GC. Specifically, we downloaded the corresponding data from the TCGA-STAD cohort, which included a total of 407 samples, including 32 adjacent non-cancer samples and 375 GC samples. After excluding samples with a survival time of less than 30 days, we analyzed the remaining 312 GC samples with complete follow-up information. Furthermore, to obtain more accurate and convincing results, we collected additional publicly available GC transcriptomic data (GSE84437). After integrating the prognostic information and expression profiles, we obtained a total of 431 GC samples with complete follow-up information. Since different data sources may have batch effects, we performed batch correction using the methods from previous studies [19, 20]. GeneCards is a comprehensive online gene information database that provides detailed information about human genes, including gene function, expression, disease associations, mutations, protein information, and drug association [21]. Users can easily search and browse thousands of genes to gain a deeper understanding of their important roles in biology, medicine, and drug development. All parthanatos-related genes in this study were obtained from the GeneCards platform, where we retrieved 32 parthanatos-related genes for subsequent in-depth analysis.

Analysis of single-cell sequencing and spatial transcriptome data

The single-cell GC dataset is derived from GSE163558, which includes 3 GC tumor samples (GSM5004180, GSM5004181, and GSM5004182), 1 tumor-adjacent sample (GSM5004183) and 6 metastasis sample (GSM5004184, GSM5004185, GSM5004186, GSM5004187, GSM5004188, GSM5004189) [22]. In order to perform quality control on our raw count data, the following criteria were set: A) min.cells = 3, min.features = 200, B) nCount_RNA >= 1000, C) 200 <= nFeature_RNA <= 8000, D) percent.mt <= 20. To address the issue of differing total counts among cells, we applied a global scaling normalization method (LogNormalize) to normalize the data. This involved normalizing the feature expression measurements of each cell by their total expression, followed by multiplication by a scaling factor (10000), and finally transforming the values using the natural logarithm. Based on the “vst” algorithm, we identified feature genes that exhibited high intercellular differences in the dataset and normalized all genes for further analysis.

To simplify computations and remove data noise, we performed a joint analysis of PCA and Harmony, aiming for batch correction and dimensionality reduction. Harmony applies principal component analysis to embed the transcriptomic expression profile into a lower-dimensional space and then iteratively removes dataset-specific effects. Based on the study of Jiang et al., we manually annotated the single-cell data after quality control and identified a total of 7 cell types, including myeloid cell, epithelial cell, T cell, stromal cell, proliferative cell, B cell, and NK cell [22]. We visualized the above cell subgroups in the form of UMAP [23]. Additionally, we employed six algorithms to score gene sets in the single-cell dataset: AUCell [24], Ucells [25], Singscore [26], ssGSEA [27], addmodulescore [28, 29], and scoring [17]. The scoring value is the sum of scores from the previous five methods and serves to assess the overall distribution of parthanatos gene set scores more stably and comprehensively.

The spatial transcriptomics data were obtained from the GSE186290 dataset. It is important to note that this dataset consists of spatial transcriptomic sequencing of tissue samples from GC mice, and there are currently no publicly available human GC spatial transcriptomics data for analysis. The Read10X_Image function was used to read the spatial distribution information of tissue images and cells, while the Load10X_Spatial function integrated the spatial expression profile with spatial localization information. Similar to the single-cell analysis, we used the aforementioned six gene set scoring methods to evaluate the parthanatos scores of each cell. Finally, we used the SpatialFeaturePlot function to visualize the spatial distribution of parthanatos scores.

Construction of parthanatos molecular classifier for GC

As research progresses, the scientific community is gradually realizing the significant molecular heterogeneity both between and within tumors. It is due to this heterogeneity that individuals with the same disease often have different treatment response strategies and clinical outcomes. In light of this, we have developed a novel molecular classifier for GC based on parthanatos-related genes. Our goal is to clearly distinguish GC patients with different parthanatos features.

First, we merged the GC transcriptomic data from the TCGA-STAD cohort and the GSE84437 cohort, resulting in a total of 743 GC samples. Molecular clustering analysis was performed using the “ConsensusClusterPlus” R package developed by Wilkerson et al. [30], with the following specific parameters: reps=50, pItem=0.8, clusterAlg=“km”, distance=“euclidean”, maxK=9. The optimal number of clusters was determined based on the consensus cumulative distribution function and delta area plot. The clustering of GC patients was based on unsupervised clustering of parthanatos-related genes. To evaluate the classification performance and clinical relevance of this classifier, we further assessed the parthanatos scores and clinical prognostic differences among different subtypes of GC patients. Additionally, we generated a heatmap to visually depict the expression characteristics of each parthanatos-related gene in this molecular classifier.

Identification of internal molecular characteristics of the parthanatos molecular classifier

To fully explain the different clinical outcomes among GC patients with different molecular subtypes, we conducted in-depth research and exploration of their intrinsic molecular features. Firstly, based on the KEGG database, we collected 42 classical metabolic pathways and 24 classical immune pathways. The “GSVA” package was used to evaluate the metabolic and immune signaling strength in the 743 GC samples. Finally, we depicted the distribution of the intensity of each metabolic and immune signal in the form of a heatmap. In general, metabolic reprogramming and immune microenvironment are classical molecular markers of tumors, and different forms of metabolic reprogramming and immune microenvironment features may be potential reasons for the different prognostic outcomes of different parthanatos subtypes.

Apart from the immune pathways, there are many other algorithms in bioinformatics that can evaluate the tumor immune microenvironment. Therefore, we have subsequently conducted various immune-related algorithms to explore the inherent relationship between the parthanatos subtype and the tumor immune microenvironment. A) The “Estimate” package is used to predict the content of stromal cells and immune cells in malignant tumor tissues based on gene expression data. This algorithm is based on enrichment analysis of individual sample gene sets and produces four scores: a) stromal score (indicating the presence of stroma in tumor tissue) b) immune score (representing the infiltration of immune cells in tumor tissue) c) ESTIMATE score (stromal score + immune score) d) tumor purity. Using the “Estimate” package, we calculated the aforementioned four scores for 743 GC patients and conducted corresponding comparisons. B) TIMER2.0 online platform provides seven algorithms for predicting immune cell infiltration in tumor samples, including: TIMER, CIBERSOFT, CIBERSOFT-ABS, QUANTISEQ, XCELL, EPIC, and MCPCOUNTER [31]. The specific procedure is as follows: based on R language, the 743 GC samples were divided into 5 groups, with each group containing 150 samples, and the last group containing 143 samples. The expression matrix of these five groups was uploaded to the Immune Estimation module of the TIMER2.0 platform for prediction, and the results of immune cell infiltration were then downloaded and merged. Heatmaps were generated to show the immune cell infiltration of each sample and calculate the corresponding statistical differences [32]. C) Immune checkpoints are the main limiting factors for immune cells to exert anti-tumor functions. Therefore, we collected classical and recognized immune checkpoints from previous literature [33], and compared the expression characteristics of immune checkpoint-related genes in different parthanatos subtypes of patients. D) The association between each parthanatos-related gene and the GC immune microenvironment was calculated. Based on previous literature reports [34], we identified a set of 29 genes related to immune cells and immune-related functions. The ssGSEA algorithm was used to evaluate the immune cell and immune function scores of 743 GC samples, and Spearman correlation analysis was performed to explore the potential association between each parthanatos-related gene and immune cell infiltration and immune function. At the same time, the correlation between parthanatos score and immune cell infiltration and immune function was calculated.

Identification of potentially sensitive drugs for GC patients based on parthanatos molecular classifier

Danielle Maeser and colleagues developed a novel tumor drug sensitivity prediction scheme, namely “oncoPredict” package [35]. This R package connects in vitro and in vivo drug screening, allowing easy prediction of tumor response to a large number of drugs screened in cancer cell lines. In this study, we utilized the “oncoPredict” package and the GDSC2 dataset to predict the response of 743 GC patients to each drug and ultimately analyzed potential beneficial drugs for different parthanatos subtypes of GC patients.

Developing a novel GC prognostic assessment model based on parthanatos molecular classifier

Based on the “limma” package, differential gene expression between different parthanatos subtypes was analyzed. The “clusterProfiler” package was used to assist in Gene Ontology (GO) enrichment analysis and KEGG pathway enrichment analysis [36–38]. Subsequently, the differentially expressed genes were used to construct a GC prognostic model. Firstly, 431 GC samples from the GSE84437 dataset were randomly divided into a 6:4 ratio. 60% of the GC samples were considered as the training cohort for the model (260 GC samples), while the remaining 40% of the GC samples were considered as an internal validation set 1 (test1 cohort, consisting of 171 GC samples). Additionally, all samples from the GSE84437 dataset were treated as an internal validation set 2 (test2 cohort, consisting of 431 GC samples). The GC samples from the TCGA-STAD dataset were considered as the external validation set for the model (test3 cohort, consisting of 312 GC samples).

Firstly, in the training cohort, the model was constructed through the following steps: A) Univariate Cox regression was used to select genes that are related to GC prognosis. B) LASSO regression was applied to filter genes in order to avoid gene collinearity and overfitting of the model [39]. C) Multivariate Cox regression was used to build the prognostic model (based on the predict function to calculate the risk value for each GC patient). D) Patients were divided into high-risk and low-risk groups based on the median risk value. E) The “survival” and “survminer” packages were used to evaluate the modeling effect of the prognostic model by comparing the survival differences between different risk groups of GC patients [40]. F) The “survivalROC” package was utilized to plot the ROC curve and calculate the AUC value to evaluate the predictive accuracy of the model [41]. Subsequently, in the internal validation set 1, internal validation set 2, and external validation set, the same model genes selected from the training cohort were subjected to multivariate Cox regression. The predict function was used to predict the risk value for each sample. Based on the median risk value from the training cohort, all samples from the different datasets were grouped accordingly. Finally, the stability of the model was verified through internal and external validation strategies. To further facilitate the wide application of the prognostic model in clinical settings, clinical features were introduced as variables in the model [42]. Based on the “rms” package, a nomogram was constructed, integrating the patient’s risk group, grade, stage, gender, and age [32]. Calibration plots and ROC curves were used to evaluate the predictive accuracy of the nomogram.

Differences in immune microenvironment between high and low risk group

To investigate the potential reasons for the different outcomes between high- and low-risk groups, we subsequently conducted immune microenvironment analysis on samples from both groups. Based on the immune-related scores obtained during molecular clustering, we compared the corresponding differences between the high- and low-risk groups. Specifically, we used the TIMER2 platform to predict the immune cell infiltration abundance using seven immune algorithms, and the statistical differences in immune cell infiltration between the high- and low-risk groups were evaluated using the wilcox.test function. Additionally, we also compared the expression of immune checkpoint genes between different risk groups.

Availability of data and materials

The datasets analyzed in this work may be found in the Supplementary Materials or the first author may be contacted.

Pan-cancer analysis based on genomics and transcriptomics

To study the variations and expression changes of parthanatos-related genes in various human tumors, we determined the mRNA expression as well as the frequency of copy number variations (CNVs) and single nucleotide variations (SNVs) of parthanatos-related genes in different tumors, using sample data from the TCGA database. Genes such as PARP1, NAMPT, AIMP2, MCL1, and TOMM20 exhibited widespread amplifications of CNVs in multiple cancers (Figure 1A), while genes like RNF146, GPX4, ESR1, CUL4A, and PTEN showed widespread deletions (Figure 1B). Figure 1C shows the SNV status of parthanatos-related genes, with genes such as PTEN, DDB1, ESR1, PARP1, and AIFM1 displaying high levels of SNVs in multiple cancers, especially in UCEC, STAD, and SKCM, providing valuable guidance for subsequent experimental research (Figure 1C). In terms of expression, genes like PARP1, FEN1, and GAS5 showed widespread high expression in cancers such as BRCA, LUAD, and LUSC, while genes like PTUD1, ESR1, and ESR2 showed widespread low expression (Figure 1D). Based on survival-related data, we also summarized the risk value of parthanatos-related genes in various cancers (Figure 1E). Additionally, Figure 1F shows the differential methylation levels of parthanatos-related genes between cancer and normal tissues. As shown in the figure, genes like RIPK1, RAB33A, and ESR1 displayed significantly higher methylation levels in cancer tissues compared to normal tissues (Figure 1F).

Figure 1. Pan-cancer analysis of parthanatos-associated genes. (A) CNV amplification data of parthanatos-related genes in different tumor types and the fan color represents the amplification frequency. (B) CNV deletion data of parthanatos-related genes in different tumor types and the fan color represents the deletion frequency. (C) SNV mutation data of parthanatos-related genes in different tumor types and the fan color represents the frequency of SNV. (D) In the expression data of parthanatos-related genes in different tumor types, the color of the squares represents the value of log2 (FC), and the size of the squares represents the value of -log2 (FC). (E) The risk profile of parthanatos-related genes in different tumor types, with pink representing risky, green representing protective, and gray representing no statistical difference. (F) The comparison of methylation of parthanatos-related genes between different tumor types and normal tissues, the color of the circle represents the methylation difference, and the size of the circle represents the statistical significance.

Single cell transcriptomic analysis of GC

Our single-cell data consists of 3 GC tumor samples (GSM5004180, GSM5004181, and GSM5004182), 1 tumor-adjacent sample (GSM5004183) and 6 metastasis sample (GSM5004184, GSM5004185, GSM5004186, GSM5004187, GSM5004188, GSM5004189). To perform quality control on our raw count data, we first filtered out unsatisfactory cells based on sequencing depth, number of genes, mitochondrial content, and ribosomal content, there were 53940 cells before quality control, and 41264 cells were left after quality control (Supplementary Figure 1). Subsequently, we conducted a joint analysis of PCA and Harmony to remove batch effects and reduce dimensionality, ultimately obtaining 20 principal components (PCs) for further analysis. During the determination of the resolution value, we observed that increasing the resolution value allowed us to clearly observe which cell clusters were continuously dividing into subclusters, revealing the relationships between cell clusters at different resolutions. When the resolution was set to 2, we observed significant interweaving between cells, so we chose a resolution of 2 and obtained 34 unknown cell clusters (Figure 2A), and the specific maker genes expressed in each cell cluster were shown in bubble charts (Figure 2B). Then we combined all the single-cell samples, and performed dimensionality reduction on the merged samples using tSNE nonlinear clustering algorithms (Supplementary Figure 2), and we displayed 34 unknown cell clusters in the UMAP map, and found that they were expressed and distributed differently in normal, tumor, and metastatic tissues (Figure 2C–2E). Based on cell annotation strategies developed by previous studies, we identified 7 cell types within the GC microenvironment: myeloid cell, epithelial cell, T cell, stromal cell, proliferative cell, B cell, and NK cell, the specific genes expressed by each cell type were shown in the form of bubble maps (Figure 3A), we also labeled the cell clusters on the UMAP map (Figure 3B). Figure 3C depicts the cellular distribution characteristics within normal tissue, tumor tissue, and metastasis tissue.

Figure 2. Single-cell dimensionality reduction. (A) Use the “clustree” package to visualize the division relationships between cell populations under different resolutions. (B) Annotate the cell characteristics of 34 clusters. (C–E) Distribution of 34 kinds of cell clusters in the perspective of UMAP dimensionality reduction algorithm.

Figure 3. Single-cell annotation and parthanatos scores prediction. (A) Expression of specific marker genes in 7 cell types. (B, C) Dimension reduction and annotation of cell clusters based on UMAP algorithm. (D) 5 algorithms scored the gene sets of seven cell clusters. (E) Comparison of the expression difference of 14 model genes in 7 cell types of normal tissue, tumor tissue and metastatic tissue. (F) Comparison of gene set scores of 7 cell clusters among normal tissue, tumor tissue and metastatic tissue.

Parthanatos gene set scoring based on single-cell and spatial transcriptome data

Cells do not rely solely on one or a few genes to carry out their functions. Many genes in the upstream and downstream pathways of the functional pathways vary in expression as the functions vary in strength. Therefore, we used five algorithms (AUCell, Ucells, Singscore, ssgsea, and addmodulescore) to score the parthanatos gene set in single-cell data. The Scoring score is the total score of the five algorithms. Figure 3D shows the expression of parthanatos-associated genes in 7 cell clusters for each algorithm. Myeloid cell, epithelial cell and proliferative cell show strong signals, while NK cells and T cells show weaker signals (Figure 3D).

To compare the differences in parthanatos-related gene expression among normal tissue, tumor tissue and metastasis tissue, we displayed the gene set scores for each cell cluster under the five algorithms in Figure 3F. The different algorithm scores indicate that compared to normal tissue, almost all cell types within the tumor tissue and metastasis tissue show higher parthanatos gene scores. Notably, the gene set scores for epithelial cell, T cell and proliferative cell in tumor tissue and metastasis tissue are significantly higher (Figure 3F). Furthermore, using the UMAP algorithm, we mapped the 7 cell clusters and their respective parthanatos gene set scores onto the merged sample tissue. By comparing the Scoring scores, we can clearly observe that epithelial cell, T cell and proliferative cell in tumor tissue and metastasis tissue have higher parthanatos gene set scores, which is consistent with our previous observations (Figure 4A). In addition, the control data from mouse GC tissue revealed, for the first time, the spatial distribution of parthanatos signals at the tissue level (Figure 4B).

Figure 4. Distribution of parthanatos scores in single-cell atlas and spatial resolution. (A) Based on the UMAP algorithm, the gene set scores of 7 cell clusters were displayed in the combined samples (including normal tissue, tumor tissue and metastatic tissue). (B) Spatial transcriptome data of GC, gene set scoring under 6 gene set scoring algorithms.

Cluster analysis based on parthanatos score

We developed a GC classification model based on parthanatos-related genes, which can clearly distinguish GC patients with different parthanatos features. Using unsupervised clustering analysis algorithm [43], we divided a total of 743 GC patient samples from TCGA and GEO databases into two subtypes (C1 and C2). The results of consistency matrix heatmap, cumulative distribution curve, and delta area curve all confirmed that k=2 is the optimal clustering number (Figure 5A–5C). Violin plots showed the enrichment scores of the two subtypes, where C2 subtype had higher parthanatos scores, indicating higher activity of parthanatos-related genes, while C1 subtype had the opposite trend (Figure 5D). The heatmap displayed the expression patterns of different parthanatos-related genes between C1 and C2 subtypes. Except for some genes, most genes, such as RIPK1, DDB1, CUL4A, AIMP2, and PARP1, were expressed at higher levels in the C2 subtype compared to the C1 subtype, further confirming the higher expression activity of parthanatos-related genes in C2 subtype patients (Figure 5F). We also studied the prognosis differences between the two subtypes using the “survival” package and “survminer” program in R Studio. It was found that patients with the C2 subtype had better prognosis and higher survival rates, suggesting that a higher parthanatos score is associated with better prognosis in GC patients (Figure 5E).

Figure 5. Parthanatos-based clustering analysis. (A) Define a consensus matrix heat map of k = 2 clusters and their associated regions. (B) Cumulative distribution function (CDF) curve under different cluster number. (C) The relative change in area under the CDF curve for different values of k. (D) The violin plot shows enrichment scores for two clusters (C1 and C2), p<0.001. (E) Survival curves of C1 and C2 clusters, purple for C1, green for C2. (F) Expression difference of parthanatos-related genes between CI and C2 subtypes. (G) Differences in the activity of metabolic pathways between CI and C2 subtypes. (H) Differences in immune pathway activity between CI and C2 subtypes.

Metabolic reprogramming and immune microenvironment characteristics of the parthanatos scoring model

In order to fully explore the intrinsic molecular features of GC patients in different subtypes, we collected 42 classic metabolic pathways and 24 classic immune pathways based on the KEGG database. We evaluated the metabolic and immune signaling intensities in 743 GC samples and depicted the distribution of each metabolic and immune signal in the form of heat maps. In the C1 subtype, signals of metabolic pathways such as retinol metabolism, drug metabolism cytochrome P450 metabolism, and metabolism of xenobiotics by cytochrome P450 metabolism were enhanced, while signals of metabolic pathways such as purine metabolism, pyrimidine metabolism, cysteine and methionine metabolism, and riboflavin metabolism were weakened. The opposite was observed in the C2 subtype (Figure 5G). Regarding immune-related pathways, signals of immune pathways such as TGF beta signaling pathway, chemokine signaling pathway, and intestinal network for IGA production were enhanced in the C1 subtype, while signals of immune pathways such as P53 signaling pathway, proteasome, progesterone mediated oocyte maturation, oocyte meiosis, RNA degradation, spliceosome, and nucleotide excision repair were enhanced in the C2 subtype (Figure 5F).

In addition, we used the ESTIMATE algorithm to analyze the proportions of immune cells, stromal cells, and tumor cells in tumor tissues, and obtained the immune score, stromal score, ESTIMATE score, and tumor purity for both C1 and C2 subtypes. From the graphs, it can be seen that C1 subtype had more stromal and immune cells, while C2 subtype had higher tumor purity (Figure 6A–6D). This phenomenon may be caused by immune function suppression and abnormal accumulation of immune cells in the tumor microenvironment of the C1 subtype. To further analyze the differences in the immune microenvironment between the two subtypes, we used seven immune cell infiltration prediction algorithms provided by the TIMER2.0 online platform (including TIMER, CIBERSOFT, CIBERSOFT-ABS, QUANTISEQ, XCELL, EPIC, and MCPCOUNTER) to analyze the extent of immune cell infiltration between the C1 and C2 subtypes. From the graphs, it can be observed that regardless of the algorithm used, the number of B cells in the C1 subtype was significantly higher than that in the C2 subtype, while the opposite was observed for neutrophils. Other immune cells also showed certain differences between the two subtypes (Figure 6E). Finally, Figure 6F shows the expression patterns of immune checkpoint-related genes in the C1 and C2 subtypes. Genes such as CD8A, CD44, NRP1, TNFSF14, TNFSF15, CD27, and CD48 were significantly higher expressed in the C1 subtype (Figure 6F). This result suggests an excessive activity of immune checkpoints in the tumor microenvironment of the C1 subtype, which may lead to tumor immune escape to a certain extent and subsequently cause excessive accumulation of immune cells. This is consistent with our previous findings on the abnormal accumulation of immune cells in the C1 subtype and the differences in immune microenvironment characteristics between the two subtypes.

Figure 6. Analysis of immune microenvironment between different subtypes. (A–D) The difference between the two clusters in immune score, stromal score, ESTIMATE score and tumor purity, P-values are represented by *. (E) Seven algorithms were used to analyze immune cell or functional differences between CI and C2 subtypes. (F) Differences in the expression of immune checkpoint related genes between CI and C2 subtypes, P-values are indicated by *. *: p<0.05, **: p<0.01, ***: p<0.001, “ns” indicates no significant difference.

Significance of parthanatos scoring model for GC-targeted drug therapy

To further explore the potential value of parthanatos in the clinical treatment of GC patients, we used the “oncoPredict” package and GDSC2 dataset to predict the response of 743 GC patients to various targeted therapies. We analyzed the potential beneficial drugs for GC patients in the C1 and C2 subtypes and plotted box plots to determine the impact of parthanatos on the sensitivity to 12 common targeted drugs in GC. Patients in the C1 subtype showed higher sensitivity to dasatinib, while patients in the C2 subtype showed higher sensitivity to 11 other drugs such as carmustine, sorafenib, rapamycin, and sepantronium (Supplementary Figure 3A–3L). This result may provide precise guidance for future individualized targeted drug therapy for GC patients.

Relationship between parthanatos and immune cell infiltration in GC tissues

The characteristics of GC include a high degree of heterogeneity in tumor cells and tumor immune microenvironment. To investigate this, we evaluated the immune cells and immune function scores of 743 GC samples based on the ssGSEA algorithm. We explored the potential associations between parthanatos-related genes and 29 types of immune cells or functions using Spearman correlation analysis. We also calculated the correlation between parthanatos scores and immune cell infiltration. We found that different parthanatos-related genes have varying degrees of correlation with various immune cells or functions. Among them, RAB33A, ESR1, and ESR2 were positively correlated with most immune cells or functions, while TOMM20, DDB1, and CUL4A were negatively correlated (Figure 7A). The bubble plot displays the correlation between classic immune infiltration-related cells and parthanatos, including both positive and negative correlations (Figure 7B). We then selected the three immune cells with the highest correlation (Treg, Th2 cell, Th1 cell) and performed correlation analysis with parthanatos scores. The R-values were 0.38, 0.28, and 0.27, respectively, indicating positive correlations (Figure 7C–7E). This finding is consistent with our previous results, where we observed a higher quantity of Th1 and Th2 cells in the C2 subtype compared to the C1 subtype using the XCELL algorithm.

Figure 7. Relationship between parthanatos score and immune microenvironment in GC. (A) Heat map depicting the relationship between parthanatos-related genes and levels of immune cell infiltration. Positive correlations are shown in red, negative correlations are shown in blue, and p-values are shown by *. (B) The bubble plot shows the correlation between parthanatos score and the levels of infiltration of various immune cells, with the size of the circle indicating the size of the correlation and the color of the circle indicating statistical significance. (C–E) Correlation between parthanatos score and Infiltration levels of three kinds of immune cells. *: p<0.05, **: p<0.01, ***: p<0.001, “ns” indicates no significant difference.

Construction and verification of survival and prognosis model of GC patients based on parthanatos-related genes

Based on the limma and “clusterProfiler” packages, we screened for differentially expressed genes between different parthanatos subtypes for GO enrichment analysis and KEGG enrichment analysis (Supplementary Figure 4). These genes were then used to construct a prognostic model for GC patients (SGCA, JAM2, SHISA3, DES, PDK4, SFRP2, GRP, TNC, PAEP, FBLN5, GLDC, CCDC80, HAND2, and PPP1R14A), Figure 3E shows the expression of 14 model genes in seven cell types of normal tissue, tumor tissue, and metastatic tissue (Figure 3E). In the train cohort, we performed univariate Cox regression, LASSO regression, and multivariate Cox regression analysis on the differentially expressed genes to build the prognostic model (Supplementary Figure 5). We divided the GC patients in the train cohort into high- and low-risk groups based on the median risk value. We then compared the prognosis differences and analyzed the ROC curve to evaluate the predictive accuracy of the model (Figure 8A and Supplementary Figure 6). Subsequently, in the internal validation set 1 (test1 cohort), internal validation set 2 (test2 cohort), and external validation set (test3 cohort), we performed multivariate Cox regression using the model genes selected from the train cohort. We used the predict function to predict the risk values for each sample and used the median risk value from the train cohort as the cut-off value to divide the three validation sets into high- and low-risk groups. We then conducted prognosis difference comparisons and analyzed the ROC curve. In different cohorts, the high-risk group consistently showed worse prognosis compared to the low-risk group, and the AUC values for 3-year and 5-year ranged from 0.6 to 0.8, indicating the stability and potential generalizability of the prognostic model (Figure 8B–8D and Supplementary Figure 6). In order to explore the differences in the immune microenvironment between the high and low-risk groups, we used the TIMER2 platform to predict the immune cell infiltration abundance using seven immune algorithms in the four cohorts. We also compared the expression differences of immune checkpoint genes between different risk groups (Supplementary Figure 7). The results showed that in all algorithms and cohorts, the low-risk group had higher numbers of plasma cells, NK cells, and certain CD4+ T cells (such as activated memory CD4+ T cells, Th1 cells, and Th2 cells), while tumor-associated fibroblasts (CAFs) showed the opposite trend. Additionally, the box plot results showed higher expression of immune checkpoint-related genes in the high-risk group compared to the low-risk group. These differences in the immune microenvironment may be one of the main intrinsic factors leading to differences in prognosis and clinical data between the high- and low-risk groups.

Figure 8. Survival curves and ROC curves for the high- and low-risk groups in the four cohorts. (A–D) Comparison of survival curve and ROC curve between train cohort, test1 cohort, test2 cohort and test3 cohort.

To facilitate the clinical application of the prognostic model, we introduced clinical features of GC as variables in the model and used the “rms” package to construct a column line chart for accurate prediction of GC patient prognosis. The column line chart consists of 10 parallel lines. Each row represents a score, with gender in the second row, risk type in the third row, grade in the fourth row, stage in the fifth row, and age in the sixth row. The total score in the seventh row is obtained by adding the scores of age, grade, stage, and risk type. With this chart, we can easily estimate the survival rates of GC patients at 1, 3, and 5 years (Figure 9A). Furthermore, Calibration plots and ROC curves were used to assess the predictive accuracy of the column line chart. The calibration curves for 1, 3, and 5-year survival closely align with the diagonal line (Figure 9B), and the AUC values were 0.691, 0.683, and 0.709, respectively (Figure 9C). This indicates that both our prognostic model and the column line chart construction are accurate and have predictive value for prognosis.

Figure 9. Establishment of clinical nomogram. (A) A nomogram for predicting prognostic survival time in GC patients. (B) Calibration curves of 1-, 3-, and 5-year survival in GC patients. (C) ROC curves of 1-, 3-, and 5-year survival in GC patients.