Characterization of Immune-Related Molecular Subtypes and a Prognostic Signature Correlating With the Response to Immunotherapy in Patients With Gastric Cancer

Gastric cancer (GC) is a disease characterized by high molecular and phenotypic heterogeneity and represents a leading cause of cancer-related death worldwide. The tumor immune microenvironment (TIME) affects the response to immunotherapy and the prognosis of patients with GC. Explorations of the TIME in GC and characterization of molecular subtypes might enhance personalized treatment and facilitate clinical decision-making. In this study, two molecular subtypes were defined through unsupervised consensus clustering based on immune-related dysregulated genes. Then, patients with different molecular subtypes of GC were shown to have distinct differences in sensitivity to immune checkpoint blockers (ICBs). The immune-related prognostic signature was established utilizing least absolute shrinkage and selection operator (LASSO)-Cox regression analysis. Three independent external cohorts and the IMvigor210 cohort were introduced to validate the robustness of IPRS. scRNA-seq data of GC samples were used to decipher the underlying mechanisms of how IPRS contributes to the TIME. GC biospecimens were collected for RT-qPCR to further validate our findings. In summary, we characterized the abnormal TIME of GC and constructed a reliable immune-related prognostic signature correlating with the response to immunotherapy. This study may provide new strategies for developing individualized treatments for patients with GC.

Introduction

Gastric cancer (GC) is a disease with high molecular and phenotypic heterogeneity (1). Although GC is not one of the top malignancies in the United States (2), it represents a leading cause of cancer-related death worldwide. Surgery is currently considered the only curative option, but recurrence is common, even after complete resection (3). Moreover, most patients with GC develop advanced-stage disease because of the lack of specific signs of early gastric cancer, and some patients have missed the optimal surgical window when receiving the diagnosis (4). The benefit of chemotherapy varies from person to person because of primary or acquired drug resistance (5). Immune checkpoint blockers (ICBs), such as cytotoxic T lymphocyte-associated antigen 4 (CTLA-4) and programmed cell death 1 (PD-1) inhibitors, exert revolutionary effects on several tumors, while the efficacy seems to be closely related to the tumor immune microenvironment (TIME) (6, 7). Noncoding RNAs, such as microRNAs, long noncoding RNAs, and circular RNAs, as well as epigenetic alterations, including DNA methylation, histone acetylation, and chromatin remodeling, were reported to play vital roles in cancer development and resistance to therapeutic reagents (8, 9), but much more work is needed before their clinical application. Therefore, the identification of new molecular biomarkers to predict and improve the prognosis of patients with GC is essential.

The tumor immune microenvironment, which contains numerous cell types and the factors they secrete, plays critical roles in tumor growth, progression, and metastasis (10). According to the presence or absence of T cell-based inflammation, solid tumors are roughly characterized as “hot” (T cell-infiltrated) and “cold” (inflamed but not infiltrated or not inflamed) tumors (11, 12). Most patients with hot tumors exhibit greater sensitivity to ICBs, potentially because immune checkpoint inhibitors relieve the exhaustion of CD8+ T cells and renew their priming (13). Among various types of cells, cancer-associated fibroblasts, M2 macrophages, and regulatory T cells prevent CD8+ T cells from killing tumor cells by creating immunologic barriers. Hence, studies exploring and targeting the predominant components of the TIME may prolong the survival of patients.

Based on accumulating evidence, the classification of GC based on molecular subtypes rapidly identifies cancer characteristics and enhances the efficacy of personalized therapy (14). In the present study, we aimed to explore the TIME of GC and identify immune-related molecular subtypes. We estimated the immune and stromal components of the TIME of GC samples using gene expression profiles. Based on solid immune-related dysregulated genes, we characterized the molecular subtypes of GC and visualized the differences. We obtained insight into the correlation between immune-related molecular subtypes and the possible response to ICBs estimated using machine learning methods to determine if the identified subtypes would help improve risk stratification and guide therapeutic strategies. For the convenience of clinical applications of the molecular subtypes, we subsequently devoted ourselves to constructing a relatively streamlined and quantified prognostic signature. We calculated the immune-related prognostic risk scores (IPRSs) of patients using the established signature and investigated the relationship between the IPRS and clinical features, tumor-infiltrating immune cells, and the potential sensitivity to chemotherapy and immunotherapy. Single-cell RNA-seq data were introduced to decipher the underlying mechanisms of how IPRS contributes to the TIME. This study may provide guidelines to improve therapeutic strategies for patients with GC.

Materials and Methods

Data Collection and Processing

The RNA-sequencing profiles of 32 adjacent normal tissues and 375 primary gastric cancer (GC) tissues in the forms of counts and fragments per kilobase million (FPKM) were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) and subsequently transformed to log2(TPM+1). Meanwhile, the corresponding clinicopathological data from patients with GC in TCGA were obtained from the cBioPortal website (https://www.cbioportal.org/). The HMU-GC cohort [GEO accession: GSE184336 (15)], which takes advantage of RNA-sequencing technology, was employed in this study to verify the reproducibility of the results of the consensus clustering and immune subtypes. Batch effects caused by non-biotech bias between different datasets were removed using ComBat function in “sva” package. In addition, another three datasets, GSE14210 (16), GSE84437 (17), the Asian Cancer Research Group (ACRG) study [GEO accession: GSE66229 (18)], were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database and were utilized for validation with the log2 transformed gene expression matrices. The detailed criteria for inclusion of patients enrolled in this study were as follows: 1) histologically confirmed GC and 2) simultaneously available information on gene expression profiles and OS. The exclusion criteria were 1) patients with other diseases except for GC and 2) the follow-up time less than one month. The IMvigor210 dataset (19) containing 298 urothelial cancer cases who received anti-PD-L1 therapy was downloaded from http://research-pub.gene.com/IMvigor210CoreBiologies/. The processed scRNA-seq data of nine GC samples were acquired from professor Ying (Genome Sequence Archive accession: HRA000051, https://ngdc.cncb.ac.cn/gsa-human/) (20). A summary of the characteristics of GC patients in the four cohorts described above are shown in Table S1.

Identification of Immune-Related Dysregulated Genes

The immune scores and stromal scores were calculated using the “Estimation of STromal and Immune cells in MAlignant Tumours using Expression data” (ESTIMATE) algorithm to infer the fractions of immune and stromal cells in primary GC samples (21). According to the median score, GC samples were assigned into high and low immune/stromal score groups. Differentially expressed genes (DEGs) between GC samples and adjacent normal tissues or high and low immune/stromal score groups were screened using the R package “DESeq2” with the thresholds of log2fold change > 1 or < -1 and adjusted p values (padj) < 0.05. The overlapping genes among the three sets of DEGs described above were considered immune-related dysregulated genes (IDGs).

Functional Enrichment Analysis

Gene set enrichment analysis (GSEA), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted using the R package “clusterProfiler” (22). Annotated gene sets in “h.all.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt” obtained from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/index.jsp) were selected as the reference gene sets in GSEA. Enriched terms with a false discovery rate (FDR) < 0.05 were considered significant and were visualized using the R packages “enrichplot” and “ggplot2”.

Characterization of Molecular Subtypes and Estimation of Tumor-Infiltrating Immune Cells

Based on IDGs, the R package “ConsensusClusterPlus” was employed to characterize the immune-related molecular subtypes of primary GC samples. Principal component analysis (PCA) was subsequently conducted to verify the differences among subtypes. Three algorithms, including single-sample GSEA (ssGSEA), Tumor Immune Estimation Resource (TIMER), and CIBERSORT, were utilized to estimate the fractions of tumor-infiltrating immune cells. The signatures of 28 immune cells were collected from TISIDB (23), which were used to calculate the abundance of tumor-infiltrating immune cells through ssGSEA.

Construction and Validation of Immune-Related Prognostic Signatures

