Cell differentiation trajectory predicts patient potential immunotherapy response and prognosis in gastric cancer

Quality control and normalization of scRNA-seq data

In this study, 402 cells from 6 GC samples were obtained from GSE112302. After quality control and normalization, 2 nonconforming cells were excluded, and 400 cells were screened for further analysis (Figure 1A). No correlation between sequencing depth and mitochondrial gene sequences was detected (Figure 1B). There was a significant positive correlation between sequencing depth and total intracellular sequences (R = 0.38, Figure 1C). A total of 16,288 genes were analyzed, of which 14,788 had low intercellular variation and 1,500 had high variation (Figure 1D).

Figure 1. Quality control and normalization of scRNA-seq data, dimensionality reduction and cell trajectory analysis. (A) After quality control and normalization, 2 nonconforming cells were excluded, and 400 cells were screened for further analysis. (B) Correlation analysis between sequencing depth and mitochondrial gene sequences was detected. (C) Correlation analysis between sequencing depth and total intracellular sequences. (D) A total of 16,288 genes were analyzed, of which 14,788 had low intercellular variation and 1,500 had high variation. (E) PCA based on scRNA-seq data. (F) Eighteen PCs with significant differences were identified with P < 0.5. (G) Four hundred GC cells were aggregated into 8 clusters and the top 10% of marker genes in each cluster are displayed on the heat map. (H, I) Pseudotime and trajectory analysis. PCA: principal component analysis, PCs: principal components, GC: gastric cancer.

Cell trajectory analysis identified three GC subsets

Principal component analysis (PCA) was used for preliminary dimensionality reduction of scRNA-seq data. The results showed that there was no significant segregation among GC cells (Figure 1E), so the first 18 principal components (PCs) with significant differences were selected for further analysis (Figure 1F). According to the t-distributed stochastic neighbor embedding (tSNE) algorithm, 400 GC cells were aggregated into 8 clusters, a total of 3339 marker genes were identified by differential analysis, and the top 10% of marker genes in each cluster are displayed on the heat map (Figure 1G). Eight clusters were annotated based on marker genes: clusters 0, 1, 2, 3, 4, 6, and 7 were all cancer cells, and cluster 5 correlated with macrophages. Pseudotime and trajectory analysis indicated that cluster 1 was distributed in subset I, containing cancer cells; clusters 0/2/3 were in subset II, consisting of cancer cells; clusters 5/7 were distributed in subset III, composed of cancer cells and macrophages; and clusters 4/6 were scattered over the three subsets and were also cancer cells (Figure 1H, 1I).

Molecular functional analysis of three subsets based on GDRGs

The study identified 402 GC differentiation-related genes (GDRGs) in subset I, 1443 GDRGs in subset II and 1016 GDRGs in subset III. Based on Gene Ontology (GO) biological process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, the GDRGs in subsets I/II were involved in mRNA catabolic processes and energy metabolic disorders (Figure 2A–2D), those in subset II were also associated with hypoxia tolerance (Figure 2C, 2D), and those in subset III were closely related to inflammation and immune-related pathways (Figure 2E, 2F).

Figure 2. Functional analysis of three subsets based on GDRGs. (A) GO BP for subset I genes. (B) KEGG enrichment analysis for subset I genes. (C) GO BP for subset II genes. (D) KEGG enrichment analysis for subset II genes. (E) GO BP of subset III genes. (F) KEGG enrichment analysis of subset III genes. GDRG: gastric cancer differentiation-related gene, GO: Gene Ontology, BP: biological process, KEGG: Kyoto Encyclopedia of Genes and Genomes.

Three GDRG-based molecular subtypes of GC patients from the GSE84437 dataset

