Demographic information, clinical data and tissue samples were obtained from 214 patients with GC who had undergone gastrectomy as the primary treatment between 2016 and 2019 at the Harbin Medical University (HMU) Cancer Hospital. These data were used to construct the HMU-GC cohort. All samples were collected after written informed consent was obtained from the patients. This study was approved by the Institutional Review Board of the HMU Cancer Hospital. RNA isolation, library construction and mRNA sequencing were performed by Novogene (Beijing, China). The data were deposited in the GEO repository (GSE184336).

We systematically searched for publicly available GC gene expression datasets in GEO and TCGA databases. Datasets missing the follow-up data were excluded. In addition, to enhance the robustness of downstream analyses, samples with survival time less than 3 months were excluded (Joung and Merkow, 2021; Resio et al., 2021). A total of 12 public treatment-naive GC cohorts (GSE62254/ACRG, GSE15459, GSE57303, GSE34942, GSE38749, GSE15456, GSE84437, GSE26901, GSE26899, GSE13861, GSE26253 and TCGA-STAD) were selected for further analysis. In addition, four GC cell datasets, namely, Cancer Cell Line Encyclopedia (CCLE), GSE22183, GSE15455 and GSE146361, were included as references.

To process microarray data, the raw CEL files obtained from Affymetrix were processed using the robust multichip average (RMA) algorithm for background correction and normalisation using the affy package (Irizarry et al., 2003; Gautier et al., 2004). The raw data from Illumina were processed using the limma package (Ritchie et al., 2015). For microarray datasets without raw data, the normalised matrix files were directly downloaded.

For high-throughput sequencing data obtained from the HMU-GC and TCGA-STAD datasets, raw read count values were converted to transcripts per kilobase million (TPM) values, which are more similar to those generated from microarrays and are comparable between samples (Wagner et al., 2012). Batch effects from non-biological technical biases among different datasets were corrected using the ComBat algorithm in the sva package (Leek et al., 2012).

The gene expression profile at the probe level (or Ensembl ID level) was converted to the official gene symbol level using the biomaRt package (Durinck et al., 2005). When multiple probes (or Ensembl IDs) were mapped to the same gene symbol, the probe (or Ensembl ID) with the largest mean expression values across samples was selected.

Tumour purity was calculated using the ESTIMATE algorithm (Yoshihara et al., 2013). In addition, GC tissues in the TCGA-STAD cohort were reviewed to infer GC purity based on visual evaluation of the whole-slide images of haematoxylin and eosin (H&E) staining. Two pathologists independently confirmed the results of histological purity.

Weighted correlation network analysis (WGCNA) was performed using the WGCNA package to screen for purity-related gene modules (Langfelder and Horvath, 2008). A scale-free topology fitting index (R²) of 0.85 was set as the threshold to construct a signed weighted gene co-expression network. The minimum co-expression module size was set as 30, and the minimum cut height for module merging was set as 0.25. A biweight midcorrelation coefficient (R) > 0.4 and p-value < 0.05 were set as thresholds to screen for gene modules significantly associated with GC purity.

Nearest template prediction (NTP) analysis was performed to classify samples based on the known clinical and molecular features (Hoshida, 2010). Univariate and multivariate Cox regression analyses were performed to calculate the hazard ratio (HR) and 95% confidence interval (CI). The Kaplan–Meier survival analysis with log-rank test was used to compare the survival rates of patients in different subgroups. Differences in prognosis between the high- and low-purity subgroups were compared using the restricted mean survival time (RMST) analysis (Kim et al., 2017). In addition, subgroup analyses were performed to examine the association between GC purity and other clinical characteristics.

Furthermore, GC purity was integrated with other independent prognostic factors to generate a composite prognostic nomogram for model visualisation and evaluation of clinical applications. The predictive value of the nomogram was compared with that of the TNM staging system in terms of concordance index (C-index). The performance of the composite model was evaluated based on calibration curves, time-dependent receiver operating characteristic (ROC) analysis, and decision curve analysis (DCA) (Vickers and Elkin, 2006).

Gene annotation enrichment analysis was performed using the clusterProfiler package (Yu et al., 2012). Gene set enrichment analysis (GSEA) was performed to screen for biological processes related to GC purity (Subramanian et al., 2005), and gene set variation analysis (GSVA) was performed to quantify the pathway enrichment scores using the GSVA package (Hänzelmann et al., 2013) and screen for significantly enriched pathways in each cluster. The well-defined ‘Hallmark gene sets’ were selected to quantify the pathway activity (Liberzon et al., 2015).