For the convenience of clinical application, we decided to develop a relatively streamlined prognostic risk model composed of gene expression levels and their respective coefficients. Briefly, a univariate Cox regression analysis was first conducted to screen the immune-related prognostic dysregulated genes with the threshold of a p value < 0.05 based on primary GC samples acquired from TCGA database, their corresponding clinical information, and IDGs previously identified. For the accuracy of model construction, patients who were followed for less than one month were excluded. A total of 338 remaining samples were randomly assigned into a training set (237 patients) and an internal validation set (101 patients) at a 7:3 ratio. The clinical features of the two sets are summarized in Table S2. Then, the least absolute shrinkage and selection operator (LASSO) regression algorithm was performed to screen the immune-related prognostic signatures with a minimum 10-fold cross-validation based on the training set and prognostic IDGs using the R package “glmnet”. Finally, prognostic risk models were constructed based on genes extracted from the LASSO regression analysis and their respective coefficients (β). Kaplan-Meier Plotter (http://kmplot.com/analysis/), a web-based survival analysis tool, was utilized to confirm the prognostic power of selected genes in this study (24). The immune-related prognostic risk score (IPRS) of each patient was calculated using the following formula: IPRS = expression level of gene₁ × β_gene1 + expression level of gene₂ × β_gene2 + …expression level of gene_n × β_genen. The internal validation set, the whole TCGA set, and four independent external sets consisting of GSE14210, GSE88437, ACRG, and IMvigor210 were used for validation. Patients in each set were divided into high-IPRS and low-IPRS groups according to the optimal IPRS cutoff. Kaplan–Meier curves, receiver operating characteristic (ROC) curves, and time-dependent area under the curve (AUC) were introduced to evaluate the robustness of immune-related prognostic signatures.

Prediction of the Sensitivity to Chemotherapy and Immune Checkpoint Blockers

Benefiting from the application of machine learning in medicine, people have predicted the potential therapeutic effects of different treatments. oncoPredict, an R package for predicting the drug response of patients with cancer that has been applied to various in vitro and in vivo contexts for drug and biomarker discovery, was employed to evaluate the sensitivity of patients with GC to common chemotherapeutic drugs (25). The candidate drugs, including 5-fluorouracil, Docetaxel, Paclitaxel, Epirubicin, Irinotecan, Cisplatin, and Oxaliplatin, were selected according to the Food and Drug Administration (https://www.cancer.gov/about-cancer/treatment/drugs/stomach) and the Chinese Society of Clinical Oncology (CSCO). The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm, which was developed by Jiang P et al. through modeling two primary mechanisms of tumor immune evasion for predicting ICB response (26), was introduced to estimate the potential response of patients with GC to ICB therapy.

Single-Cell RNA Sequencing Data Analyses

The processed scRNA-seq data of nine GC samples were acquired from professor Ying (20). The raw gene expression matrices were processed using the R package “Seurat” (version 4.1.0) with a standard workflow including quality control, dimension reduction, clustering, and cell type annotation. Quality control criteria were as follows: 1) cells had either fewer than 500 RNA features, 2) over 60,000 or less than 500 RNA counts, 3) over 50% RNA features derived from the mitochondrial genome, or 4) less than 3% RNA features derived from the ribosome genes were removed. Gene expression matrices of the remaining 23,511 cells were normalized and subsequently scaled by regressing out the cell-cycle signature scores. Data integration and batch effect removal steps were conducted using R package “harmony” after PCA. Top 20 dimensions from harmony were selected to perform UMAP dimensionality reduction, followed by cell clustering with a resulotion value of 0.3. Clusters with cell numbers less than 50 were deserted. A total of 14 main clusters were generated and annotated to 9 cell types according to marker genes from the origin article (20). Cell-cell communication were evaluated using CellChat and iTALK tools. Top ligand-receptor pairs and signaling pathways were displayed.

Real-Time Quantitative Polymerase Chain Reaction

Seven pairs of GC samples and adjacent normal tissues were obtained from Shanghai East Hospital Biobank. This research was approved by the Ethics Committee of East Hospital Affiliated Tongji University, Tongji University School of Medicine (2020-053). All patients signed informed consent forms before donating their specimens. Total RNA was extracted using TRIpure Total RNA Extraction Reagent (ELK Biotechnology, EP013) according to the manufacturer’s instructions. Reverse transcription was performed using M-MLV Reverse Transcriptase (ELK Biotechnology, EQ002). RT–qPCR was performed with QuFast SYBR Green PCR Master Mix (ELK Biotechnology, EQ001) and a StepOne™ Real-Time PCR System (Life Technologies). GAPDH was selected as the internal reference. The relative expression level was calculated using the 2^-ΔΔCT method. The primer sequences are presented in Table S3.

Statistical Analysis

Statistical tests were conducted using R software (version 4.1.0). Categorical variables were compared using chi-square tests. Comparisons between two groups were performed using the Wilcoxon rank-sum test. Univariate and multivariate Cox regression analyses were conducted to screen the independent prognostic factors. Correlation analyses were performed with the Pearson correlation test. Kaplan–Meier curves for overall survival (OS) were plotted, and the difference between groups was compared using the log-rank test. The random-effects meta-analysis model was employed to calculate the pooled hazard ratio (HR) with the R package “meta”. A P value < 0.05 was considered statistically significant.

Results

Identification of Tumor Immune Microenvironment-Related Dysregulated Genes

A flow diagram of overall analyses was displayed in Figure 1. The tumor immune microenvironment (TIME) favors the growth and progression of cancer cells, affecting the clinical outcomes of patients (27, 28). In the present study, we assessed the TIME of each GC sample by calculating the immune and stromal scores. Subsequently, DEGs between the low and high immune/stromal score groups were screened. Compared with the low immune score group, 2092 upregulated DEGs and 1630 downregulated DEGs were identified in the high immune score group. Meanwhile, 3104 upregulated and 875 downregulated DEGs were identified in the high stromal score group compared to the low stromal score group. In addition, DEGs between GC samples and normal tissues were also calculated, and 10724 dysregulated DEGs were screened (Figure S1 and Table S4). We took the intersection among the three sets of DEGs described above to investigate the immune-related dysregulated genes (IDGs). As shown in Figure 2A, 1083 overlapping genes were identified and considered IDGs (Table S4). We performed functional enrichment analyses to further explore the function of IDGs. As a result, the main enriched GO terms were associated with immune responses, including humoral immune response, immunoglobulin complex, and antigen binding. KEGG pathways were mainly enriched in neuroactive ligand–receptor interactions and calcium signaling pathways (Figure 2B, C, and Table S5). Based on these results, the 1083 overlapping genes were robust IDGs.

FIGURE 1

Figure 1 The flow diagram of overall analysis.

FIGURE 2

Figure 2 Identification of IDGs and immune-related clusters. (A) Screening the IDGs by Venn diagram. (B, C) GO (B) and KEGG (C) enrichment analyses of IDGs. (D) The cumulative distribution function (CDF) curves for k = 2 to 6. (E) Three clusters were identified through consensus clustering. (F) PCA analysis revealed the dissimilarity among three clusters. BP, biological process; CC, cellular component; MF, molecular function.

Characterization of Immune-Related Molecular Subtypes

Previously, the tumor microenvironment subtypes have been shown to be correlated with the responses of patients with various cancers to immunotherapy, and patients with immune-favorable subtypes tend to benefit more from immunotherapy (29). Emerging evidence indicates that differences in the molecular pathology of indistinguishable cancers affect the clinical features of the disease, and molecular subtypes now guide clinical therapeutic strategies for multiple cancers (30, 31). Based on IDGs, GC samples were assigned to three clusters using consensus clustering with an optimal k of 3 (Figures 2D, E). PCA revealed distinct differences among these three clusters (Figure 2F).