GDRG-based consensus clustering analysis was completed in GSE84437, and three molecular subtypes that contained all the GC samples were identified at a clustering threshold of maxK = 9 (Figure 3A–3C). The Kaplan-Meier analysis determined the statistical significance of the consensus clustering results for GC, with subtype I (C1) having the best overall survival (OS), subtype II (C2) having the second best OS, and subtype III (C3) having the worst OS (Figure 3D, P = 0.004). From subtype I (C1) to subtype III (C3), patient age decreased, and patients tend to have advanced T stage and N stage, but there was no significant difference in the sex distribution among the three subtypes (Figure 3E and Supplementary Figure 1). In addition, the up/downregulated GDRGs in subsets I/II/III showed the same expression trend in subtypes I/II/III (C 1/2/3) (Figure 3F–3H and Supplementary Figure 1), which indicated that subtypes I/II/III (C 1/2/3) were separately composed of subsets I/II/III.

Figure 3. GDRG-based consensus clustering analysis of GC patients from the GSE84437 dataset. (A–C) Three molecular subtypes were identified at a clustering threshold of K = 9. (D) Kaplan-Meier analysis between the three molecular subtypes. (E) Proportion of clinicopathologic features among the three molecular subtypes. (F–H) The up/downregulated GDRGs in subsets I/II/III showed the same expression trends as subtypes I/II/III (C1/2/3). GC: gastric cancer, GDRG: GC differentiation-related gene.

Comprehensive analysis of tumor microenvironment scores and immune cell infiltration across three molecular subtypes

Based on tumor microenvironment scores, the current study found that immune/stromal scores increased in turn in subtypes I/II/III (C 1/2/3) (Figure 4A, 4B, P < 0.05), while tumor purity gradually decreased in subtypes I/II/III (C 1/2/3) (Figure 4C, P < 0.001). The content of 22 immune cells in each sample was calculated according to the CIBERSORT algorithm and displayed visually such that different colors represent different cell types (Figure 4D). Differential analysis demonstrated that subtype I (C1) contained more memory-activated CD4+ T cells, M0 macrophages, and activated NK cells (Figure 4E, all P < 0.05), which was associated with better OS (Figure 4F). Clusters 2/3 had more memory resting CD4+ T cells, resting mast cells and naive B cells (Figure 4E, all P < 0.05), and the higher infiltration density of these cells was correlated with worse OS (Figure 4G). In addition, the infiltration density of the remaining immune cells (e.g., neutrophils) also showed statistically significant differences in the three subtypes (Figure 4E, all P < 0.05). However, due to the small numbers of these cells in the tumor microenvironment, they were not considered in this study.

Figure 4. Comprehensive analysis of tumor microenvironment scores and immune cell infiltration across three molecular subtypes. (A–C) Tumor microenvironment scores across three molecular subtypes. (D) The contents of 22 immune cells in each sample from the GES84437 dataset. (E) Differential analysis showed that the contents of 14 kinds of immune cells were different in the three subtypes. (F) Kaplan-Meier analysis of memory-activated CD4⁺ T cells, activated NK cells (P < 0.05) and M0 macrophages. (G) Kaplan-Meier analysis of resting mast cells, memory resting CD4⁺ T cells and naive B cells.

Expression levels of 38 immune checkpoint genes (ICGs) across three molecular subtypes and prognostic analysis

Thirty-eight confirmed ICGs were obtained from previous studies [7–15]. Through differential expression analysis, we observed significantly high expression of 10 ICGs (namely, CD80, CTLA4, IFNG, IL23A, LDHA, PD-L1, PVR, TNFRSF9, TNFSF9 and YTHDF1) in subtype I (C1) (Figure 5A, all P < 0.05). Importantly, upregulated CD80 (P < 0.05), CTLA4 (P < 0.05), TNFRSF9 (P < 0.05), IFNG, LDHA and TNFSF9 correlated with better OS (Figure 5B). Five ICGs (namely, FGL1, JAK1, LAMA3, LGALS9 and VTCN1) were upregulated in subtype II (C2) (Figure 5A, all P < 0.05), and Kaplan-Meier analysis indicated that highly expressed VTCN1 (P < 0.05), JAK1 and LAMA3 predicted worse OS (Figure 5C). Nine ICGs (namely, B2M, CD28, CD40LG, CD8A, HAVCR2, LDHB, PTPRC, TNFSF18 and TNFSF4) were significantly overexpressed in subtype III (C3) (Figure 5A, all P < 0.05), and Kaplan-Meier analysis showed that LDHB and TNFSF4 expression was associated with poor OS, while PTPRC predicted better OS (Figure 5D). These results not only provide a reasonable molecular basis for the different OS in the three subtypes but also may guide immunotherapy.