To quantify the infiltration level of stromal and immune cells in each GC sample, we calculated the stromal and immune scores using the ESTIMATE algorithm (Yoshihara et al., 2013). Based on the guidelines of transcriptome-based cell-type quantification methods (Sturm et al., 2019), we used the MCPcounter and xCell algorithms to quantify specific immune and stromal cells in the GC samples (Becht et al., 2016; Aran et al., 2017). Based on the whole-slide images of H&E staining obtained from the TCGA-STAD cohort, we further used the convolutional neural network, a supervised deep-learning approach, to identify the proportion of tumour-infiltrating lymphocytes (TILs) in digitised H&E-stained tissue specimens (Saltz et al., 2018).

Seven common chemotherapeutic agents (5-fluorouracil, cisplatin, oxaliplatin, capecitabine, paclitaxel, docetaxel and irinotecan), which are approved for GC treatment, were selected to predict the chemotherapeutic response. Based on two public drug sensitivity databases, namely, CTRP (Basu et al., 2013) and PRISM (Corsello et al., 2020), the chemotherapeutic response was predicted via the pRRophetic package using ridge regression to estimate the area under the curve (AUC) value of each sample (Geeleher et al., 2014). Lower AUC values indicated increased sensitivity to treatment. The prediction accuracy was evaluated via 10-fold cross-validation based on each training set. Default values were selected for all parameters, including the combat algorithm for removing batch effects and the mean value for summarising duplicate gene expression.

Differences between groups for continuous variables were evaluated using the Kolmogorov–Smirnov, Mann–Whitney or Kruskal–Wallis tests. The two-sided Pearson’s chi-squared or Fisher exact test was used to analyse the categorical data. The association between continuous variables was tested using Spearman correlation analysis. Restricted cubic splines were used to test potential non-linear associations. All statistical analyses were conducted using the R software, and p-values were two-sided. A p-value < 0.05 was considered statistically significant. The Benjamini–Hochberg method was used to control the false discovery rate (FDR) for multiple hypothesis testing.

More than 40 GC cell lines involving 215 cell samples were analysed as references (Supplementary Table S1). A high degree of purity, i.e. median purity of 0.999 (inter-quartile range [IQR], 0.996–1), was observed in GC cell lines. Furthermore, no significant difference was found between the microarray and RNA-seq datasets (Kolmogorov–Smirnov test, p = 0.789), which verified the validity and robustness of the ESTIMATE algorithm.

In the GC meta-dataset (Supplementary Figure S1A), a total of 2,599 GC tissues from 13 cohorts were initially selected. After excluding samples with survival time less than 3 months (Joung and Merkow, 2021; Resio et al., 2021), 2,259 samples were retained for further analyses. As demonstrated in Figure 1A, GC samples were arranged in order of increasing purity, and a wide range (0.215–0.964) of purity was observed among the 2,259 samples, with a median purity of 0.681 (IQR, 0.608–0.759). In addition, the morphological features of GC cells were assessed in 179 patients in the TCGA cohort (Figure 1B; Supplementary Table S2). Consistent with the findings of the previous study (Aran et al., 2015), GC purity evaluated based on the ESTIMATE algorithm was significantly correlated with morphology (Figure 1C). However, significant differences were observed between the two evaluation methods (Figures 1D,E), which suggested that morphology evaluation can only provide a qualitative estimation of GC purity.

Furthermore, we investigated different samples from the same patient to determine their consistency. In the GSE14209 dataset, 22 patients were analysed twice (pre- and post-chemotherapy). We observed high concordance between the pre- and post-chemotherapy samples (Supplementary Figure S1B). The results suggest that tumour purity levels among patients with cancer are robust and consistent, implying that tumour purity may be an intrinsic property of GC.

Based on the speculation that tumour purity is an intrinsic property of GC, we assessed the association between GC purity and clinical and molecular features. We found a negative correlation between GC purity and age at diagnosis (Spearman R = 0.122; p < 0.001; Supplementary Figure S1C). According to the TNM staging system and histopathological features, GC samples were classified into the corresponding clinicopathological and molecular features. A consistent decrease in GC purity was observed as the TNM stages advanced (Figure 2A). Based on the Lauren classification, intestinal-type GC had the highest purity level, whereas diffuse-type GC was accompanied had the lowest purity level (Figure 2B). These results suggest that low tumour purity is closely associated with the malignant progression of GC.

Based on clinical features, low-purity GC samples were more likely to belong to the mesenchymal subtype, whereas high-purity GC samples were enriched in the epithelial subtype (Figure 2C). Based on the other three molecular classifications, namely, the ACRG (Figure 2D), Singapore (Figure 2E) and TCGA (Figure 2F) classification systems, decreased tumour purity levels were associated with epithelial–mesenchymal transition (EMT) and invasive or genomically stable (GS) molecular subtypes, which were usually correlated with the invasion and metastasis of GC. These findings indicate that GC purity is closely related to specific clinicopathological and molecular features.