To further explore the relationship between clusters and immune activities, we firstly applied TIMER algorithm to estimate the levels of tumor-infiltrating immune cells. Compared with GC samples in Cluster_1, samples in Cluster_2 and Cluster_3 exhibited significantly higher fractions of all the estimated six types of immune cells, while the discrepancies in levels of immune cells between Cluster_2 and Cluster_3 were not so distinct (Figure 3A). In addition, the abundance of CD8+ T cells and myeloid dendritic cells between Cluster_2 and Cluster_3 showed no statistical differences. This tendency was more obvious in the HMU-GC cohort (Figures S2A–E). On the other hand, we took advantage of TIDE algorithm and inferred the possible sensitivities to ICB to compare the potential differences in immunotherapeutic efficiency among the identified three clusters. Patients in Cluster_1 showed distinctly lower TIDE scores than any other clusters, meaning that patients in Cluster_1 may benefit more from ICB (Figure 3B). Agreeing with the results of immune cell analyses, the TIDE scores between Cluster_2 and Cluster_3 also showed no significant differences (Figure 3B), suggesting that patients in the Cluster_2 and Cluster_3 were in the similar immune states. Therefore, we manually classified the GC patients in Cluster_1 as Immune_L subtype other ones belonged to the Immune_H subtype. Significant higher fractions of all the six kinds of immune cells were observed in the Immune_H subtype compared to the Immune_L subtype (Figure 3C). The TIDE scores between the two immune subtypes were also obviously discrepant (Figure 3D). The characterized immune subtypes were further confirmed by the abundance of 28 immune cells calculated through the ssGSEA algorithm (Figure 3E). To ensure the reproducibility of the results of consensus clustering and immune subtypes, we introduced the HMU-GC cohort and found the almost same results (Figures S2A–H). Next, we compared the identified immune subtype with previously reported subtypes in GC, such as the TCGA molecular subtype (32) and the Lauren class. As displayed in Figure 3F, the alluvial diagram indicated that the immune subtype was associated with the TCGA molecular subtype, the Lauren class, and tumor grade, while having less connection to the TNM stage. Taken together, GC patients exhibited a certain degree of heterogeneity in immune status, and the immune subtype we identified could partly reflect the levels of tumor-infiltrating immune cells and sensitivities to ICB, which may help guide clinical decision-making.

FIGURE 3

Figure 3 Exploration of the TIME and characterization of immune-related molecular subtypes. (A, B) The fractions of tumor-infiltrating cells estimated by TIMER (A) and TIDE scores (B) among three clusters. (C, D) The fractions of tumor-infiltrating cells estimated by TIMER (C) and TIDE scores (D) between the two immune subtypes. (E) The levels of 28 immune cells of GC samples in the two immune subtypes were calculated through ssGSEA. (F) The relationships between immune subtypes and other previously reported subtypes in GC. NS. or ns, no statistical significance, *p < 0.05, **p < 0.01, ***p < 0.001.

Construction of a Prognostic Signature Based on Prognostic IDGs

For a better clinical application of identified immune-related molecular subtypes, we constructed a relatively simple and quantified prognostic gene signature. We first conducted a univariate Cox regression analysis on 1083 previously identified IDGs, and 63 prognostic IDGs were screened (Table S6). Then, we fit the aforementioned prognostic IDGs into a LASSO regression model with the training set to identify the optimal prognostic signature (Figure 4A). Subsequently, 10-fold cross-validation was performed to overcome the overfitting effect. A log lambda of -2.742 was selected with the lowest partial likelihood deviance (Figure 4B). A panel of five genes remained with nonzero coefficients (Figures 4C, D and Table S7). The survival analysis revealed that patients with higher expression levels of each of these five genes experienced worse clinical outcomes (Figure 4E). The results from Kaplan-Meier Plotter also confirmed the prognostic power of these five genes in GC (Figure S3). In addition, the multivariate Cox regression analysis indicated that four of five genes, namely, GCG, ITIH2, CYTL1, and MAGEA11, were independent risk factors for the OS of patients with GC (Figure 4C). The formula for calculating the immune-related prognostic risk score (IPRS) from the prognostic signature is as follows: IPRS = (0.12733) × Exp_CDH10 + (0.07203) × Exp_GCG + (0.08738) × Exp_ITIH2 + (0.17213) × Exp_CYTL1 + (0.05973) × Exp_MAGEA11.

FIGURE 4

Figure 4 Construction of an immune-related prognostic gene signature. (A) The changing trajectory of each variable in LASSO-Cox regression analysis. (B) Selection of the optimal lambda value with the minimum partial likelihood deviance. (C) Forest plot of multivariate Cox analysis of five selected genes based on the training set. (D) The coefficients of five selected genes. (E) Kaplan-Meier curves with log-rank p-values of the five selected genes. *p < 0.05, **p < 0.01.

Validation and Evaluation of the Prognostic Signature

Patients in each set were divided into high and low IPRS groups according to the optimal cutoff value to validate the robustness of the prognostic signature. Kaplan–Meier curves revealed that the IPRS performed well in distinguishing patients with good or unfavorable overall survival not only in the training set (Figure 5A) but also in the internal validation set and the whole TCGA set (Figures 5B, C). We determined the accuracy of the IPRS and traditional clinical features by plotting ROC curves and calculating the AUC value for OS based on the whole TCGA set. The AUC value for IPRS was 0.745 at 3 years and 0.771 at 5 years (Figure 5D), which suggested promising predictive power. Compared to clinical characteristics, the IPRS exhibited the highest accuracy for predicting OS at 5 years (Figure 5E). In addition, we calculated the time-dependent AUC value of each factor to evaluate the overall predictive ability. As shown in Figure 5F, the AUC value of IPRS was approximately 0.75 and ranked at the top most of the time, indicating that the IPRS might serve as a supplement to the AJCC staging system for improving the prognosis of patients with GC. On the other hand, we conducted univariate and multivariate Cox regression analyses to assess the prognostic value for the OS of patients with GC. The hazard ratio (HR) and 95% confidence interval (CI) of the IPRS were 6.07 and 3.28-11.22, respectively, in the univariate Cox regression analysis (Figure S4A). The IPRS was identified as the only independent risk factor in the multivariate Cox regression analysis (HR = 5.11, 95% CI = 2.35-11.09, Figure S4B). Taken together, the IPRS exhibited great prognostic value and might help improve the clinical outcomes of patients with GC.

FIGURE 5

Figure 5 Validation and assessment of the immune-related prognostic gene signature. (A–C) Kaplan-Meier curves of low and high IPRS groups in the training set (A), internal validation set (B), and the whole TCGA set (C). (D) ROC curves of the IPRS at 3 years and 5 years. (E) ROC curves of the IPRS and traditional clinical characteristics at 5-year-survival. (F) Time-dependent AUC curves of the IPRS and clinical factors.

Correlation Between the IPRS and Clinicopathological Features

We generated a series of boxplots based on the whole TCGA set to investigate and illustrate the relationships between the IPRS and clinicopathological characteristics (Figure 6). Among the various clinical features, a higher IPRS was associated with male sex (p < 0.05), distant metastasis (p < 0.01), and poor clinical outcomes (living status and disease-free status, both p < 0.001) but was unrelated to age and lymph node metastasis. Moreover, the IPRS exhibited a certain discriminatory power for the tumor grade (G3 vs. G2, p < 0.01), tumor stage (stage IV vs. lower stages, all p < 0.05), and pathological T stage (T1 vs. higher T stage levels, all p < 0.05). Besides, IPRS in the Immune_H subtype was significantly higher than that in the Immune_L subtype (p < 0.01, Figure S5A). Based on these results, the IPRS that was established from gene signatures partially reflects the clinical features of patients with GC and has the ability to distinguish patients with relatively early-stage GC from advanced-stage tumors, which might assist with clinical decision-making.

FIGURE 6

Figure 6 The relationships between IPRS and clinicopathological features.

IPRS May Predict the Responses to Chemotherapy and ICB

Because the IPRS was developed based on prognostic IDGs, we explored the association of the IPRS with tumor-infiltrating immune cells. The IPRS exhibited a positive correlation with the fractions of naïve B cells, M2 macrophages, and resting mast cells (Figure 7A). In contrast, the IPRS was negatively correlated with the fractions of CD8+ T cells and activated CD4+ memory T cells (Figure 7A). We determined whether the IPRS predicted the potential sensitivity of patients to chemotherapy and ICBs based on the immune-related molecular subtypes previously identified by comparing the IC50 values of common chemical drugs and TIDE scores between the low and high IPRS groups. As we expected, patients with GC presenting a lower IPRS tended to be more sensitive to various chemotherapeutic agents for GC, such as 5-fluorouracil, Paclitaxel, Docetaxel, and Oxaliplatin (Figure 7B). Setting the default value of 0 as the threshold (26), patients with GC were considered responsive or nonresponsive to ICBs according to the individual TIDE scores. The distribution of the TIDE score of each patient in the low and high IPRS groups is shown in Figure 7C. In addition, the chi-square test indicated that patients with GC in the low IPRS group were more likely to respond to ICBs (chi-square test, p < 0.001). Thus, the IPRS might predict the responses to chemotherapy and ICBs.