Figure 5. Expression levels of 38 ICGs across three molecular subtypes and prognostic analysis. (A) Differential expression analysis of 38 ICGs. (B) Kaplan-Meier analysis of CD80, CTLA4, TNFRSF9, IFNG, LDHA and TNFSF9. (C) Kaplan-Meier analysis of VTCN1, JAK1 and LAMA3. (D) Kaplan-Meier analysis of LDHB, TNFSF4 and PTPRC. ICGs: immune checkpoint genes.

Generation, evaluation and validation of a prognostic risk scoring signature

After the intersection of GDRGs in The Cancer Genome Atlas (TCGA) and GSE84437 cohorts, 1729 GDRGs were enrolled in weighted correlation network analysis (WGCNA). Nine modules were accessed with soft threshold = 8 (Figure 6A–6D), of which 6 modules (namely, blue module, yellow module, pink module, red module, turquoise module and gray module) were closely related to GC grade (Figure 6E). A total of 258 differentially expressed GDRGs (Figure 6F and Supplementary Figure 2) were obtained from 6 modules, and 33 prognostic GDRGs (Figure 6G) were further screened by univariate analysis and incorporated into multivariate Cox regression analysis. Finally, a prognostic risk scoring (RS) signature consisting of 8 GDRGs was generated. The RS of each sample could be calculated according to the relative coefficient and gene expression. RS = (0.00048 * expression of VCAN) + (-0.00096 * expression of TNFAIP2) + (0.00125 * expression of STMN2) + (0.00027 * expression of RNASE1) + (0.00013 * expression of DUSP1) + (0.00082 * expression of AQP2) + (0.00205 * expression of ADAM8) + (0.00006 * expression of TFF1). According to the prognostic risk scoring signature, the RS of each GC sample in the TCGA and GSE84437 cohorts was calculated, and the OS of the low-risk group was significantly better than that of the high-risk group (Figure 7A, TCGA: P = 2.681e-05; Figure 7B, GSE84437: P = 1.282e-03). In addition, the areas under the receiver operating characteristic (ROC) curves for predicting 1-year, 3-year and 5-year OS were 0.695, 0.666 and 0.701, respectively, in the TCGA cohort (Figure 7C), and 0.605, 0.593 and 0.604, respectively, in the GSE94437 dataset (Figure 7D).

Figure 6. Weighted correlation network analysis, differential expression analysis and univariate analysis of GDRGs. (A–D) Based on weighted correlation network analysis, 9 modules were accessed with a soft threshold = 8. (E) Correlation analysis between modules and clinicopathological data. (F) Differential expression analysis identified 258 differentially expressed GDRGs in 6 modules, with |log2(FC)| > 1 and FDR < 0.05. (G) Univariate analysis of differentially expressed GDRGs. GC: gastric cancer, GDRGs: GC differentiation-related gene, FDR: false discovery rate.

Figure 7. Generation, evaluation and validation of a prognostic risk scoring signature. (A) Kaplan-Meier analysis between the low-risk group and the high-risk group in the TCGA cohort. (B) Kaplan-Meier analysis between the low-risk group and the high-risk group in the GSE84437 dataset. (C) In the TCGA cohort, the areas under the ROC curves for predicting 1-year, 3-year and 5-year OS. (D) In the GSE84437 dataset, the areas under the ROC curves for predicting 1-year, 3-year and 5-year OS. (E) Expression levels of eight GDRGs in eight clusters. OS: overall survival, TCGA: The Cancer Genome Atlas, ROC: receiver operating characteristic, GDRGs: gastric cancer differentiation-related genes.