We found that low purity was significantly related to the malignant phenotype. To further evaluate the prognostic value of GC purity in combination with other clinical variables, we performed univariate and multivariate Cox proportional hazard regression analyses based on combined clinicopathologic variables in the GC meta-dataset cohort. In addition to chemotherapy history and TNM stage, which is a well-known prognostic factor, GC purity was identified as a significant predictor of overall survival (OS) and recurrence-free survival (RFS) based on the univariate analysis (Table 1).

Although the Epstein–Barr virus (EBV) subtype is characterised by lower purity scores (Figure 2F), it is associated with a good prognosis (Sohn et al., 2017). Therefore, the relationship between GC purity and prognosis may not necessarily be linear. We used restricted cubic splines to assess the possible non-linear association between GC purity and prognosis. However, a non-linear relationship (a U-shaped curve) was observed between GC purity and OS (Figure 3A) and RFS (Figure 3B).

Furthermore, patients with GC were stratified into two (high or low purity) and three (high, moderate or low purity) subgroups using the optimal cut-off values determined using the X-tile software (Camp et al., 2004). In both subgroup analyses, the OS (Figures 3C,D) and RFS (Figures 3E,F) rates were lower in patients with low purity scores than in patients with high purity scores. The RMST difference also determined a benefit of high purity, and the benefit seemed to increase over time (Table 2). For example, the RMST differences between the two groups were 2 months for OS and 4 months for RFS after 5 years of follow-up, which increased to 7 months for OS and 12 months for RFS after 10 years. These results suggested that low GC purity was associated with a shorter survival time and faster recurrence.

When all relevant clinical variables were included in the multivariate Cox regression analysis, GC purity was identified as a significant prognostic factor (Table 1). Subgroup analyses were further performed to assess the association between GC purity and other prognostic factors. No significant association was observed between GC purity and both OS and RFS (Table 3), indicating that GC purity retained its prognostic relevance even after classic clinicopathological prognostic features were considered.

Furthermore, we integrated these independent factors (tumour purity, TNM stage and chemotherapy history) to establish a prognostic nomogram for predicting the OS of patients with GC (Supplemenatry Figure S2A). The calibration plot revealed an optimal agreement between nomogram prediction and actual observation at different time points (Figure 3G). Moreover, we compared the OS prediction power between the nomogram and clinical model. The C-index of the nomogram was significantly higher than that of the clinical model for predicting OS (Figure 3H). In addition, the time-dependent ROC curve confirmed that the nomogram exhibited better performance in predicting the prognosis of GC (Figure 3I).

Lastly, DCA curves were used to assess the clinical usefulness and net benefit of these two models. The composite nomogram demonstrated a larger net benefit than that exhibited by the clinical model within most of the threshold probabilities (Figure 3J). Similar results were also found for RFS (Supplemenatry Figure S2B; Figures 3K–N). These results indicated that integrating GC purity into a prognostic model can significantly improve its predictive power. In addition, the nomogram exhibited better clinical utility for predicting prognosis, emphasising the requirement of considering GC purity in clinical practice.

Given that tumour purity is closely related to the clinicopathological, molecular and prognostic features, we used RNA expression profiles to examine the underlying biological processes associated with GC purity. We performed principal components analysis (PCA) to identify the transcriptomic features related to GC purity and found a strong association between mRNA expression profiles and GC purity (Figure 4A), implying that distinct biological phenotypes are attributed to varied GC cell percentages.

Furthermore, we performed GSEA to assess the biological features associated with GC purity. The results suggested that samples with different GC purity exhibited distinct biological processes (Figure 4B). Consistent with the clinical and molecular features, low tumour purity was significantly related to pathways associated with invasion and metastasis and multiple immune-related pathways. However, high purity was considerably associated with metabolism and proliferation-related pathways.

Subsequently, we used WGCNA to obtain purity-related modules. We constructed a cluster dendrogram according to the soft threshold power (Supplemenatry Figure S3A) and identified 11 colour modules (Supplemenatry Figure S3B). Among the identified modules, four were highly associated (|R| > 0.4) with tumour purity (Figure 4C; Supplemenatry Table S3). In addition, gene significance significantly correlated with module membership in each module (Supplemenatry Figure S3C), suggesting that genes in these modules might play an essential role associated with GC purity.

Furthermore, module enrichment analyses were performed to explore the biological features of purity-related modules (Figure 4D). Consistent with the results of GSEA, genes in the red module were enriched in pathways associated with metabolic activation (Supplemenatry Table S4). Genes in the blue module were significantly enriched in proliferation-specific pathways (Supplemenatry Table S5). In addition, genes in the turquoise (Supplemenatry Table S6) and yellow (Supplemenatry Table S7) modules were prominently related to stromal and immune activation pathways, respectively.