FIGURE 7

Figure 7 Investigation of the correlation between IPRS and tumor-infiltrating cells and the sensitivity to chemotherapies and ICBs. (A) IPRS was positively or negatively correlated with some tumor-infiltrating cell types. (B) The estimated IC50 for various chemotherapeutic drugs. (C) The distributions and proportions of potential responders and non-responders were estimated by TIDE scores between low and high IPRS groups. NS., no statistical significance, **p < 0.01, ***p < 0.001.

GSEA was performed by setting the gene signatures in “h.all.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt” as references to investigate the potential regulatory mechanisms responsible for the differences between the low and high IPRS groups. Hallmarks, including angiogenesis, coagulation, the epithelial-mesenchymal transition, and myogenesis, were significantly enriched in the high IPRS group (Figure 8A and Table S8). KEGG pathways such as neuroactive ligand-receptor interaction, complement, coagulation cascades, and adipocytokine signaling pathways were positively correlated with the IPRS, while pathways such as mismatch repair and protein export were negatively regulated in the high IPRS group (Figure 8B and Table S8). Since the IPRS was capable of predicting the possible sensitivities to ICB, we wondered if the IPRS had connections to the classical predictive markers of immunotherapy, such as immune checkpoints, TMB (Tumor Mutation Burden), and MSI (Microsatellite Instable) status. As a result, most of the immune checkpoints were low expressed (Figure S5B), and the scores of immune checkpoints exhibited no obvious differences in different risk groups (Figure S5C). While, the TMB was negatively correlated with the IPRS (R = -0.3, p < 0.001, Figure 8C). In addition, the IPRS showed significant discrepancies among different MSI statuses (Figure 8D). The relationships between IPRS and predictive markers may be responsible for the capability of IPRS to predict the possible response to ICB in GC.

FIGURE 8

Figure 8 Exploration of the regulatory mechanisms resulting in the differences between low and high IPRS groups. (A, B) Hallmarks of cancers (A) and KEGG pathways (B) were investigated between low and high IPRS groups through GSEA. (C) IPRS was negatively correlated with the TMB. (D) IPRS showed significant discrepancies among different MSI statuses. MSI, Microsatellite Instable; MSS, Microsatellite Stable.

External Validation of the Immune-Related Prognostic Signature

Three independent external GC cohorts (GSE14210, GSE84437, and ACRG) were utilized to further validate the established prognostic signature. Patients in each cohort were divided into high and low IPRS groups according to the optimal cutoff value. Kaplan–Meier curves showed that the immune-related prognostic signature performed well in discriminating patients with favorable and poor OS (Figures 9A–C). For the DFS (disease-free survival), the IPRS also exhibited a certain distinguishing power (Figures 9D–F). Along with TCGA set, the meta-analysis revealed that the pooled HR and 95% CI of the IPRS were 2.02 and 1.61-2.53, respectively (Figure 9G). Besides, we compared the efficacy of the immune-related prognostic signature with other three prognostic signatures, including a six-gene pyroptosis-related signature (33), a two-gene signature (34), and a five-gene glycolysis-related signature (35). Among the four prognostic signatures, the IPRS showed the highest accuracy in predicting the 3-year-survival and 5-year-survival (Figures S6A, B). By generating a time-dependent AUC plot, we found that the efficacy of IPRS ranked at the top in the first five years and was surpassed by glycolysis-related signature at the seventh year (Figure S6C). To examine the utility of IPRS in speculating the benefits from immunotherapy, we employed the IMvigor210 cohort consisting of patients who received anti-PD-L1 therapy. It is as expected that patients with low IPRS were more likely to achieve CR/RR and survived longer than patients with high IPRS (Figures 9H, I). The objective response rate of anti-PD-L1 therapy was significantly higher in the low IPRS group than that in the high IPRS group (chi-square test, p = 0.002, Figure 9J). Therefore, our results revealed that the immune-related prognostic signature was robust and reliable.

FIGURE 9

Figure 9 External validation of the immune-related prognostic gene signature. (A–C) Kaplan-Meier analyses of the low and high IPRS groups in three independent external cohorts on OS, including GSE14210 (A), GSE84437 (B), ACRG (C). (D–F) Kaplan-Meier analyses of the low and high IPRS groups in three cohorts on DFS, including TCGA-STAD (D), GSE14210 (E), ACRG (F). (G) Meta-analysis was conducted to evaluate the pooled HR of the immune-related prognostic signature. (H) IPRS in groups with a different anti-PD-L1 response status. (I) Kaplan-Meier analysis of patients with high and low IPRS in the IMvigor210 cohort. (J) The objective rate of clinical response to anti-PD-L1 immunotherapy in high and low IPRS groups in the IMvigor210 cohort. CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease. *p < 0.05.

Experimental Verification of the Aberrant Expression of Screened IDGs in Biospecimens

As mentioned above, a panel of five genes was screened to construct the prognostic signature. Among them, CDH10, GCG, and CYTL1 were downregulated in GC samples compared to adjacent normal tissues, while TITH2 and MAGEA11 were upregulated (Figure 10A and Table S4). We conducted experimental validation using biospecimens to confirm the expression patterns of these five genes. All genes showed the same patterns of expression as TCGA gene expression profiles except GCG (Figure 10B), which interested us and suggested further research. By analyzing the gastric cancer samples and paired adjacent normal tissues in TCGA, we found that these 27 paired samples exhibited different GCG expression patterns, which in turn reflected the high heterogeneity of GC (Figure S7).

FIGURE 10

Figure 10 Validation of the aberrant expression of the five selected IDGs in biospecimens. (A) The expression levels of five selected IDGs in adjacent tissues and gastric cancers obtained from TCGA. (B) RT-qPCR was performed using biospecimens to validate the expression pattern of five selected IDGs. *p < 0.05, **p < 0.01, ***p < 0.001.

Decipher the Underlying Mechanisms of How IPRS Contributes to TIME at the Single-Cell Level

The scRNA-seq data of nine GC samples were analyzed to decipher the underlying mechanisms of how IPRS might contribute to the TIME (Figure 11A). A total of 23,511 cells remained after quality control and were used for subsequent analyses. According to canonical marker genes, 14 main clusters with cell numbers over 50 were annotated to nine cell types (Figures 11B–D): B cell (4,094 cells, 17.4%, marked with CD79A and MS4A1), T cell (6,071 cells, 25.9%, marked with CD2, CD3D, and CD3E), endocrine cell (234 cells, 1.0%, markered with CHGA), endothelial cell (231 cells, 0.9%, marked with ENG and VWF), epithelium (9,992 cells, 42.6%, marked with EPCAM, KRT18, and KRT8), fibroblast (327 cells, 1.4%, marked with ACTA2 and COL1A2), macrophage (1,015 cells, 4.3%, marked with CD14 and CD68), mast cell (500 cells, 2.1%, marked with CPA3 and KIT), plasma cell (1,012 cells, 4.3%, marked with CD79A and MZB1). T cells were subclustered into CD4+ T cells or CD8+ T cells according to CD4 or CD8A expression levels. To reveal the potential mechanisms of IPRS, we firstly checked the expression levels of five selected genes, including CYTL1, GCG, CDH10, ITIH2, and MAGEA11. CYTL1 was predominantly expressed by endothelial cells, while GCG was mainly expressed by endocrine cells (Figure 11E). Other genes were either with extremely low expression patterns or not detectable. Hence, the functions of endothelial cells and endocrine cells in the TIME were the key links to the underlying mechanisms of IPRS. Cell-cell communication indicated that endocrine cells had impacts on cancer cells through strong MDK-SDC4 and MDK-NCL interactions, while endothelial cells affected tumor cells mainly through COL4A2-SDC4 and COL4A1-SDC4 interactions (Figures 11F, H). SDC4 was reported to promote cancer cell progression and angiogenesis (36). NCL was found to play a role in promoting neuroblastoma tumorigenesis (37). These results revealed that endothelial cells and endocrine cells in the TIME could foster cancer development. Besides, IGFBP4-LRP6 interactions between endothelial cells and tumor cells (Figure 11G) might also promote tumor formation and progression (38). More importantly, high expression levels of HAVCR2 on immune cells including CD8+ T cells, CD4+ T cells, and macrophages plus high expression levels of LGALS9 on endothelial cells and endocrine cells (Figure 11I) indicated that endothelial cells and endocrine cells may suppress immune cells activation and induce immunosuppressive tumor microenvironment. Therefore, the underlying mechanisms of IPRS could be summarized as promoting tumor growth and suppressing the anti-tumor functions of tumor-infiltrating immune cells.