In addition, the expression levels of eight GDRGs in eight clusters are shown in Figure 7E. VCAN increased in cluster 7 (belonging to subset III); TNFAIP2, DUSP1 and ADAM8 in cluster 5 (belonging to subset III); STMN2 in cluster 3 (belonging to subset II); RNASE1 and TFF1 in cluster 2 (belonging to subset II); and AQP2 in cluster 1 (belonging to subset I).

Establishment and evaluation of a nomogram for predicting patient 3-year and 5-year OS

In the TCGA cohort, univariate analysis revealed that age, TNM stage, N stage and RS jointly affected patient prognosis (Figure 8A, all P < 0.05). Older patients tended to have later clinicopathological stage disease, and higher RS corresponded to poorer prognosis. Subsequently, multivariate analysis demonstrated that age and RS were independent prognostic factors for GC patients (Figure 8B, all P < 0.05). Four prognostic factors were combined to construct a nomogram for predicting 3-year and 5-year OS based on the TCGA cohort (Figure 8C). The areas under the ROC curves for predicting 3-year and 5-year OS were 0.727 and 0.730, respectively (Figure 8D). The calibration curves for predicting 3-year and 5-year OS were in good agreement with the observed values (Figure 8E).

Figure 8. Construction and evaluation of a nomogram based on the TCGA cohort. (A) Univariate analysis of risk score and clinicopathological characteristics. (B) Multivariate analysis of risk score and clinicopathological characteristics. (C) A nomogram for predicting 3-year and 5-year OS. (D) The areas under the ROC curves for predicting 3-year and 5-year OS. (E) The calibration curves for predicting 3-year and 5-year OS. TCGA: The Cancer Genome Atlas, OS: overall survival, ROC: receiver operating characteristic.

Acquisition and processing of scRNA-seq and bulk RNA-seq