Eventually, we analysed proteomic data generated using reverse-phase protein arrays in the TCGA cohort (Figure 4E). Consistent with the clinical and gene expression data, α-catenin, E-cadherin and other crucial epithelial adhesion proteins in epithelial cells were positively correlated with tumour purity. Low tumour purity was associated with high expression of MYH11, RICTOR and CAV1, which are markers for mesenchymal lineage or EMT. In addition, low GC purity was associated with increased expression of numerous stromal- and immune-associated proteins. High GC purity was, however, associated with increased expression of multiple proliferation- and metabolism-associated proteins.

Low GC purity was markedly associated with stromal and immune activation pathways (Figure 4). Therefore, we further examined the relationship between TME features and GC purity to characterise TME heterogeneity. Both stromal and immune scores, representing the overall infiltration of stromal and immune cells, respectively, in tumour tissues, were inversely correlated with GC purity (Figure 5A), suggesting that TME with decreased tumour purity had a significantly increased infiltration of stromal and immune cells. In addition, TILs were confirmed to be negatively correlated with GC purity (Figures 5B,C) (Saltz et al., 2018).

Furthermore, we used the MCPcounter algorithm to evaluate the relative abundance of infiltrating stromal and immune cell subpopulations against GC purity (Figure 5A). T cells, cytotoxic lymphocytes, B lineage, myeloid dendritic cells, monocytic lineage, fibroblasts and endothelial cells exhibited a consistently negative correlation (R < −0.4) with GC purity. We further characterised TME heterogeneity using the xCell algorithm, which can evaluate as many as 39 different immuno-oncology cell types (Figure 5A). Similar to the previous result, fibroblasts, endothelial cells, myeloid dendritic cells and several T and B cell types were found to be negatively correlated with GC purity (R < −0.4). Although numerous immune cell types (CD8⁺ naive T cells, Th₁ cells and Th₂ cells) were positively associated with GC purity, the correlation was relatively weak.

Lastly, we attempted to assess the prognostic implications of these purity-related stromal and immune cells in GC (Figure 5D). Survival analyses revealed that fibroblasts and endothelial cells were consistently associated with worse OS and RFS. However, T cells, particularly memory T cells, were consistently associated with better OS and RFS. The results of the xCell algorithm revealed that monocytes might correlate with worse OS and RFS. Therefore, T cells, fibroblasts, endothelial cells and monocytes can be a cluster of nontumour cells contributing to TME and specific clinical outcomes for patients with GC with varying purity levels.

Chemotherapy is the mainstay of treatment for patients with advanced GC. However, chemotherapy resistance is the primary cause of treatment failure. Given that distinct clinical, biological and microenvironmental features are associated with GC purity, we further investigated the relationship between GC purity and chemotherapy response to promote individualised treatment decisions. The analysis was performed in a subset of patients in the meta-data cohort with data on chemotherapy (n = 863). Of the 863 patients, 732 received adjuvant chemotherapy. We found that GC purity could serve as an independent predictor for adjuvant chemotherapy benefit (Table 4). Patients with GC with a high purity level who received adjuvant chemotherapy had a significant OS benefit (HR, 0.101; 95% CI, 0.014–0.719; p = 0.022).

Furthermore, we determined whether varied purity levels were associated with differences in clinical benefit from adjuvant chemotherapy. Among patients who received adjuvant chemotherapy (n = 732), we found significant OS differences among three subgroups (Figures 6A,B). Adjuvant chemotherapy was associated with significantly increased OS and RFS rates in patients with GC in the high- and middle-purity subgroups (Figures 6C,D). However, no benefit from adjuvant chemotherapy was observed among patients with GC in the low-purity subset (Figure 6E).

Eventually, we assessed the value of GC purity in facilitating individualised chemotherapy regimens. We performed ridge regression to predict the drug susceptibility results for each sample (Figure 6F) using the CTRP and PRISM-derived drug response data. Spearman correlation analysis between the AUC value and GC purity was used to select agents with a significant correlation coefficient (Figure 6G). In addition, a differential analysis among the three subtypes was conducted to identify agents with lower estimated AUC values in each subgroup (Figures 6H,I). These analyses yielded two CTRP-derived compounds (fluorouracil and paclitaxel) and four PRISM-derived compounds (fluorouracil, capecitabine, cisplatin, and oxaliplatin). These compounds had lower estimated AUC values in the high-purity subgroup and a negative correlation with GC purity, suggesting that chemotherapy resistance to these five agents was observed in patients with GC with low tumour purity.