FIGURE 11

Figure 11 Decipher the underlying mechanisms of how IPRS contributes to the TIME at the single-cell level. (A) UMAP plot showing sample origin after batch effect removal. (B) UMAP plot showing the expression levels of canonical marker genes for nine cell types. (C) Violin plot showing the expression distribution of canonical marker genes for nine cell types. (D) UMAP plot showing the distribution of nine cell types. (E) Bubble chart showing the expression patterns of genes involved in the prognostic model. (F) Bubble chart reflecting top ligand-receptor interactions of endothelial cells and endocrine cells communicating with other cell types. (G–I) Circle plots showing top ligand-receptor pairs in growth factor module (G), cytokine module (H), and immune checkpoint module (I) of cell-cell communication networks. Line width positively correlates with the expression levels of ligands; arrow width positively correlates with the expression levels of receptors.

Discussion

Gastric cancer is considered a highly molecularly heterogeneous disease with an extremely poor prognosis (1). Recently, molecular subtype classification has emerged as a new strategy for the treatment of tumors (39–41). Benefiting from the development of next-generation sequencing and other genomic technologies, substantial progress has been achieved in the molecular classification of GC in the past few years (32, 42). Four molecular subtypes, including Epstein–Barr virus-associated, microsatellite instable, chromosomal instable, and genomically stable carcinomas, were identified by the Cancer Genome Atlas Research Network (40), which facilitated the exploration of novel targeted therapeutics. In addition, a number of potential biomarkers or signatures have been identified or constructed based on sequencing profiles for predicting clinical outcomes and therapeutic effects to improve the prognosis of patients with GC (43–45). However, novel molecular biomarkers are still urgently needed to aid in the clinical management of patients with GC.

The tumor immune microenvironment fosters tumor progression and metastasis and is considered a novel potential therapeutic target in various cancers, including GC (46). Based on accumulating evidence, the TIME of GC is very specific and suitable for promoting the growth and expansion of cancer cells (47). We first evaluated the levels of immune and stromal components and quantified them using the immune/stromal score with the ESTIMATE algorithm to obtain insights into the TIME of GC (21). A total of 1083 dysregulated genes were screened and validated as robust immune-related genes via functional annotation. Taking advantage of consensus clustering, we characterized three clusters and two immune-related molecular subtypes of GC based on IDGs. Estimates of tumor-infiltrating cells using multiple algorithms all indicated that gastric tumors in Immune_H subtype showed higher immune response activity and a higher percentage of CD8+ T cells, which were considered relatively “hot” tumors (11). Surprisingly, the relatively “hot” tumors were evaluated to have higher TIDE scores, indicating that their sensitivity to ICBs was lower (Figures 3C, D and Figures S2F, G). Previous studies discovered two main mechanisms of tumor immune evasion (11, 48). In some tumors, T cells are blocked by immunosuppressive factors, leading to a low level of infiltrating T cells and insensitivity to ICBs. Other tumors may have adequate cytotoxic T cells, but they might also exhibit a lower response to ICBs because of the dysfunction of these T cells. Although the fractions of CD8+ T cells in relatively “hot” tumors were high, most of them were in a dysfunctional state, resulting in a lower TIDE score (Figures 3C, D and Figures S2F, G). Hence, a classification solely based on the level of infiltrating cytotoxic T cells is too rough to guide immunotherapy, at least in patients with gastric tumors. While the two immune-related molecular subtypes may serve as a potential biomarker for clinical decision-making.

However, the stratification of patients with GC into the two identified molecular subtypes is costly since hundreds of genes must be detected. Researchers explored a lot of popular approaches to construct practical prognostic models for patients with other cancer types, including immune-related lncRNA signatures (49), EMT-related lncRNA signatures (50), m5C-related lncRNA signatures (51), and mutation-derived genome instability-related lncRNA signatures (52). For the convenience of application, we subsequently developed a quantified model to evaluate the risk of individual patients with GC using as few genes as possible. Ultimately, a panel of five IDGs was screened through variable selection with the Cox regression model (53, 54). Genes included in the immune-related prognostic signature all showed positive or negative regulation of tumor progression, such as CDH10, which encodes a type II classical cadherin of the cadherin superfamily that was shown to be frequently mutated and associated with various cancers, such as pancreatic ductal adenocarcinomas (55), gastric and colorectal cancers (56), and lung squamous cell carcinomas (57). The protein encoded by the GCG gene is glucagon, which regulates blood glucose levels by increasing gluconeogenesis and reducing glycolysis (58). Glucagon-induced hyperglycemia promotes tumor growth and angiogenesis in mice (59). As shown in the current study, patients with GC presenting high GCG expression experienced worse clinical outcomes, agreeing with a previous study conducted using mice (59). It’s surprising to find that the expression pattern of GCG in GC samples from TCGA was inconsistent with our biospecimens. To explain this phenomenon, we analyzed the GC samples and paired adjacent normal tissues in TCGA, finding that a number of pairs of samples did show the different GCG expression patterns, which may be due to the high heterogeneity of GC. ITIH2 plays important roles in extracellular matrix stabilization and in the prevention of tumor metastasis (60). In addition, overexpression of the deregulated ITIH2 protein in glioma cells not only inhibits cancer cell invasion but also suppresses cell proliferation and promotes cell-cell adhesion (61). However, ITIH2 seemed to be a risk factor for patients with GC in our study, suggesting the need for more experiments. Previous studies have revealed that CYTL1 plays opposite roles in distinct cancer types (62). CYTL1 is considered a tumor suppressor in breast cancer by inhibiting metabolic reprogramming (63). Meanwhile, CYTL1 was also identified to be associated with the growth and metastasis of neuroblastoma cells, together with its role in vessel formation (64). MAGEA11, a cancer germline antigen, is correlated with tumor progression, drug resistance, and poor prognosis in human cancers (65, 66).

Integrating the five genes above and their individual coefficients, we calculated the IPRS of patients with GC and assigned them into high and low IPRS groups. Patients with a high IPRS experienced a significantly shorter OS. The ROC curve analysis revealed that IPRS was superior to traditional clinical features, with an AUC value of 0.771 at five years. Time-dependent AUC values showed that the IPRS and TNM stage were the top two accurate predictors and that IPRS ranked first most of the time. Univariate and multivariate Cox analyses indicated that IPRS was the only independent risk factor for the OS of patients with GC. In addition, we employed three independent external cohorts to validate the reliability of the prognostic signature and found that the IPRS still worked well. Taken together, IPRS showed promising prognostic performance and might represent a supplement for the TNM stage in the management of patients with GC.

Since the IPRS originated from immune-related molecular subtypes, we explored the relationships between the IPRS and TIME. Correlation analyses revealed that the IPRS was positively associated with the fractions of infiltrating naïve B cells, M2 macrophages, and resting mast cells. In contrast, the IPRS seemed to negatively correlate with the fractions of CD8+ T cells and activated memory CD4+ T cells. Mature B cells are considered to suppress tumor progression through various mechanisms, such as secreting immunoglobulins, activating the T cell response, and killing tumor cells directly (67), while naïve B cells have not been exposed to antigens and perform limited functions in the TIME. Macrophages are a plastic cell type that adopts either pro- (M1-like) or anti-inflammatory (M2-like) phenotypes in response to signals from the TIME (68). Existing studies show that M2-like tumor-associated macrophages play a central role in tumor development through their contributions to basement membrane breakdown, angiogenesis, and immune suppression (69). Mast cells interact with infiltrated immune cells and tumor cells through cell-to-cell interactions and promote neovascularization and tumor invasion (70). Infiltrated CD8+ T cells are the primary target cells of immunotherapy and are capable of directly killing cancer cells (13). Additionally, memory CD4+ T cells also participate in antitumor immunity (71). Considering the correlation between the IPRS and infiltrated immune cells, the finding that the immune-related signature performed well in risk stratification and predicting immunotherapy sensitivity is not surprising.