The scRNA-seq data of 402 GC cells in 6 samples were obtained from the GSE112302 dataset in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. Afterwards, scRNA-seq data were initially processed by the ‘Seurat’ package, the percentage of mitochondrial genes was calculated through the PercentageFeatureSet function, and the relationship between sequencing depth and mitochondrial gene sequences and/or total intracellular sequences was elucidated by correlation analysis. Cells with a gene number < 100, sequencing number < 50 and mitochondrial gene content > 5% were excluded. After data filtering, scRNA-seq data were normalized by the LogNormalize method, and the top 1500 genes with highly variable features were identified by variance analysis. Moreover, the bulk RNA-seq data of 32 normal and 406 GC samples were accessed from TCGA (http://cancergenome.nih.gov/) database, and an additional 433 GC samples were obtained from the GSE84437 dataset (https://www.ncbi.nlm.nih.gov/geo/). The basic clinical information for each sample is detailed in Table 1. In this study, GC samples from patients with a survival time < 30 days, ambiguous survival status or unclear clinicopathological characteristics were excluded.

Dimensionality reduction and cell annotation

Under the condition of a false discovery rate (FDR) < 0.05, dimensions with significant separation were screened out through PCA [40], and then dimension reduction for the top 18 principal components (PCs) through the tSNE algorithm was employed to obtain the major clusters [41]. With the criteria of log₂ [fold change (FC)] > 0.5 and FDR < 0.05, marker genes in each cluster were accessed, and the top 10% of marker genes from clusters were laid out in the heatmap. Clusters were annotated through the ‘SingleR’ package based on marker genes.

Pseudotime and trajectory analysis

Pseudotime and trajectory analyses of GC cells were carried out by the ‘Monocle’ package [42]. Then, intracellular differentially expressed genes in cells with distinct differentiation states with |log2 (FC)| > 0.5 and FDR < 0.05 were designated GDRGs. The ‘clusterProfiler’, ‘org.Hs.eg.db’, ‘enrichplot’ and ‘ggplot2’ packages were used for GO annotation and KEGG enrichment analysis; only BPs were extracted in the GO annotation.

GDRG-based molecular subtypes of GC patients from the GSE84437 dataset

After log2-scale transformation of GDRG expression in the GSE84437 dataset, 1729 GDRGs were retained for GC molecular typing. The ‘ConsensusClusterPlus’ package, which provided quantitative and visual evidence of stability for measuring the number of unsupervised subtypes, was used for consensus clustering of GC [43]. The K-means algorithm and cumulative distribution function (CDF) curve were utilized to determine the best number of subtypes, and 50 iterations with maxK = 9 were conducted for stable subtypes. Kaplan-Meier analysis was completed to evaluate the survival results, and the ‘ggplot2’ package was used to display the proportion of clinicopathological features in each molecular subtype. Additionally, the expression levels of GDRGs in specific molecular subtypes were investigated within different cell differentiation trajectories and visualized in box plots.

Tumor microenvironment scores, immune cell infiltration and ICG expression across molecular subtypes

The study calculated the immune/stromal scores and tumor purity of each sample through the ‘ESTIMATE’ package. The content of 22 immune cells in each sample from GSE84437 was identified by CIBERSORT software. Meanwhile, the infiltration density of immune cells in different GC subtypes was compared by the ‘limma’ package, and a histogram of immune cells with significant differences was presented. Furthermore, 38 validated ICGs were summarized through extensive reading of the literature [7–15], and then their expression across GC subtypes was examined by differential expression analysis based on the GSE84437 dataset. Kaplan-Meier analysis was used to investigate the prognostic value of immune cells and ICGs.

Prognostic risk scoring signature generation and validation

In this study, the TCGA cohort was taken as the training set, and the GSE84437 dataset was taken as the validation set for prognostic RS signature generation and validation. First, the GDRGs in TCGA and GSE84437 cohorts were intersected, and transcription profiles were simultaneously normalized and corrected with log2-scale transformation. Subsequently, GDRGs were included in WGCNA, and key modules related to GC differentiation were identified based on correlation analysis. The genes in key modules were analyzed by differential expression analysis (with |log2(FC)| > 1 and FDR < 0.05) and univariate analysis (P < 0.05), and the remaining GDRGs were enrolled in multivariate Cox regression analysis to generate a GDRG-based prognostic RS signature. The mean value of the GDRG-based RS signature can be used to divide patients into high-risk and low-risk groups. The RS can be calculated as the sum of the products of GDRG expression levels and coefficients via the following formula:

$$RS={\displaystyle \sum _i^k\left(Ex{p}_i\times Co{e}_i\right)}$$

where ‘i’ and ‘k’ represent the ‘i’th gene and total number of genes, respectively. Kaplan-Meier analysis and ROC curves were used to evaluate the accuracy and prediction efficiency.

Nomogram construction based on TCGA cohort

Based on the TCGA cohort, clinicopathologic variables (e.g., age, sex, tumor grade and TNM staging) and RS were separately incorporated into univariate and multivariate analyses. Subsequently, prognostic variables were combined into a nomogram for predicting 3-year and 5-year OS. ROC curves and calibration curves were used to evaluate the predictive performance and accuracy of the nomogram.

Statistical analysis

Some continuous variables are expressed as the means ± standard deviation. The χ2 test and T-test/variance analysis were separately used to compare the distribution of dichotomous variables and continuous variables. Kaplan-Meier statistics and log-rank tests were used for survival analysis. All statistical analyses were carried out in R and Perl, and P <0.05 was considered to be statistically significant.

Data accessibility

The scRNA-seq data of GC samples were accessed from GEO database (GSE112302, https://www.ncbi.nlm.nih.gov/geo/). The bulk RNA-seq data of GC samples were obtained from GEO database (GSE84437, https://www.ncbi.nlm.nih.gov/geo/) and TCGA database (http://cancergenome.nih.gov/).

Ethics approval

Ethical approval was not required because the data came from publicly available databases.