Angiogenesis and the epithelial-mesenchymal transition, two representative hallmarks of cancers, were enriched in the high IPRS groups. Angiogenesis is induced by tumors to satisfy the needs for continuous nutrients and evacuating metabolic wastes (72). Bevacizumab, a humanized anti-VEGF-A monoclonal antibody, significantly improves the prognosis of patients with GC when administered in combination with chemotherapy (73). The epithelial-mesenchymal transition, a process in which epithelial cells gain mesenchymal features, is correlated with tumor invasion, metastasis, and resistance to drugs and apoptotic stimuli (74, 75). These findings may explain why patients with a high IPRS experienced poor clinical outcomes and were insensitive to chemotherapeutic reagents. Immune checkpoints, TMB, and MSI status are typical predictive markers of immunotherapy. IPRS showed no correlation with checkpoints but exhibited significant connections to TMB and MSI status, which may be responsible for the ability of IPRS to predict the potential response to ICB.

scRNA-seq data analyses revealed that CYTL1 and GCG are mainly expressed by endothelial cells and endocrine cells, respectively. Tumor endothelial cells could contribute to cancer progression by facilitating angiogenesis (76), which agrees with our bioinformatic analyses. In addition, the results of cell-cell communication indicated that endothelial cells and endocrine cells in the tumor microenvironment may foster tumor cell proliferation by targeting the SDC4 receptor while suppressing immune cell activation through the HAVCR2 receptor. These results further explained the underlying mechanisms of IPRS.

However, our research also had some limitations. First, our analysis was mainly based on public retrospective datasets, and the results require further validation by prospective studies in the future. Second, we failed to validate whether ICB-treated patients with different molecular subtypes of GC or IPRS benefit differently due to the lack of expression data for patients with GC undergoing ICB treatment.

In summary, we identified two immune-related molecular subtypes of GC with distinct tumor microenvironments. We also estimated the potential responses to immunotherapy between patients with the two molecular subtypes. In addition, we constructed a robust and relatively streamlined prognostic gene signature for the convenience of clinical application. Our study may provide new insights into personalized treatment and contribute to improving the prognosis of patients with GC.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics Statement

This research was approved by the Ethics Committee of East Hospital Affiliated Tongji University, Tongji University School of Medicine(2020-053). All patients signed informed consent forms before donating their specimens.

Author Contributions

GW and LY designed the study, analyzed data, and wrote the manuscript. XJ and RC designed the study and provided funding acquisition. YW, RH, KZ, and TG performed the experiments and analyzed the data. BC, XJ, and RC supervised the research and wrote the manuscript. All authors read and approved the final submitted manuscript.

Funding

Leading talent training program of Pudong New Area Health Committee (PWR12021-04); the Research project of Shanghai Municipal Health Commission (20204Y0302); Shanghai Science and Technology Committee (20Y11912100); the Research project of Pudong Science and Technology Commission (PKJ2020-Y20).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to thank the TCGA and GEO databases for the availability of the data. The authors would also like to thank professor Xiaomin Ying (Academy of Military Medical Sciences, Beijing, China) for providing the processed scRNA-seq data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.939836/full#supplementary-material

Supplementary Table 1 | The summarized clinical characteristics of GC cohorts used in this study.

Supplementary Table 2 | The baselines of the training and the internal validation sets.

Supplementary Table 3 | The sequence and detailed information of primers used in this research.

Supplementary Table 4 | The details of three sets of DEGs and overlapped IDGs.

Supplementary Table 5 | The results of GO and KEGG enrichment analyses on IDGs.

Supplementary Table 6 | The detailed results of univariate Cox regression analyses on IDGs.

Supplementary Table 7 | The selected five IDGs and their individual coefficients.

Supplementary Table 8 | | Hallmarks and KEGG pathways that significantly enriched in the low and high IPRS groups.

Abbreviations

ACRG, Asian Cancer Research Group; AJCC, American Joint Committee on Cancer; AUC, area under the curve; CI, confidence interval; CIBERSORT, cell type identification by estimating relative subsets of RNA transcripts; DEG, differentially expressed gene; DFS, disease-free survival; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumours using Expression data; FDR, false discovery rate; FPKM, Fragments Per Kilobase of exon model per Million mapped fragments; GC, gastric cancer; GDSC, Genomics of Drug Sensitivity in Cancer; GEO, Gene Expression Omnibus; GO, Gene Ontology; GSEA, gene set enrichment analysis; HR, hazard ratio; IC50, the half-maximal inhibitory concentration; ICB, immune checkpoint blocker; IDG, immune-related dysregulated gene; IPRS, immune-related prognostic risk score; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; MSI: Microsatellite Instable; MSigDB: Molecular Signatures Database; OS, overall survival; PCA, principal component analysis; ROC, receiver operating characteristic curve; RT–qPCR, real-time quantitative polymerase chain reaction; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; TIDE, Tumor Immune Dysfunction and Exclusion; TIME, tumor immune microenvironment; TIMER, Tumor Immune Estimation Resource; TMB, Tumor Mutation Burden; TPM, Transcripts Per kilobase of exon model per Million mapped reads.

References

1. Smyth EC, Nilsson M, Grabsch HI, van Grieken NC, Lordick F. Gastric Cancer. Lancet (2020) 396(10251):635–48. doi: 10.1016/S0140-6736(20)31288-5

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin (2021) 71(1):7–33. doi: 10.3322/caac.21654

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Johnston FM, Beckman M. Updates on Management of Gastric Cancer. Curr Oncol Rep (2019) 21(8):67. doi: 10.1007/s11912-019-0820-4

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Song Z, Wu Y, Yang J, Yang D, Fang X. Progress in the Treatment of Advanced Gastric Cancer. Tumour Biol (2017) 39(7):1010428317714626. doi: 10.1177/1010428317714626

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Biagioni A, Skalamera I, Peri S, Schiavone N, Cianchi F, Giommoni E, et al. Update on Gastric Cancer Treatments and Gene Therapies. Cancer Metastasis Rev (2019) 38(3):537–48. doi: 10.1007/s10555-019-09803-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Ribas A, Wolchok JD. Cancer Immunotherapy Using Checkpoint Blockade. Science (2018) 359(6382):1350–5. doi: 10.1126/science.aar4060

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Coutzac C, Pernot S, Chaput N, Zaanan A. Immunotherapy in Advanced Gastric Cancer, is it the Future? Crit Rev Oncol Hematol (2019) 133:25–32. doi: 10.1016/j.critrevonc.2018.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Wei L, Sun J, Zhang N, Zheng Y, Wang X, Lv L, et al. Noncoding RNAs in Gastric Cancer: Implications for Drug Resistance. Mol Cancer (2020) 19(1):62. doi: 10.1186/s12943-020-01185-7

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Choi SJ, Jung SW, Huh S, Chung YS, Cho H, Kang H, et al. Alteration of DNA Methylation in Gastric Cancer With Chemotherapy. J Microbiol Biotechnol (2017) 27(8):1367–78. doi: 10.4014/jmb.1704.04035

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Pansy K, Uhl B, Krstic J, Szmyra M, Fechter K, Santiso A, et al. Immune Regulatory Processes of the Tumor Microenvironment Under Malignant Conditions. Int J Mol Sci (2021) 22(24):13311. doi: 10.3390/ijms222413311

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Gajewski TF, Schreiber H, Fu YX. Innate and Adaptive Immune Cells in the Tumor Microenvironment. Nat Immunol (2013) 14(10):1014–22. doi: 10.1038/ni.2703

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Galon J, Bruni D. Approaches to Treat Immune Hot, Altered and Cold Tumours With Combination Immunotherapies. Nat Rev Drug Discov (2019) 18(3):197–218. doi: 10.1038/s41573-018-0007-y

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Farhood B, Najafi M, Mortezaee K. CD8(+) Cytotoxic T Lymphocytes in Cancer Immunotherapy: A Review. J Cell Physiol (2019) 234(6):8509–21. doi: 10.1002/jcp.27782

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Joshi SS, Badgwell BD. Current Treatment and Recent Progress in Gastric Cancer. CA Cancer J Clin (2021) 71(3):264–79. doi: 10.3322/caac.21657

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lou S, Wang Y, Zhang J, Yin X, Zhang Y, Wang Y, et al. Patient-Level DNA Damage Repair Pathway Profiles and Anti-Tumor Immunity for Gastric Cancer. Front Immunol (2021) 12:806324. doi: 10.3389/fimmu.2021.806324

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kim HK, Choi IJ, Kim CG, Kim HS, Oshima A, Yamada Y, et al. Three-Gene Predictor of Clinical Outcome for Gastric Cancer Patients Treated With Chemotherapy. Pharmacogenomics J (2012) 12(2):119–27. doi: 10.1038/tpj.2010.87

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Yoon SJ, Park J, Shin Y, Choi Y, Park SW, Kang SG, et al. Deconvolution of Diffuse Gastric Cancer and the Suppression of CD34 on the BALB/c Nude Mice Model. BMC Cancer (2020) 20(1):314. doi: 10.1186/s12885-020-06814-4

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Oh SC, Sohn BH, Cheong JH, Kim SB, Lee JE, Park KC, et al. Clinical and Genomic Landscape of Gastric Cancer With a Mesenchymal Phenotype. Nat Commun (2018) 9(1):1777. doi: 10.1038/s41467-018-04179-8

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. Tgfβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zhang M, Hu S, Min M, Ni Y, Lu Z, Sun X, et al. Dissecting Transcriptional Heterogeneity in Primary Gastric Adenocarcinoma by Single Cell RNA Sequencing. Gut (2021) 70(3):464–75. doi: 10.1136/gutjnl-2019-320368

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. Clusterprofiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. Innovation (N Y) (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141

CrossRef Full Text | Google Scholar

23. Ru B, Wong CN, Tong Y, Zhong JY, Zhong SSW, Wu WC, et al. TISIDB: An Integrated Repository Portal for Tumor-Immune System Interactions. Bioinformatics (2019) 35(20):4200–2. doi: 10.1093/bioinformatics/btz210

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Lánczky A, Győrffy B. Web-Based Survival Analysis Tool Tailored for Medical Research (KMplot): Development and Implementation. J Med Internet Res (2021) 23(7):e27633. doi: 10.2196/27633

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Maeser D, Gruener RF, Huang RS. Oncopredict: An R Package for Predicting In Vivo or Cancer Patient Drug Response and Biomarkers From Cell Line Screening Data. Brief Bioinform (2021) 22(6):bbab260. doi: 10.1093/bib/bbab260

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Oya Y, Hayakawa Y, Koike K. Tumor Microenvironment in Gastric Cancers. Cancer Sci (2020) 111(8):2696–707. doi: 10.1111/cas.14521

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Becerril-Rico J, Alvarado-Ortiz E, Toledo-Guzmán ME, Pelayo R, Ortiz-Sánchez E. The Cross Talk Between Gastric Cancer Stem Cells and the Immune Microenvironment: A Tumor-Promoting Factor. Stem Cell Res Ther (2021) 12(1):498. doi: 10.1186/s13287-021-02562-9

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bagaev A, Kotlov N, Nomie K, Svekolkin V, Gafurov A, Isaeva O, et al. Conserved Pan-Cancer Microenvironment Subtypes Predict Response to Immunotherapy. Cancer Cell (2021) 39(6):845–65.e7. doi: 10.1016/j.ccell.2021.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Collisson EA, Bailey P, Chang DK, Biankin AV. Molecular Subtypes of Pancreatic Cancer. Nat Rev Gastroenterol Hepatol (2019) 16(4):207–20. doi: 10.1038/s41575-019-0109-y

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Tan TZ, Rouanne M, Tan KT, Huang RY, Thiery JP. Molecular Subtypes of Urothelial Bladder Cancer: Results From a Meta-Cohort Analysis of 2411 Tumors. Eur Urol (2019) 75(3):423–32. doi: 10.1016/j.eururo.2018.08.027

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Chia NY, Tan P. Molecular Classification of Gastric Cancer. Ann Oncol (2016) 27(5):763–9. doi: 10.1093/annonc/mdw040

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Shao W, Yang Z, Fu Y, Zheng L, Liu F, Chai L, et al. The Pyroptosis-Related Signature Predicts Prognosis and Indicates Immune Microenvironment Infiltration in Gastric Cancer. Front Cell Dev Biol (2021) 9:676485. doi: 10.3389/fcell.2021.676485

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Han S, Zhu W, Yang W, Guan Q, Chen C, He Q, et al. A Prognostic Signature Constructed by CTHRC1 and LRFN4 in Stomach Adenocarcinoma. Front Genet (2021) 12:646818. doi: 10.3389/fgene.2021.646818

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Xia R, Tang H, Shen J, Xu S, Liang Y, Zhang Y, et al. Prognostic Value of a Novel Glycolysis-Related Gene Expression Signature for Gastrointestinal Cancer in the Asian Population. Cancer Cell Int (2021) 21(1):154. doi: 10.1186/s12935-021-01857-4

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Beauvais DM, Nelson SE, Adams KM, Stueven NA, Jung O, Rapraeger AC. Plasma Membrane Proteoglycans Syndecan-2 and Syndecan-4 Engage With EGFR and RON Kinase to Sustain Carcinoma Cell Cycle Progression. J Biol Chem (2022) 298(6):102029. doi: 10.1016/j.jbc.2022.102029

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Wang J, Wang Z, Lin W, Han Q, Yan H, Yao W, et al. LINC01296 Promotes Neuroblastoma Tumorigenesis via the NCL-SOX11 Regulatory Complex. Mol Ther Oncolytics (2022) 24:834–48. doi: 10.1016/j.omto.2022.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Hyun SY, Min HY, Lee HJ, Cho J, Boo HJ, Noh M, et al. Ninjurin1 Drives Lung Tumor Formation and Progression by Potentiating Wnt/β-Catenin Signaling Through Frizzled2-LRP6 Assembly. J Exp Clin Cancer Res (2022) 41(1):133. doi: 10.1186/s13046-022-02323-3

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Waks AG, Winer EP. Breast Cancer Treatment: A Review. Jama (2019) 321(3):288–300. doi: 10.1001/jama.2018.19323

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Cancer Genome Atlas Research Network. Comprehensive Molecular Characterization of Gastric Adenocarcinoma. Nature (2014) 513(7517):202–9. doi: 10.1038/nature13480

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Guinney J, Dienstmann R, Wang X, de Reyniès A, Schlicker A, Soneson C, et al. The Consensus Molecular Subtypes of Colorectal Cancer. Nat Med (2015) 21(11):1350–6. doi: 10.1038/nm.3967

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Röcken C. Molecular Classification of Gastric Cancer. Expert Rev Mol Diagn (2017) 17(3):293–301. doi: 10.1080/14737159.2017.1286985

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Xiao Z, Hu L, Yang L, Wang S, Gao Y, Zhu Q, et al. Tgfβ2 is a Prognostic-Related Biomarker and Correlated With Immune Infiltrates in Gastric Cancer. J Cell Mol Med (2020) 24(13):7151–62. doi: 10.1111/jcmm.15164

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Shi M, Gu Y, Jin K, Fang H, Chen Y, Cao Y, et al. CD47 Expression in Gastric Cancer Clinical Correlates and Association With Macrophage Infiltration. Cancer Immunol Immunother (2021) 70(7):1831–40. doi: 10.1007/s00262-020-02806-2

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang G, Yang L, Hu M, Hu R, Wang Y, Chen B, et al. Comprehensive Analysis of the Prognostic Significance of Hsa-miR-100-5p and Its Related Gene Signature in Stomach Adenocarcinoma. Front Cell Dev Biol (2021) 9:736274. doi: 10.3389/fcell.2021.736274

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Rihawi K, Ricci AD, Rizzo A, Brocchi S, Marasco G, Pastore LV, et al. Tumor-Associated Macrophages and Inflammatory Microenvironment in Gastric Cancer: Novel Translational Implications. Int J Mol Sci (2021) 22(8):3805. doi: 10.3390/ijms22083805

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Rojas A, Araya P, Gonzalez I, Morales E. Gastric Tumor Microenvironment. Adv Exp Med Biol (2020) 1226:23–35. doi: 10.1007/978-3-030-36214-0_2

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Joyce JA, Fearon DT. T Cell Exclusion, Immune Privilege, and the Tumor Microenvironment. Science (2015) 348(6230):74–80. doi: 10.1126/science.aaa6204

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Liu J, Mei J, Wang Y, Chen X, Pan J, Tong L, et al. Development of a Novel Immune-Related lncRNA Signature as a Prognostic Classifier for Endometrial Carcinoma. Int J Biol Sci (2021) 17(2):448–59. doi: 10.7150/ijbs.51207

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Gao Y, Liu J, Cai B, Chen Q, Wang G, Lu Z, et al. Development of Epithelial-Mesenchymal Transition-Related lncRNA Signature for Predicting Survival and Immune Microenvironment in Pancreatic Cancerwithexperiment Validation. Bioengineered (2021) 12(2):10553–67. doi: 10.1080/21655979.2021.2000197

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Yuan H, Liu J, Zhao L, Wu P, Chen G, Chen Q, et al. Prognostic Risk Model and Tumor Immune Environment Modulation of M5c-Related LncRNAs in Pancreatic Ductal Adenocarcinoma. Front Immunol (2021) 12:800268. doi: 10.3389/fimmu.2021.800268

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Liu J, Cui G, Ye J, Wang Y, Wang C, Bai J, et al. Comprehensive Analysis of the Prognostic Signature of Mutation-Derived Genome Instability-Related lncRNAs for Patients With Endometrial Cancer. Front Cell Dev Biol (2022) 10:753957. doi: 10.3389/fcell.2022.753957

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Engebretsen S, Bohlin J. Statistical Predictions With Glmnet. Clin Epigenet (2019) 11(1):123. doi: 10.1186/s13148-019-0730-1

CrossRef Full Text | Google Scholar

54. Garcia RI, Ibrahim JG, Zhu H. Variable Selection in the Cox Regression Model With Covariates Missing at Random. Biometrics (2010) 66(1):97–104. doi: 10.1111/j.1541-0420.2009.01274.x

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Jinawath N, Shiao MS, Norris A, Murphy K, Klein AP, Yonescu R, et al. Alterations of Type II Classical Cadherin, Cadherin-10 (CDH10), is Associated With Pancreatic Ductal Adenocarcinomas. Genes Chromosomes Cancer (2017) 56(5):427–35. doi: 10.1002/gcc.22447

PubMed Abstract | CrossRef Full Text | Google Scholar

56. An CH, Je EM, Yoo NJ, Lee SH. Frameshift Mutations of Cadherin Genes DCHS2, CDH10 and CDH24 Genes in Gastric and Colorectal Cancers With High Microsatellite Instability. Pathol Oncol Res (2015) 21(1):181–5. doi: 10.1007/s12253-014-9804-8

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Li C, Gao Z, Li F, Li X, Sun Y, Wang M, et al. Whole Exome Sequencing Identifies Frequent Somatic Mutations in Cell-Cell Adhesion Genes in Chinese Patients With Lung Squamous Cell Carcinoma. Sci Rep (2015) 5:14237. doi: 10.1038/srep14237

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Inculet RI, Peacock JL, Gorschboth CM, Norton JA. Gluconeogenesis in the Tumor-Influenced Rat Hepatocyte: Importance of Tumor Burden, Lactate, Insulin, and Glucagon. J Natl Cancer Inst (1987) 79(5):1039–46.

PubMed Abstract | Google Scholar

59. Wang Y, Zhu YD, Gui Q, Wang XD, Zhu YX. Glucagon-Induced Angiogenesis and Tumor Growth Through the HIF-1-VEGF-Dependent Pathway in Hyperglycemic Nude Mice. Genet Mol Res (2014) 13(3):7173–83. doi: 10.4238/2014.September.5.3

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Cuvelier A, Muir JF, Martin JP, Sesboüé R. [Proteins of the Inter-Alpha Trypsin Inhibitor (ITI) Family. A Major Role in the Biology of the Extracellular Matrix]. Rev Mal Respir (2000) 17(2):437–46.

PubMed Abstract | Google Scholar

61. Tripathi PH, Akhtar J, Arora J, Saran RK, Mishra N, Polisetty RV, et al. Quantitative Proteomic Analysis of GnRH Agonist Treated GBM Cell Line LN229 Revealed Regulatory Proteins Inhibiting Cancer Cell Proliferation. BMC Cancer (2022) 22(1):133. doi: 10.1186/s12885-022-09218-8

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Wang X, et al. CYTL1 Inhibits Tumor Metastasis With Decreasing STAT3 Phosphorylation. Oncoimmunology (2019) 8(5):e1577126. doi: 10.1080/2162402X.2019.1577126

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Wang X, Li T, Cheng Y, Wang P, Yuan W, Liu Q, et al. Intracellular CYTL1, a Novel Tumor Suppressor, Stabilizes NDUFV1 to Inhibit Metabolic Reprogramming in Breast Cancer. Signal Transduct Target Ther (2022) 7(1):35. doi: 10.1038/s41392-021-00856-1

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Zhu S, Kuek V, Bennett S, Xu H, Rosen V, Xu J. Protein Cytl1: Its Role in Chondrogenesis, Cartilage Homeostasis, and Disease. Cell Mol Life Sci (2019) 76(18):3515–23. doi: 10.1007/s00018-019-03137-x

PubMed Abstract | CrossRef Full Text | Google Scholar

65. James SR, Cedeno CD, Sharma A, Zhang W, Mohler JL, Odunsi K, et al. DNA Methylation and Nucleosome Occupancy Regulate the Cancer Germline Antigen Gene MAGEA11. Epigenetics (2013) 8(8):849–63. doi: 10.4161/epi.25500

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Hartmann S, Zwick L, Scheurer MJJ, Fuchs AR, Brands RC, Seher A, et al. MAGE-A11 Expression Contributes to Cisplatin Resistance in Head and Neck Cancer. Clin Oral Investig (2018) 22(3):1477–86. doi: 10.1007/s00784-017-2242-8

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Tokunaga R, Naseem M, Lo JH, Battaglin F, Soni S, Puccini A, et al. B Cell and B Cell-Related Pathways for Novel Cancer Treatments. Cancer Treat Rev (2019) 73:10–9. doi: 10.1016/j.ctrv.2018.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Kim J, Bae JS. Tumor-Associated Macrophages and Neutrophils in Tumor Microenvironment. Mediators Inflamm (2016) 2016:6058147. doi: 10.1155/2016/6058147

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Ngambenjawong C, Gustafson HH, Pun SH. Progress in Tumor-Associated Macrophage (TAM)-Targeted Therapeutics. Adv Drug Deliv Rev (2017) 114:206–21. doi: 10.1016/j.addr.2017.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Komi DEA, Redegeld FA. Role of Mast Cells in Shaping the Tumor Microenvironment. Clin Rev Allergy Immunol (2020) 58(3):313–25. doi: 10.1007/s12016-019-08753-w

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Caserta S, Borger JG, Zamoyska R. Central and Effector Memory CD4 and CD8 T-Cell Responses to Tumor-Associated Antigens. Crit Rev Immunol (2012) 32(2):97–126. doi: 10.1615/CritRevImmunol.v32.i2.10

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Hanahan D, Weinberg RA. Hallmarks of Cancer: The Next Generation. Cell (2011) 144(5):646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Viallard C, Larrivée B. Tumor Angiogenesis and Vascular Normalization: Alternative Therapeutic Targets. Angiogenesis (2017) 20(4):409–26. doi: 10.1007/s10456-017-9562-9

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Pastushenko I, Blanpain C. EMT Transition States During Tumor Progression and Metastasis. Trends Cell Biol (2019) 29(3):212–26. doi: 10.1016/j.tcb.2018.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Mittal V. Epithelial Mesenchymal Transition in Tumor Metastasis. Annu Rev Pathol (2018) 13:395–412. doi: 10.1146/annurev-pathol-020117-043854

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Hida K, Maishi N, Annan DA, Hida Y. Contribution of Tumor Endothelial Cells in Cancer Progression. Int J Mol Sci (2018) 19(5):1272. doi: 10.3390/ijms19051272

CrossRef Full Text | Google Scholar